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