1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkDataSetGradient.cxx
5
6 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7 All rights reserved.
8 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9
10 This software is distributed WITHOUT ANY WARRANTY; without even
11 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12 PURPOSE. See the above copyright notice for more information.
13
14 =========================================================================*/
15 // .SECTION Thanks
16 // This file is part of the generalized Youngs material interface reconstruction algorithm contributed by
17 // CEA/DIF - Commissariat a l'Energie Atomique, Centre DAM Ile-De-France <br>
18 // BP12, F-91297 Arpajon, France. <br>
19 // Implementation by Thierry Carrard (CEA)
20
21 #include "vtkDataSetGradient.h"
22 #include "vtkDataSetGradientPrecompute.h"
23 #include "vtkDataSet.h"
24 #include "vtkDataArray.h"
25 #include "vtkDoubleArray.h"
26 #include "vtkPointData.h"
27 #include "vtkCellData.h"
28 #include "vtkCell.h"
29 #include "vtkMath.h"
30 #include "vtkInformation.h"
31 #include "vtkInformationVector.h"
32 #include "vtkObjectFactory.h"
33
34
35 // utility macros
36 #define ADD_VEC(a,b) a[0]+=b[0];a[1]+=b[1];a[2]+=b[2]
37 #define SCALE_VEC(a,b) a[0]*=b;a[1]*=b;a[2]*=b
38 #define ZERO_VEC(a) a[0]=0;a[1]=0;a[2]=0
39 #define COPY_VEC(a,b) a[0]=b[0];a[1]=b[1];a[2]=b[2];
40 #define MAX_CELL_POINTS 128
41 #define VTK_CQS_EPSILON 1e-12
42
43
44 // standard constructors and factory
45 vtkStandardNewMacro(vtkDataSetGradient);
46
47 /*!
48 The default constructor
49 \sa ~vtkDataSetGradient()
50 */
vtkDataSetGradient()51 vtkDataSetGradient::vtkDataSetGradient()
52 : ResultArrayName(0)
53 {
54 this->SetResultArrayName("gradient");
55 }
56
57 /*!
58 The destrcutor
59 \sa vtkDataSetGradient()
60 */
~vtkDataSetGradient()61 vtkDataSetGradient::~vtkDataSetGradient()
62 {
63 this->SetResultArrayName( 0 );
64 }
65
PrintSelf(ostream & os,vtkIndent indent)66 void vtkDataSetGradient::PrintSelf(ostream& os, vtkIndent indent)
67 {
68 this->Superclass::PrintSelf(os,indent);
69 os << indent << "Result array name: " << this->ResultArrayName << "\n";
70 }
71
RequestData(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)72 int vtkDataSetGradient::RequestData(vtkInformation * vtkNotUsed(request),
73 vtkInformationVector **inputVector,
74 vtkInformationVector *outputVector)
75 {
76 // get the info objects
77 vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
78 vtkInformation *outInfo = outputVector->GetInformationObject(0);
79
80 // get connected input & output
81 vtkDataSet *_output = vtkDataSet::SafeDownCast( outInfo->Get(vtkDataObject::DATA_OBJECT()) );
82 vtkDataSet* _input = vtkDataSet::SafeDownCast( inInfo->Get(vtkDataObject::DATA_OBJECT()) );
83
84 if( _input==0 || _output==0 )
85 {
86 vtkErrorMacro(<<"Missing input or output \n");
87 return 0;
88 }
89
90 // get array to compute gradient from
91 vtkDataArray* inArray = this->GetInputArrayToProcess( 0, _input );
92 if( inArray==0 )
93 {
94 inArray = _input->GetPointData()->GetScalars();
95 }
96 if( inArray==0 )
97 {
98 inArray = _input->GetCellData()->GetScalars();
99 }
100
101 if( inArray==0 )
102 {
103 vtkErrorMacro(<<"no input array to process\n");
104 return 0;
105 }
106
107 vtkDebugMacro(<<"Input array to process : "<<inArray->GetName()<<"\n");
108
109 bool pointData;
110 if( _input->GetCellData()->GetArray(inArray->GetName()) == inArray )
111 {
112 pointData = false;
113 vtkDebugMacro(<<"cell data to point gradient\n");
114 }
115 else if( _input->GetPointData()->GetArray(inArray->GetName()) == inArray )
116 {
117 pointData = true;
118 vtkDebugMacro(<<"point data to cell gradient\n");
119 }
120 else
121 {
122 vtkErrorMacro(<<"input array must be cell or point data\n");
123 return 0;
124 }
125
126 // we're just adding a scalar field
127 _output->ShallowCopy( _input );
128
129 vtkDataArray* cqsArray = _output->GetFieldData()->GetArray("GradientPrecomputation");
130 vtkDataArray* sizeArray = _output->GetCellData()->GetArray("CellSize");
131 if( cqsArray==0 || sizeArray==0 )
132 {
133 vtkDebugMacro(<<"Couldn't find field array 'GradientPrecomputation', computing it right now.\n");
134 vtkDataSetGradientPrecompute::GradientPrecompute(_output);
135 cqsArray = _output->GetFieldData()->GetArray("GradientPrecomputation");
136 sizeArray = _output->GetCellData()->GetArray("CellSize");
137 if( cqsArray==0 || sizeArray==0 )
138 {
139 vtkErrorMacro(<<"Computation of field array 'GradientPrecomputation' or 'CellSize' failed.\n");
140 return 0;
141 }
142 }
143
144 vtkIdType nCells = _input->GetNumberOfCells();
145 vtkIdType nPoints = _input->GetNumberOfPoints();
146
147 vtkDoubleArray* gradientArray = vtkDoubleArray::New();
148 gradientArray->SetName( this->ResultArrayName );
149 gradientArray->SetNumberOfComponents(3);
150
151 if( pointData ) // compute cell gradient from point data
152 {
153 gradientArray->SetNumberOfTuples( nCells );
154 vtkIdType cellPoint = 0;
155 for(vtkIdType i=0;i<nCells;i++)
156 {
157 vtkCell* cell = _input->GetCell(i);
158 int np = cell->GetNumberOfPoints();
159 double gradient[3] = {0,0,0};
160 for(int p=0;p<np;p++)
161 {
162 double cqs[3], scalar;
163 cqsArray->GetTuple( cellPoint++ , cqs );
164 scalar = inArray->GetTuple1( cell->GetPointId(p) );
165 SCALE_VEC( cqs , scalar );
166 ADD_VEC( gradient , cqs );
167 }
168 SCALE_VEC( gradient , ( 1.0 / sizeArray->GetTuple1(i) ) );
169 gradientArray->SetTuple( i , gradient );
170 }
171
172 _output->GetCellData()->AddArray( gradientArray );
173 //_output->GetCellData()->SetVectors( gradientArray );
174 }
175 else // compute point gradient from cell data
176 {
177 gradientArray->SetNumberOfTuples( nPoints );
178 gradientArray->FillComponent(0, 0.0);
179 gradientArray->FillComponent(1, 0.0);
180 gradientArray->FillComponent(2, 0.0);
181 double * gradient = gradientArray->WritePointer(0,nPoints*3);
182 double * gradientDivisor = new double [nPoints];
183 for(vtkIdType i=0;i<nPoints;i++)
184 {
185 gradientDivisor[i] = 0.0;
186 }
187 vtkIdType cellPoint = 0;
188 for(vtkIdType i=0;i<nCells;i++)
189 {
190 vtkCell* cell = _input->GetCell(i);
191 int np = cell->GetNumberOfPoints();
192 double scalar = inArray->GetTuple1( i );
193 for(int p=0;p<np;p++)
194 {
195 double cqs[3];
196 double pointCoord[3];
197 vtkIdType pointId = cell->GetPointId(p);
198 cqsArray->GetTuple( cellPoint++ , cqs );
199 _input->GetPoint( cell->GetPointId(p), pointCoord );
200 scalar *= cell->GetCellDimension();
201 SCALE_VEC( (gradient+pointId*3) , scalar );
202 gradientDivisor[pointId] += vtkMath::Dot(cqs,pointCoord);
203 }
204 }
205 for(vtkIdType i=0;i<nPoints;i++)
206 {
207 SCALE_VEC( (gradient+i*3) , (1.0/gradientDivisor[i]) );
208 }
209 delete [] gradientDivisor;
210 _output->GetPointData()->AddArray( gradientArray );
211 //_output->GetPointData()->SetVectors( gradientArray );
212 }
213 gradientArray->Delete();
214
215 vtkDebugMacro(<<_output->GetClassName()<<" @ "<<_output<<" :\n");
216
217 return 1;
218 }
219
220