1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkBiQuadraticQuad.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
16 //Thanks to Soeren Gebbert who developed this class and
17 //integrated it into VTK 5.0.
18
19 #include "vtkBiQuadraticQuad.h"
20
21 #include "vtkObjectFactory.h"
22 #include "vtkDoubleArray.h"
23 #include "vtkMath.h"
24 #include "vtkPointData.h"
25 #include "vtkQuad.h"
26 #include "vtkQuadraticEdge.h"
27 #include "vtkPoints.h"
28
29 vtkStandardNewMacro(vtkBiQuadraticQuad);
30
31 //----------------------------------------------------------------------------
32 // Construct the quad with nine points.
vtkBiQuadraticQuad()33 vtkBiQuadraticQuad::vtkBiQuadraticQuad()
34 {
35 this->Edge = vtkQuadraticEdge::New();
36 this->Quad = vtkQuad::New();
37 this->Points->SetNumberOfPoints(9);
38 this->PointIds->SetNumberOfIds(9);
39 for (int i = 0; i < 9; i++)
40 {
41 this->Points->SetPoint(i, 0.0, 0.0, 0.0);
42 this->PointIds->SetId(i,0);
43 }
44 this->Scalars = vtkDoubleArray::New();
45 this->Scalars->SetNumberOfTuples(4);
46 }
47
48 //----------------------------------------------------------------------------
~vtkBiQuadraticQuad()49 vtkBiQuadraticQuad::~vtkBiQuadraticQuad()
50 {
51 this->Edge->Delete();
52 this->Quad->Delete();
53
54 this->Scalars->Delete();
55 }
56
57 //----------------------------------------------------------------------------
GetEdge(int edgeId)58 vtkCell *vtkBiQuadraticQuad::GetEdge(int edgeId)
59 {
60 edgeId = (edgeId < 0 ? 0 : (edgeId > 3 ? 3 : edgeId));
61 int p = (edgeId + 1) % 4;
62
63 // load point id's
64 this->Edge->PointIds->SetId (0, this->PointIds->GetId(edgeId));
65 this->Edge->PointIds->SetId (1, this->PointIds->GetId(p));
66 this->Edge->PointIds->SetId (2, this->PointIds->GetId(edgeId + 4));
67
68 // load coordinates
69 this->Edge->Points->SetPoint (0, this->Points->GetPoint(edgeId));
70 this->Edge->Points->SetPoint (1, this->Points->GetPoint(p));
71 this->Edge->Points->SetPoint (2, this->Points->GetPoint(edgeId + 4));
72
73 return this->Edge;
74 }
75
76 //----------------------------------------------------------------------------
77 static int LinearQuads[4][4] = { {0, 4, 8, 7}, {4, 1, 5, 8},
78 {8, 5, 2, 6}, {7, 8, 6, 3} };
79
80 //----------------------------------------------------------------------------
EvaluatePosition(const double x[3],double * closestPoint,int & subId,double pcoords[3],double & minDist2,double * weights)81 int vtkBiQuadraticQuad::EvaluatePosition (const double x[3],
82 double *closestPoint,
83 int &subId, double pcoords[3],
84 double &minDist2, double *weights)
85 {
86 double pc[3], dist2;
87 int ignoreId, i, returnStatus = 0, status;
88 double tempWeights[4];
89 double closest[3];
90
91 //four linear quads are used
92 for (minDist2 = VTK_DOUBLE_MAX, i = 0; i < 4; i++)
93 {
94 this->Quad->Points->SetPoint (0,
95 this->Points->GetPoint (LinearQuads[i][0]));
96 this->Quad->Points->SetPoint (1,
97 this->Points->GetPoint (LinearQuads[i][1]));
98 this->Quad->Points->SetPoint (2,
99 this->Points->GetPoint (LinearQuads[i][2]));
100 this->Quad->Points->SetPoint (3,
101 this->Points->GetPoint (LinearQuads[i][3]));
102
103 status = this->Quad->EvaluatePosition (x, closest, ignoreId, pc, dist2,
104 tempWeights);
105 if (status != -1 && dist2 < minDist2)
106 {
107 returnStatus = status;
108 minDist2 = dist2;
109 subId = i;
110 pcoords[0] = pc[0];
111 pcoords[1] = pc[1];
112 }
113 }
114
115 // adjust parametric coordinates
116 if (returnStatus != -1)
117 {
118 if (subId == 0)
119 {
120 pcoords[0] /= 2.0;
121 pcoords[1] /= 2.0;
122 }
123 else if (subId == 1)
124 {
125 pcoords[0] = 0.5 + (pcoords[0] / 2.0);
126 pcoords[1] /= 2.0;
127 }
128 else if (subId == 2)
129 {
130 pcoords[0] = 0.5 + (pcoords[0] / 2.0);
131 pcoords[1] = 0.5 + (pcoords[1] / 2.0);
132 }
133 else
134 {
135 pcoords[0] /= 2.0;
136 pcoords[1] = 0.5 + (pcoords[1] / 2.0);
137 }
138 pcoords[2] = 0.0;
139 if(closestPoint!=nullptr)
140 {
141 // Compute both closestPoint and weights
142 this->EvaluateLocation(subId,pcoords,closestPoint,weights);
143 }
144 else
145 {
146 // Compute weigths only
147 this->InterpolationFunctionsPrivate(pcoords,weights);
148 }
149 }
150
151 return returnStatus;
152 }
153
154 //----------------------------------------------------------------------------
EvaluateLocation(int & vtkNotUsed (subId),const double pcoords[3],double x[3],double * weights)155 void vtkBiQuadraticQuad::EvaluateLocation (int& vtkNotUsed(subId),
156 const double pcoords[3],
157 double x[3], double *weights)
158 {
159 int i, j;
160 double pt[3];
161
162 this->InterpolationFunctionsPrivate(pcoords,weights);
163
164 x[0] = x[1] = x[2] = 0.0;
165 for (i=0; i<9; i++)
166 {
167 this->Points->GetPoint(i, pt);
168 for (j=0; j<3; j++)
169 {
170 x[j] += pt[j] * weights[i];
171 }
172 }
173 }
174
175 //----------------------------------------------------------------------------
CellBoundary(int subId,const double pcoords[3],vtkIdList * pts)176 int vtkBiQuadraticQuad::CellBoundary(int subId, const double pcoords[3], vtkIdList * pts)
177 {
178 return this->Quad->CellBoundary (subId, pcoords, pts);
179 }
180
181
182 //----------------------------------------------------------------------------
183 void
Contour(double value,vtkDataArray * cellScalars,vtkIncrementalPointLocator * locator,vtkCellArray * verts,vtkCellArray * lines,vtkCellArray * polys,vtkPointData * inPd,vtkPointData * outPd,vtkCellData * inCd,vtkIdType cellId,vtkCellData * outCd)184 vtkBiQuadraticQuad::Contour (double value,
185 vtkDataArray *cellScalars,
186 vtkIncrementalPointLocator * locator,
187 vtkCellArray * verts,
188 vtkCellArray * lines,
189 vtkCellArray * polys,
190 vtkPointData * inPd,
191 vtkPointData * outPd, vtkCellData * inCd, vtkIdType cellId, vtkCellData * outCd)
192 {
193 //contour each linear quad separately
194 for (int i=0; i<4; i++)
195 {
196 for (int j=0; j<4; j++)
197 {
198 this->Quad->Points->SetPoint(j,this->Points->GetPoint(LinearQuads[i][j]));
199 this->Quad->PointIds->SetId(j,this->PointIds->GetId(LinearQuads[i][j]));
200 this->Scalars->SetValue(j,cellScalars->GetTuple1(LinearQuads[i][j]));
201 }
202
203 this->Quad->Contour(value,this->Scalars,locator,verts,lines,polys,
204 inPd,outPd,inCd,cellId,outCd);
205 }
206 }
207
208 //----------------------------------------------------------------------------
209 // Clip this quadratic quad using scalar value provided. Like contouring,
210 // except that it cuts the quad to produce other quads and triangles.
211 void
Clip(double value,vtkDataArray * cellScalars,vtkIncrementalPointLocator * locator,vtkCellArray * polys,vtkPointData * inPd,vtkPointData * outPd,vtkCellData * inCd,vtkIdType cellId,vtkCellData * outCd,int insideOut)212 vtkBiQuadraticQuad::Clip (double value, vtkDataArray * cellScalars,
213 vtkIncrementalPointLocator * locator, vtkCellArray * polys,
214 vtkPointData * inPd, vtkPointData * outPd,
215 vtkCellData * inCd, vtkIdType cellId, vtkCellData * outCd, int insideOut)
216 {
217 //contour each linear quad separately
218 for (int i=0; i<4; i++)
219 {
220 for ( int j=0; j<4; j++) //for each of the four vertices of the linear quad
221 {
222 this->Quad->Points->SetPoint(j,this->Points->GetPoint(LinearQuads[i][j]));
223 this->Quad->PointIds->SetId(j,this->PointIds->GetId(LinearQuads[i][j]));
224 this->Scalars->SetValue(j,cellScalars->GetTuple1(LinearQuads[i][j]));
225 }
226
227 this->Quad->Clip(value,this->Scalars,locator,polys,inPd,
228 outPd,inCd,cellId,outCd,insideOut);
229 }
230 }
231
232 //----------------------------------------------------------------------------
233 // Line-line intersection. Intersection has to occur within [0,1] parametric
234 // coordinates and with specified tolerance.
235 int
IntersectWithLine(const double * p1,const double * p2,double tol,double & t,double * x,double * pcoords,int & subId)236 vtkBiQuadraticQuad::IntersectWithLine (const double *p1,
237 const double *p2, double tol, double &t, double *x, double *pcoords, int &subId)
238 {
239 int subTest, i;
240 subId = 0;
241
242 //intersect the four linear quads
243 for (i = 0; i < 4; i++)
244 {
245 this->Quad->Points->SetPoint (0,
246 this->Points->GetPoint(LinearQuads[i][0]));
247 this->Quad->Points->SetPoint (1,
248 this->Points->GetPoint(LinearQuads[i][1]));
249 this->Quad->Points->SetPoint (2,
250 this->Points->GetPoint(LinearQuads[i][2]));
251 this->Quad->Points->SetPoint (3,
252 this->Points->GetPoint(LinearQuads[i][3]));
253
254 if (this->Quad->IntersectWithLine (p1, p2, tol, t, x, pcoords, subTest))
255 {
256 return 1;
257 }
258 }
259
260 return 0;
261 }
262
263 //----------------------------------------------------------------------------
Triangulate(int vtkNotUsed (index),vtkIdList * ptIds,vtkPoints * pts)264 int vtkBiQuadraticQuad::Triangulate (int vtkNotUsed(index),
265 vtkIdList *ptIds,
266 vtkPoints *pts)
267 {
268 pts->SetNumberOfPoints(24);
269 ptIds->SetNumberOfIds(24);
270
271 // First the corner vertices
272 ptIds->SetId(0,this->PointIds->GetId(0));
273 ptIds->SetId(1,this->PointIds->GetId(4));
274 ptIds->SetId(2,this->PointIds->GetId(7));
275 pts->SetPoint(0,this->Points->GetPoint(0));
276 pts->SetPoint(1,this->Points->GetPoint(4));
277 pts->SetPoint(2,this->Points->GetPoint(7));
278
279 ptIds->SetId(3,this->PointIds->GetId(4));
280 ptIds->SetId(4,this->PointIds->GetId(1));
281 ptIds->SetId(5,this->PointIds->GetId(5));
282 pts->SetPoint(3,this->Points->GetPoint(4));
283 pts->SetPoint(4,this->Points->GetPoint(1));
284 pts->SetPoint(5,this->Points->GetPoint(5));
285
286 ptIds->SetId(6,this->PointIds->GetId(5));
287 ptIds->SetId(7,this->PointIds->GetId(2));
288 ptIds->SetId(8,this->PointIds->GetId(6));
289 pts->SetPoint(6,this->Points->GetPoint(5));
290 pts->SetPoint(7,this->Points->GetPoint(2));
291 pts->SetPoint(8,this->Points->GetPoint(6));
292
293 ptIds->SetId(9,this->PointIds->GetId(6));
294 ptIds->SetId(10,this->PointIds->GetId(3));
295 ptIds->SetId(11,this->PointIds->GetId(7));
296 pts->SetPoint(9,this->Points->GetPoint(6));
297 pts->SetPoint(10,this->Points->GetPoint(3));
298 pts->SetPoint(11,this->Points->GetPoint(7));
299
300 //Now the triangles in the middle
301 ptIds->SetId(12,this->PointIds->GetId(4));
302 ptIds->SetId(13,this->PointIds->GetId(8));
303 ptIds->SetId(14,this->PointIds->GetId(7));
304 pts->SetPoint(12,this->Points->GetPoint(4));
305 pts->SetPoint(13,this->Points->GetPoint(8));
306 pts->SetPoint(14,this->Points->GetPoint(7));
307
308 ptIds->SetId(15,this->PointIds->GetId(4));
309 ptIds->SetId(16,this->PointIds->GetId(5));
310 ptIds->SetId(17,this->PointIds->GetId(8));
311 pts->SetPoint(15,this->Points->GetPoint(4));
312 pts->SetPoint(16,this->Points->GetPoint(5));
313 pts->SetPoint(17,this->Points->GetPoint(8));
314
315 ptIds->SetId(18,this->PointIds->GetId(5));
316 ptIds->SetId(19,this->PointIds->GetId(6));
317 ptIds->SetId(20,this->PointIds->GetId(8));
318 pts->SetPoint(18,this->Points->GetPoint(5));
319 pts->SetPoint(19,this->Points->GetPoint(6));
320 pts->SetPoint(20,this->Points->GetPoint(8));
321
322 ptIds->SetId(21,this->PointIds->GetId(6));
323 ptIds->SetId(22,this->PointIds->GetId(7));
324 ptIds->SetId(23,this->PointIds->GetId(8));
325 pts->SetPoint(21,this->Points->GetPoint(6));
326 pts->SetPoint(22,this->Points->GetPoint(7));
327 pts->SetPoint(23,this->Points->GetPoint(8));
328
329 return 1;
330 }
331
332 //----------------------------------------------------------------------------
333 void
Derivatives(int vtkNotUsed (subId),const double pcoords[3],const double * values,int dim,double * derivs)334 vtkBiQuadraticQuad::Derivatives (int vtkNotUsed (subId),
335 const double pcoords[3], const double *values,
336 int dim, double *derivs)
337 {
338 double sum[3], weights[9];
339 double functionDerivs[18];
340 double elemNodes[9][3];
341 double *J[3], J0[3], J1[3], J2[3];
342 double *JI[3], JI0[3], JI1[3], JI2[3];
343
344 for(int i = 0; i<9; i++)
345 {
346 this->Points->GetPoint(i, elemNodes[i]);
347 }
348
349 this->InterpolationFunctionsPrivate(pcoords,weights);
350 this->InterpolationDerivsPrivate(pcoords,functionDerivs);
351
352 // Compute transposed Jacobian and inverse Jacobian
353 J[0] = J0; J[1] = J1; J[2] = J2;
354 JI[0] = JI0; JI[1] = JI1; JI[2] = JI2;
355 for(int k = 0; k<3; k++)
356 {
357 J0[k] = J1[k] = 0.0;
358 }
359
360 for(int i = 0; i<9; i++)
361 {
362 for(int j = 0; j<2; j++)
363 {
364 for(int k = 0; k<3; k++)
365 {
366 J[j][k] += elemNodes[i][k] * functionDerivs[j*9+i];
367 }
368 }
369 }
370
371 // Compute third row vector in transposed Jacobian and normalize it, so that Jacobian determinant stays the same.
372 vtkMath::Cross(J0,J1,J2);
373 if ( vtkMath::Normalize(J2) == 0.0 || !vtkMath::InvertMatrix(J,JI,3)) //degenerate
374 {
375 for (int j=0; j < dim; j++ )
376 {
377 for (int i=0; i < 3; i++ )
378 {
379 derivs[j*dim + i] = 0.0;
380 }
381 }
382 return;
383 }
384
385
386 // Loop over "dim" derivative values. For each set of values,
387 // compute derivatives
388 // in local system and then transform into modelling system.
389 // First compute derivatives in local x'-y' coordinate system
390 for (int j=0; j < dim; j++ )
391 {
392 sum[0] = sum[1] = sum[2] = 0.0;
393 for (int i=0; i < 9; i++) //loop over interp. function derivatives
394 {
395 sum[0] += functionDerivs[i] * values[dim*i + j];
396 sum[1] += functionDerivs[9 + i] * values[dim*i + j];
397 }
398 // dBydx = sum[0]*JI[0][0] + sum[1]*JI[0][1];
399 // dBydy = sum[0]*JI[1][0] + sum[1]*JI[1][1];
400
401 // Transform into global system (dot product with global axes)
402 derivs[3*j] = sum[0]*JI[0][0] + sum[1]*JI[0][1];
403 derivs[3*j + 1] = sum[0]*JI[1][0] + sum[1]*JI[1][1];
404 derivs[3*j + 2] = sum[0]*JI[2][0] + sum[1]*JI[2][1];
405 }
406 }
407
408
409
410 //----------------------------------------------------------------------------
411 // Compute interpolation functions. The first four nodes are the corner
412 // vertices; the others are mid-edge nodes, the last one is the mid-center
413 // node.
InterpolationFunctionsPrivate(const double pcoords[3],double weights[9])414 void vtkBiQuadraticQuad::InterpolationFunctionsPrivate (const double pcoords[3],
415 double weights[9])
416 {
417 //Normally these coordinates are named r and s, but I chose x and y,
418 //because you can easily mark and paste these functions to the
419 //gnuplot splot function. :D
420 double x = pcoords[0];
421 double y = pcoords[1];
422
423 //midedge weights
424 weights[0] = 4.0 * (1.0 - x) * (x - 0.5) * (1.0 - y) * (y - 0.5);
425 weights[1] = -4.0 * (x) * (x - 0.5) * (1.0 - y) * (y - 0.5);
426 weights[2] = 4.0 * (x) * (x - 0.5) * (y) * (y - 0.5);
427 weights[3] = -4.0 * (1.0 - x) * (x - 0.5) * (y) * (y - 0.5);
428 //corner weights
429 weights[4] = 8.0 * (x) * (1.0 - x) * (1.0 - y) * (0.5 - y);
430 weights[5] = -8.0 * (x) * (0.5 - x) * (1.0 - y) * (y);
431 weights[6] = -8.0 * (x) * (1.0 - x) * (y) * (0.5 - y);
432 weights[7] = 8.0 * (1.0 - x) * (0.5 - x) * (1.0 - y) * (y);
433 //surface center weight
434 weights[8] = 16.0 * (x) * (1.0 - x) * (1.0 - y) * (y);
435 }
436
437 //----------------------------------------------------------------------------
438 // Derivatives in parametric space.
InterpolationDerivsPrivate(const double pcoords[3],double derivs[18])439 void vtkBiQuadraticQuad::InterpolationDerivsPrivate (const double pcoords[3],
440 double derivs[18])
441 {
442 // Coordinate conversion
443 double x = pcoords[0];
444 double y = pcoords[1];
445
446 // Derivatives in the x-direction
447 // edge
448 derivs[0] = 4.0 * (1.5 - 2.0 * x) * (1.0 - y) * (y - 0.5);
449 derivs[1] =-4.0 * (2.0 * x - 0.5) * (1.0 - y) * (y - 0.5);
450 derivs[2] = 4.0 * (2.0 * x - 0.5) * (y) * (y - 0.5);
451 derivs[3] =-4.0 * (1.5 - 2.0 * x) * (y) * (y - 0.5);
452 // midedge
453 derivs[4] = 8.0 * (1.0 - 2.0 * x) * (1.0 - y) * (0.5 - y);
454 derivs[5] =-8.0 * (0.5 - 2.0 * x) * (1.0 - y) * (y);
455 derivs[6] =-8.0 * (1.0 - 2.0 * x) * (y) * (0.5 - y);
456 derivs[7] = 8.0 * (2.0 * x - 1.5) * (1.0 - y) * (y);
457 // center
458 derivs[8] =16.0 * (1.0 - 2.0 * x) * (1.0 - y) * (y);
459
460 // Derivatives in the y-direction
461 // edge
462 derivs[9] = 4.0 * (1.0 - x) * (x - 0.5) * (1.5 - 2.0 * y);
463 derivs[10]=-4.0 * (x) * (x - 0.5) * (1.5 - 2.0 * y);
464 derivs[11]= 4.0 * (x) * (x - 0.5) * (2.0 * y - 0.5);
465 derivs[12]=-4.0 * (1.0 - x) * (x - 0.5) * (2.0 * y - 0.5);
466 // midedge
467 derivs[13]= 8.0 * (x) * (1.0 - x) * (2.0 * y - 1.5);
468 derivs[14]=-8.0 * (x) * (0.5 - x) * (1.0 - 2.0 * y);
469 derivs[15]=-8.0 * (x) * (1.0 - x) * (0.5 - 2.0 * y);
470 derivs[16]= 8.0 * (1.0 - x) * (0.5 - x) * (1.0 - 2.0 * y);
471 // center
472 derivs[17]=16.0 * (x) * (1.0 - x) * (1.0 - 2.0 * y);
473
474 }
475
476 //----------------------------------------------------------------------------
477 static double vtkQQuadCellPCoords[27] = { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0,
478 1.0, 1.0, 0.0, 0.0, 1.0, 0.0,
479 0.5, 0.0, 0.0, 1.0, 0.5, 0.0,
480 0.5, 1.0, 0.0, 0.0, 0.5, 0.0,
481 0.5, 0.5, 0.0 };
482
GetParametricCoords()483 double * vtkBiQuadraticQuad::GetParametricCoords ()
484 {
485 return vtkQQuadCellPCoords;
486 }
487
488 //----------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)489 void vtkBiQuadraticQuad::PrintSelf(ostream & os, vtkIndent indent)
490 {
491 this->Superclass::PrintSelf(os, indent);
492
493 os << indent << "Edge:\n";
494 this->Edge->PrintSelf (os, indent.GetNextIndent ());
495 os << indent << "Quad:\n";
496 this->Quad->PrintSelf (os, indent.GetNextIndent ());
497 os << indent << "Scalars:\n";
498 this->Scalars->PrintSelf (os, indent.GetNextIndent ());
499 }
500