1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkDataSetGradientPrecompute.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 "vtkDataSetGradientPrecompute.h"
22
23 #include "vtkMath.h"
24 #include "vtkTriangle.h"
25 #include "vtkInformationVector.h"
26 #include "vtkInformation.h"
27 #include "vtkDataSet.h"
28 #include "vtkDoubleArray.h"
29 #include "vtkCell.h"
30 #include "vtkCell3D.h"
31 #include "vtkTetra.h"
32 #include "vtkFieldData.h"
33 #include "vtkCellData.h"
34 #include "vtkObjectFactory.h"
35
36 #define VTK_DATASET_GRADIENT_TETRA_OPTIMIZATION
37 #define VTK_DATASET_GRADIENT_TRIANGLE_OPTIMIZATION
38 //#define DEBUG
39
40 vtkStandardNewMacro(vtkDataSetGradientPrecompute);
41
vtkDataSetGradientPrecompute()42 vtkDataSetGradientPrecompute::vtkDataSetGradientPrecompute()
43 {
44 }
45
~vtkDataSetGradientPrecompute()46 vtkDataSetGradientPrecompute::~vtkDataSetGradientPrecompute()
47 {
48 }
49
PrintSelf(ostream & os,vtkIndent indent)50 void vtkDataSetGradientPrecompute::PrintSelf(ostream& os, vtkIndent indent)
51 {
52 this->Superclass::PrintSelf(os,indent);
53 }
54
55 #define ADD_VEC(a,b) a[0]+=b[0];a[1]+=b[1];a[2]+=b[2]
56 #define SCALE_VEC(a,b) a[0]*=b;a[1]*=b;a[2]*=b
57 #define ZERO_VEC(a) a[0]=0;a[1]=0;a[2]=0
58 #define MAX_CELL_POINTS 128
59 #define VTK_CQS_EPSILON 1e-12
60
TETRA_CQS_VECTOR(double v0[3],double v1[3],double v2[3],double p[3],double cqs[3])61 static inline void TETRA_CQS_VECTOR( double v0[3], double v1[3], double v2[3], double p[3], double cqs[3] )
62 {
63 double surface = fabs( vtkTriangle::TriangleArea (v0,v1,v2) );
64
65 vtkTriangle::ComputeNormal( v0,v1,v2 , cqs );
66
67 // inverse face normal if not toward opposite vertex
68 double edge[3];
69 edge[0] = p[0] - v0[0];
70 edge[1] = p[1] - v0[1];
71 edge[2] = p[2] - v0[2];
72 if( vtkMath::Dot(edge,cqs) < 0 )
73 {
74 cqs[0] = - cqs[0];
75 cqs[1] = - cqs[1];
76 cqs[2] = - cqs[2];
77 }
78
79 SCALE_VEC( cqs , surface / 2.0 );
80 }
81
TRIANGLE_CQS_VECTOR(double v0[3],double v1[3],double p[3],double cqs[3])82 static inline void TRIANGLE_CQS_VECTOR( double v0[3], double v1[3], double p[3], double cqs[3] )
83 {
84 double length = sqrt(vtkMath::Distance2BetweenPoints(v0,v1));
85 double a[3], b[3], c[3];
86 for(int i=0;i<3;i++)
87 {
88 a[i] = v1[i] - v0[i];
89 b[i] = p[i] - v0[i];
90 }
91 vtkMath::Cross( a, b, c );
92 vtkMath::Cross( c , a , cqs );
93 vtkMath::Normalize(cqs);
94 SCALE_VEC( cqs , length / 2.0 );
95 }
96
LINE_CQS_VECTOR(double v0[3],double p[3],double cqs[3])97 static inline void LINE_CQS_VECTOR(double v0[3], double p[3], double cqs[3])
98 {
99 cqs[0] = p[0] - v0[0];
100 cqs[1] = p[1] - v0[1];
101 cqs[2] = p[2] - v0[2];
102 vtkMath::Normalize(cqs);
103 }
104
GradientPrecompute(vtkDataSet * ds)105 int vtkDataSetGradientPrecompute::GradientPrecompute(vtkDataSet* ds)
106 {
107 vtkIdType nCells = ds->GetNumberOfCells();
108 vtkIdType nCellNodes = 0;
109 for(vtkIdType i=0;i<nCells;i++)
110 {
111 nCellNodes += ds->GetCell(i)->GetNumberOfPoints();
112 }
113
114 vtkDoubleArray* cqs = vtkDoubleArray::New();
115 cqs->SetName("GradientPrecomputation");
116 cqs->SetNumberOfComponents(3);
117 cqs->SetNumberOfTuples(nCellNodes);
118 cqs->FillComponent(0, 0.0);
119 cqs->FillComponent(1, 0.0);
120 cqs->FillComponent(2, 0.0);
121
122 // The cell size determines the amount of space the cell takes up. For 3D
123 // cells this is the volume. For 2D cells this is the area. For 1D cells
124 // this is the length. For 0D cells this is undefined, but we set it to 1 so
125 // as not to get invalid results when normalizing something by the cell size.
126 vtkDoubleArray* cellSize = vtkDoubleArray::New();
127 cellSize->SetName("CellSize");
128 cellSize->SetNumberOfTuples(nCells);
129
130 vtkIdType curPoint = 0;
131 for(vtkIdType c=0;c<nCells;c++)
132 {
133 vtkCell* cell = ds->GetCell(c);
134 int np = cell->GetNumberOfPoints();
135
136 double cellCenter[3] = {0,0,0};
137 double cellPoints[MAX_CELL_POINTS][3];
138 double cellVectors[MAX_CELL_POINTS][3];
139 double tmp[3];
140 double size = 0.0;
141
142 for(int p=0;p<np;p++)
143 {
144 ds->GetPoint( cell->GetPointId(p), cellPoints[p] );
145 ADD_VEC( cellCenter , cellPoints[p] );
146 ZERO_VEC( cellVectors[p] );
147 }
148 SCALE_VEC(cellCenter,1.0/np);
149
150 // -= 3 D =-
151 if( cell->GetCellDimension() == 3 )
152 {
153 #ifdef VTK_DATASET_GRADIENT_TETRA_OPTIMIZATION
154 if( np == 4 ) // cell is a tetrahedra
155 {
156 //vtkWarningMacro(<<"Tetra detected\n");
157 size = fabs( vtkTetra::ComputeVolume(cellPoints[0], cellPoints[1], cellPoints[2], cellPoints[3]) ) *1.5 ;
158
159 TETRA_CQS_VECTOR( cellPoints[0], cellPoints[1], cellPoints[2], cellPoints[3] , tmp );
160 ADD_VEC(cellVectors[3],tmp);
161
162 TETRA_CQS_VECTOR( cellPoints[1], cellPoints[2], cellPoints[3], cellPoints[0] , tmp );
163 ADD_VEC(cellVectors[0],tmp);
164
165 TETRA_CQS_VECTOR( cellPoints[2], cellPoints[3], cellPoints[0], cellPoints[1] , tmp );
166 ADD_VEC(cellVectors[1],tmp);
167
168 TETRA_CQS_VECTOR( cellPoints[3], cellPoints[0], cellPoints[1], cellPoints[2] , tmp );
169 ADD_VEC(cellVectors[2],tmp);
170 }
171 else if( np > 4 )
172 #endif
173 {
174 vtkCell3D* cell3d = static_cast<vtkCell3D*>( cell );
175 int nf = cell->GetNumberOfFaces();
176 for(int f=0;f<nf;f++)
177 {
178 int* faceIds = 0;
179 int nfp = cell->GetFace(f)->GetNumberOfPoints();
180 cell3d->GetFacePoints(f,faceIds);
181 #ifdef VTK_DATASET_GRADIENT_TRIANGLE_OPTIMIZATION
182 if( nfp == 3 ) // face is a triangle
183 {
184 //vtkWarningMacro(<<"triangular face detected\n");
185 size+=fabs(vtkTetra::ComputeVolume(cellCenter,cellPoints[faceIds[0]],cellPoints[faceIds[1]],cellPoints[faceIds[2]]))*1.5;
186
187 TETRA_CQS_VECTOR( cellCenter, cellPoints[faceIds[0]], cellPoints[faceIds[1]], cellPoints[faceIds[2]] , tmp );
188 ADD_VEC(cellVectors[faceIds[2]],tmp);
189
190 TETRA_CQS_VECTOR( cellCenter, cellPoints[faceIds[1]], cellPoints[faceIds[2]], cellPoints[faceIds[0]] , tmp );
191 ADD_VEC(cellVectors[faceIds[0]],tmp);
192
193 TETRA_CQS_VECTOR( cellCenter, cellPoints[faceIds[2]], cellPoints[faceIds[0]], cellPoints[faceIds[1]] , tmp );
194 ADD_VEC(cellVectors[faceIds[1]],tmp);
195 }
196 else if( nfp > 3 ) // generic case
197 #endif
198 {
199 double faceCenter[3] = {0,0,0};
200 for(int p=0;p<nfp;p++)
201 {
202 ADD_VEC( faceCenter , cellPoints[faceIds[p]] );
203 }
204 SCALE_VEC( faceCenter, 1.0/nfp );
205 for(int p=0;p<nfp;p++)
206 {
207 int p2 = (p+1) % nfp ;
208 size += fabs( vtkTetra::ComputeVolume(cellCenter, faceCenter, cellPoints[faceIds[p]] , cellPoints[faceIds[p2]]) ) ;
209
210 TETRA_CQS_VECTOR( cellCenter, faceCenter, cellPoints[faceIds[p]] , cellPoints[faceIds[p2]] , tmp );
211 ADD_VEC( cellVectors[faceIds[p2]] , tmp );
212
213 TETRA_CQS_VECTOR( cellCenter, faceCenter, cellPoints[faceIds[p2]] , cellPoints[faceIds[p]] , tmp );
214 ADD_VEC( cellVectors[faceIds[p]] , tmp );
215 }
216 }
217 }
218 }
219 }
220
221 // -= 2 D =-
222 else if (cell->GetCellDimension() == 2)
223 {
224 if( np == 3 ) // cell is a triangle
225 {
226 size = fabs(vtkTriangle::TriangleArea(cellPoints[0], cellPoints[1], cellPoints[2]));
227
228 TRIANGLE_CQS_VECTOR( cellPoints[0] , cellPoints[1] , cellPoints[2] , tmp );
229 ADD_VEC( cellVectors[2] , tmp );
230
231 TRIANGLE_CQS_VECTOR( cellPoints[1] , cellPoints[2] , cellPoints[0] , tmp );
232 ADD_VEC( cellVectors[0] , tmp );
233
234 TRIANGLE_CQS_VECTOR( cellPoints[2] , cellPoints[0] , cellPoints[1] , tmp );
235 ADD_VEC( cellVectors[1] , tmp );
236 }
237 else if( np > 3) // generic case
238 {
239 for(int f=0;f<np;f++)
240 {
241 const int e0 = f;
242 const int e1 = (f+1)%np;
243 size += fabs(vtkTriangle::TriangleArea(cellCenter, cellPoints[e0], cellPoints[e1]));
244 TRIANGLE_CQS_VECTOR( cellCenter , cellPoints[e0] , cellPoints[e1] , tmp );
245 ADD_VEC( cellVectors[e1] , tmp );
246
247 TRIANGLE_CQS_VECTOR( cellCenter , cellPoints[e1] , cellPoints[e0] , tmp );
248 ADD_VEC( cellVectors[e0] , tmp );
249 }
250 }
251 else
252 {
253 //vtkWarningMacro(<<"Can't process 2D cells with less than 3 points.");
254 //return 0;
255 }
256 }
257
258 // -= 1 D =-
259 else if (cell->GetCellDimension() == 1)
260 {
261 if (np == 2) // cell is a single line segment
262 {
263 size
264 = sqrt(vtkMath::Distance2BetweenPoints(cellPoints[0], cellPoints[1]));
265
266 LINE_CQS_VECTOR(cellPoints[0], cellPoints[1], tmp);
267 ADD_VEC(cellVectors[1], tmp);
268
269 LINE_CQS_VECTOR(cellPoints[1], cellPoints[0], tmp);
270 ADD_VEC(cellVectors[0], tmp);
271 }
272 else if (np > 2) // generic case, a poly line
273 {
274 for (int p = 0; p < np; p++)
275 {
276 size
277 += sqrt(vtkMath::Distance2BetweenPoints(cellCenter, cellPoints[p]));
278 LINE_CQS_VECTOR(cellCenter, cellPoints[p], tmp);
279 ADD_VEC(cellVectors[p], tmp);
280 }
281 }
282 }
283
284 // -= 0 D =-
285 else
286 {
287 // For vertex cells, estimate gradient as weighted sum of vectors from
288 // centroid.
289 size = 1.0;
290 for (int p = 0; p < np; p++)
291 {
292 cellVectors[p][0] = cellPoints[p][0] - cellCenter[0];
293 cellVectors[p][1] = cellPoints[p][1] - cellCenter[1];
294 cellVectors[p][2] = cellPoints[p][2] - cellCenter[2];
295 }
296 }
297
298 cellSize->SetTuple1(c,size);
299
300 // check cqs consistency
301 #ifdef DEBUG
302 double checkZero[3] = {0,0,0};
303 double checkVolume = 0;
304 for(int p=0;p<np;p++)
305 {
306 checkVolume += vtkMath::Dot( cellPoints[p] , cellVectors[p] );
307 ADD_VEC(checkZero,cellVectors[p]);
308 cqs->SetTuple( curPoint + p , cellVectors[p] );
309 }
310 checkVolume /= (double) cell->GetCellDimension();
311
312 if( vtkMath::Norm(checkZero)>VTK_CQS_EPSILON || fabs(size-checkVolume)>VTK_CQS_EPSILON )
313 {
314 cout<<"Bad CQS sum at cell #"<<c<<", Sum="<<vtkMath::Norm(checkZero)<<", volume="<<size<<", ratio Vol="<<size/checkVolume<<"\n";
315 }
316 #endif
317
318 curPoint += np;
319 }
320
321 ds->GetFieldData()->AddArray( cqs );
322 ds->GetCellData()->AddArray( cellSize );
323 cqs->Delete();
324 cellSize->Delete();
325
326 return 1;
327 }
328
RequestData(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)329 int vtkDataSetGradientPrecompute::RequestData(vtkInformation *vtkNotUsed(request),
330 vtkInformationVector **inputVector,
331 vtkInformationVector *outputVector)
332 {
333 // get the info objects
334 vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
335 vtkInformation *outInfo = outputVector->GetInformationObject(0);
336
337 // get connected input & output
338 vtkDataSet* _output = vtkDataSet::SafeDownCast( outInfo->Get(vtkDataObject::DATA_OBJECT()) );
339 vtkDataSet* _input = vtkDataSet::SafeDownCast( inInfo->Get(vtkDataObject::DATA_OBJECT()) );
340
341 if( _input==0 || _output==0 )
342 {
343 vtkErrorMacro(<<"missing input/output connection\n");
344 return 0;
345 }
346
347 _output->ShallowCopy(_input);
348 return vtkDataSetGradientPrecompute::GradientPrecompute(_output);
349 }
350
351