1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkQuadraticLinearQuad.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 "vtkQuadraticLinearQuad.h"
20
21 #include "vtkObjectFactory.h"
22 #include "vtkDoubleArray.h"
23 #include "vtkMath.h"
24 #include "vtkQuad.h"
25 #include "vtkLine.h"
26 #include "vtkQuadraticEdge.h"
27 #include "vtkPoints.h"
28
29 vtkStandardNewMacro (vtkQuadraticLinearQuad);
30
31 //----------------------------------------------------------------------------
32 // Construct the quadratic linear quad with six points.
vtkQuadraticLinearQuad()33 vtkQuadraticLinearQuad::vtkQuadraticLinearQuad ()
34 {
35 this->Edge = vtkQuadraticEdge::New();
36 this->LinEdge = vtkLine::New ();
37 this->Quad = vtkQuad::New ();
38 this->Scalars = vtkDoubleArray::New ();
39 this->Scalars->SetNumberOfTuples (4); // vertices of a linear quad
40 this->Points->SetNumberOfPoints (6);
41 this->PointIds->SetNumberOfIds (6);
42 for (int i = 0; i < 6; i++)
43 {
44 this->Points->SetPoint(i, 0.0, 0.0, 0.0);
45 this->PointIds->SetId(i,0);
46 }
47 }
48
49 //----------------------------------------------------------------------------
~vtkQuadraticLinearQuad()50 vtkQuadraticLinearQuad::~vtkQuadraticLinearQuad ()
51 {
52 this->Edge->Delete();
53 this->LinEdge->Delete();
54 this->Quad->Delete();
55 this->Scalars->Delete();
56 }
57
58 //----------------------------------------------------------------------------
59 static int LinearQuads[2][4] = { {0, 4, 5, 3}, {4, 1, 2, 5} };
60
61 static int LinearQuadEdges[4][3] = { {0, 1, 4}, {1, 2,-1},
62 {2, 3, 5}, {3, 0,-1}};
63
64 //----------------------------------------------------------------------------
GetEdgeArray(int edgeId)65 int *vtkQuadraticLinearQuad::GetEdgeArray(int edgeId)
66 {
67 return LinearQuadEdges[edgeId];
68 }
69
70 //----------------------------------------------------------------------------
GetEdge(int edgeId)71 vtkCell *vtkQuadraticLinearQuad::GetEdge(int edgeId)
72 {
73 edgeId = (edgeId < 0 ? 0 : (edgeId > 3 ? 3 : edgeId));
74
75 // We have 2 linear edges
76 if (edgeId == 1 || edgeId == 3)
77 {
78 this->LinEdge->PointIds->SetId(0,this->PointIds->GetId(LinearQuadEdges[edgeId][0]));
79 this->LinEdge->PointIds->SetId(1,this->PointIds->GetId(LinearQuadEdges[edgeId][1]));
80
81 this->LinEdge->Points->SetPoint(0,this->Points->GetPoint(LinearQuadEdges[edgeId][0]));
82 this->LinEdge->Points->SetPoint(1,this->Points->GetPoint(LinearQuadEdges[edgeId][1]));
83
84 return this->LinEdge;
85 }
86 // and two quadratic edges
87 else // (edgeId == 0 || edgeId == 2)
88 {
89 this->Edge->PointIds->SetId(0,this->PointIds->GetId(LinearQuadEdges[edgeId][0]));
90 this->Edge->PointIds->SetId(1,this->PointIds->GetId(LinearQuadEdges[edgeId][1]));
91 this->Edge->PointIds->SetId(2,this->PointIds->GetId(LinearQuadEdges[edgeId][2]));
92
93 this->Edge->Points->SetPoint(0,this->Points->GetPoint(LinearQuadEdges[edgeId][0]));
94 this->Edge->Points->SetPoint(1,this->Points->GetPoint(LinearQuadEdges[edgeId][1]));
95 this->Edge->Points->SetPoint(2,this->Points->GetPoint(LinearQuadEdges[edgeId][2]));
96
97 return this->Edge;
98 }
99
100
101 }
102
103 //----------------------------------------------------------------------------
EvaluatePosition(double * x,double * closestPoint,int & subId,double pcoords[3],double & minDist2,double * weights)104 int vtkQuadraticLinearQuad::EvaluatePosition (double *x,
105 double *closestPoint,
106 int &subId, double pcoords[3],
107 double &minDist2, double *weights)
108 {
109 double pc[3], dist2;
110 int ignoreId, i, returnStatus = 0, status;
111 double tempWeights[4];
112 double closest[3];
113
114 // two linear quads are used
115 for (minDist2 = VTK_DOUBLE_MAX, i = 0; i < 2; i++)
116 {
117 this->Quad->Points->SetPoint(
118 0,this->Points->GetPoint(LinearQuads[i][0]));
119 this->Quad->Points->SetPoint(
120 1,this->Points->GetPoint(LinearQuads[i][1]));
121 this->Quad->Points->SetPoint(
122 2,this->Points->GetPoint(LinearQuads[i][2]));
123 this->Quad->Points->SetPoint(
124 3,this->Points->GetPoint(LinearQuads[i][3]));
125
126 status = this->Quad->EvaluatePosition(x,closest,ignoreId,pc,dist2,
127 tempWeights);
128 if ( status != -1 && dist2 < minDist2 )
129 {
130 returnStatus = status;
131 minDist2 = dist2;
132 subId = i;
133 pcoords[0] = pc[0];
134 pcoords[1] = pc[1];
135 }
136 }
137
138 // adjust parametric coordinates
139 if ( returnStatus != -1 )
140 {
141 if ( subId == 0 )
142 {
143 pcoords[0] /= 2.0;
144 }
145 else if ( subId == 1 )
146 {
147 pcoords[0] = 0.5 + (pcoords[0]/2.0);
148 }
149 pcoords[2] = 0.0;
150 if(closestPoint!=0)
151 {
152 // Compute both closestPoint and weights
153 this->EvaluateLocation(subId,pcoords,closestPoint,weights);
154 }
155 else
156 {
157 // Compute weigths only
158 this->InterpolationFunctions(pcoords,weights);
159 }
160 }
161
162 return returnStatus;
163 }
164
165 //----------------------------------------------------------------------------
EvaluateLocation(int & vtkNotUsed (subId),double pcoords[3],double x[3],double * weights)166 void vtkQuadraticLinearQuad::EvaluateLocation(int& vtkNotUsed(subId),
167 double pcoords[3],
168 double x[3], double *weights)
169 {
170 int i, j;
171 double pt[3];
172
173 this->InterpolationFunctions(pcoords,weights);
174
175 x[0] = x[1] = x[2] = 0.0;
176 for (i=0; i<6; i++)
177 {
178 this->Points->GetPoint(i, pt);
179 for (j=0; j<3; j++)
180 {
181 x[j] += pt[j] * weights[i];
182 }
183 }
184 }
185
186 //----------------------------------------------------------------------------
CellBoundary(int subId,double pcoords[3],vtkIdList * pts)187 int vtkQuadraticLinearQuad::CellBoundary(int subId, double pcoords[3], vtkIdList *pts)
188 {
189 return this->Quad->CellBoundary(subId, pcoords, pts);
190 }
191
192 //----------------------------------------------------------------------------
Contour(double value,vtkDataArray * cellScalars,vtkIncrementalPointLocator * locator,vtkCellArray * verts,vtkCellArray * lines,vtkCellArray * polys,vtkPointData * inPd,vtkPointData * outPd,vtkCellData * inCd,vtkIdType cellId,vtkCellData * outCd)193 void vtkQuadraticLinearQuad::Contour (double value,
194 vtkDataArray * cellScalars,
195 vtkIncrementalPointLocator * locator,
196 vtkCellArray * verts,
197 vtkCellArray * lines,
198 vtkCellArray * polys,
199 vtkPointData * inPd,
200 vtkPointData * outPd,
201 vtkCellData * inCd,
202 vtkIdType cellId,
203 vtkCellData * outCd)
204 {
205 for (int i = 0; i < 2; i++)
206 {
207 for (int j = 0; j < 4; j++)
208 {
209 this->Quad->Points->SetPoint(j,this->Points->GetPoint(LinearQuads[i][j]));
210 this->Quad->PointIds->SetId(j,this->PointIds->GetId(LinearQuads[i][j]));
211 this->Scalars->SetValue(j,cellScalars->GetTuple1(LinearQuads[i][j]));
212 }
213 this->Quad->Contour (value, this->Scalars, locator, verts, lines, polys,
214 inPd, outPd, inCd, cellId, outCd);
215 }
216 }
217
218 //----------------------------------------------------------------------------
219 // Clip this quadratic quad using scalar value provided. Like contouring,
220 // except that it cuts the quad to produce other quads and triangles.
Clip(double value,vtkDataArray * cellScalars,vtkIncrementalPointLocator * locator,vtkCellArray * polys,vtkPointData * inPd,vtkPointData * outPd,vtkCellData * inCd,vtkIdType cellId,vtkCellData * outCd,int insideOut)221 void vtkQuadraticLinearQuad::Clip (double value, vtkDataArray * cellScalars,
222 vtkIncrementalPointLocator * locator, vtkCellArray * polys,
223 vtkPointData * inPd, vtkPointData * outPd,
224 vtkCellData * inCd, vtkIdType cellId,
225 vtkCellData * outCd, int insideOut)
226 {
227 for (int i = 0; i < 2; i++)
228 {
229 for (int j = 0; j < 4; j++) //for each of the four vertices of the linear quad
230 {
231 this->Quad->Points->SetPoint(j,this->Points->GetPoint(LinearQuads[i][j]));
232 this->Quad->PointIds->SetId(j,this->PointIds->GetId(LinearQuads[i][j]));
233 this->Scalars->SetTuple(j,cellScalars->GetTuple(LinearQuads[i][j]));
234 }
235 this->Quad->Clip (value, this->Scalars, locator, polys, inPd, outPd, inCd,
236 cellId, outCd, insideOut);
237 }
238 }
239
240 //----------------------------------------------------------------------------
241 // Line-line intersection. Intersection has to occur within [0,1] parametric
242 // coordinates and with specified tolerance.
IntersectWithLine(double * p1,double * p2,double tol,double & t,double * x,double * pcoords,int & subId)243 int vtkQuadraticLinearQuad::IntersectWithLine (double *p1,
244 double *p2,
245 double tol,
246 double &t,
247 double *x,
248 double *pcoords,
249 int &subId)
250 {
251 int subTest, i;
252 subId = 0;
253
254 //intersect the two linear quads
255 for (i = 0; i < 2; i++)
256 {
257 this->Quad->Points->SetPoint(0,this->Points->GetPoint(LinearQuads[i][0]));
258 this->Quad->Points->SetPoint(1,this->Points->GetPoint(LinearQuads[i][1]));
259 this->Quad->Points->SetPoint(2,this->Points->GetPoint(LinearQuads[i][2]));
260 this->Quad->Points->SetPoint(3,this->Points->GetPoint(LinearQuads[i][3]));
261
262 if (this->Quad->IntersectWithLine (p1, p2, tol, t, x, pcoords, subTest))
263 {
264 return 1;
265 }
266 }
267
268 return 0;
269 }
270
271 //----------------------------------------------------------------------------
Triangulate(int vtkNotUsed (index),vtkIdList * ptIds,vtkPoints * pts)272 int vtkQuadraticLinearQuad::Triangulate (int vtkNotUsed (index),
273 vtkIdList * ptIds, vtkPoints * pts)
274 {
275 pts->Reset();
276 ptIds->Reset();
277
278 // Create four linear triangles:
279 // Choose the triangulation that minimizes the edge length
280 // across the cell.
281 double x0[3], x1[3], x2[3], x3[3], x4[3], x5[3];
282 this->Points->GetPoint (0, x0);
283 this->Points->GetPoint (1, x1);
284 this->Points->GetPoint (2, x2);
285 this->Points->GetPoint (3, x3);
286 this->Points->GetPoint (4, x4);
287 this->Points->GetPoint (5, x5);
288
289 if (vtkMath::Distance2BetweenPoints (x0, x5) <=
290 vtkMath::Distance2BetweenPoints (x3, x4))
291 {
292 ptIds->InsertId (0, this->PointIds->GetId (0));
293 ptIds->InsertId (1, this->PointIds->GetId (4));
294 ptIds->InsertId (2, this->PointIds->GetId (5));
295 pts->InsertPoint (0, this->Points->GetPoint (0));
296 pts->InsertPoint (1, this->Points->GetPoint (4));
297 pts->InsertPoint (2, this->Points->GetPoint (5));
298
299 ptIds->InsertId (3, this->PointIds->GetId (0));
300 ptIds->InsertId (4, this->PointIds->GetId (5));
301 ptIds->InsertId (5, this->PointIds->GetId (3));
302 pts->InsertPoint (3, this->Points->GetPoint (0));
303 pts->InsertPoint (4, this->Points->GetPoint (5));
304 pts->InsertPoint (5, this->Points->GetPoint (3));
305 }
306 else
307 {
308 ptIds->InsertId (0, this->PointIds->GetId (0));
309 ptIds->InsertId (1, this->PointIds->GetId (4));
310 ptIds->InsertId (2, this->PointIds->GetId (3));
311 pts->InsertPoint (0, this->Points->GetPoint (0));
312 pts->InsertPoint (1, this->Points->GetPoint (4));
313 pts->InsertPoint (2, this->Points->GetPoint (3));
314
315 ptIds->InsertId (3, this->PointIds->GetId (4));
316 ptIds->InsertId (4, this->PointIds->GetId (5));
317 ptIds->InsertId (5, this->PointIds->GetId (3));
318 pts->InsertPoint (3, this->Points->GetPoint (4));
319 pts->InsertPoint (4, this->Points->GetPoint (5));
320 pts->InsertPoint (5, this->Points->GetPoint (3));
321 }
322
323 if (vtkMath::Distance2BetweenPoints (x4, x2) <=
324 vtkMath::Distance2BetweenPoints (x5, x1))
325 {
326 ptIds->InsertId (6, this->PointIds->GetId (4));
327 ptIds->InsertId (7, this->PointIds->GetId (1));
328 ptIds->InsertId (8, this->PointIds->GetId (2));
329 pts->InsertPoint (6, this->Points->GetPoint (4));
330 pts->InsertPoint (7, this->Points->GetPoint (1));
331 pts->InsertPoint (8, this->Points->GetPoint (2));
332
333 ptIds->InsertId (9, this->PointIds->GetId (4));
334 ptIds->InsertId (10, this->PointIds->GetId (2));
335 ptIds->InsertId (11, this->PointIds->GetId (5));
336 pts->InsertPoint (9, this->Points->GetPoint (4));
337 pts->InsertPoint (10, this->Points->GetPoint (2));
338 pts->InsertPoint (11, this->Points->GetPoint (5));
339 }
340 else
341 {
342 ptIds->InsertId (6, this->PointIds->GetId (4));
343 ptIds->InsertId (7, this->PointIds->GetId (1));
344 ptIds->InsertId (8, this->PointIds->GetId (5));
345 pts->InsertPoint (6, this->Points->GetPoint (4));
346 pts->InsertPoint (7, this->Points->GetPoint (1));
347 pts->InsertPoint (8, this->Points->GetPoint (5));
348
349 ptIds->InsertId (9, this->PointIds->GetId (1));
350 ptIds->InsertId (10, this->PointIds->GetId (2));
351 ptIds->InsertId (11, this->PointIds->GetId (5));
352 pts->InsertPoint (9, this->Points->GetPoint (1));
353 pts->InsertPoint (10, this->Points->GetPoint (2));
354 pts->InsertPoint (11, this->Points->GetPoint (5));
355 }
356
357 return 1;
358 }
359
360 //----------------------------------------------------------------------------
Derivatives(int vtkNotUsed (subId),double pcoords[3],double * values,int dim,double * derivs)361 void vtkQuadraticLinearQuad::Derivatives (int vtkNotUsed (subId),
362 double pcoords[3], double *values, int dim, double *derivs)
363 {
364 double x0[3], x1[3], x2[3], deltaX[3], weights[6];
365 int i, j;
366 double functionDerivs[12];
367
368 this->Points->GetPoint (0, x0);
369 this->Points->GetPoint (1, x1);
370 this->Points->GetPoint (2, x2);
371
372 this->InterpolationFunctions (pcoords, weights);
373 this->InterpolationDerivs (pcoords, functionDerivs);
374
375 for (i = 0; i < 3; i++)
376 {
377 deltaX[i] = x1[i] - x0[i] - x2[i];
378 }
379 for (i = 0; i < dim; i++)
380 {
381 for (j = 0; j < 3; j++)
382 {
383 if (deltaX[j] != 0)
384 {
385 derivs[3 * i + j] = (values[2 * i + 1] - values[2 * i]) / deltaX[j];
386 }
387 else
388 {
389 derivs[3 * i + j] = 0;
390 }
391 }
392 }
393 }
394
395
396
397 //----------------------------------------------------------------------------
398 // Compute interpolation functions. The first four nodes are the corner
399 // vertices; the others are mid-edge nodes.
InterpolationFunctions(double pcoords[3],double weights[6])400 void vtkQuadraticLinearQuad::InterpolationFunctions(double pcoords[3],
401 double weights[6])
402 {
403 double x = pcoords[0];
404 double y = pcoords[1];
405
406 //corners
407 weights[0] = -1.0 *(2.0 *x - 1.0) * (x - 1.0) * (y - 1.0);
408 weights[1] = -1.0 *(2.0 *x - 1.0) * (x) * (y - 1.0);
409 weights[2] = (2.0 *x - 1.0) * (x) * (y);
410 weights[3] = (2.0 *x - 1.0) * (x - 1.0 ) * (y);
411
412 //Edge middle nodes
413 weights[4] = 4.0 * (x) * (1.0 - x) * (1.0 - y);
414 weights[5] = 4.0 * (x) * (1.0 - x) * (y);
415
416
417 }
418
419 //----------------------------------------------------------------------------
420 // Derivatives in parametric space.
InterpolationDerivs(double pcoords[3],double derivs[12])421 void vtkQuadraticLinearQuad::InterpolationDerivs(double pcoords[3], double derivs[12])
422 {
423 double x = pcoords[0];
424 double y = pcoords[1];
425
426 // Derivatives in the x-direction
427 //corners
428 derivs[0] = -1.0 * (4.0 * x - 3.0) * (y - 1.0);
429 derivs[1] = -1.0 * (4.0 * x - 1.0) * (y - 1.0);
430 derivs[2] = (4.0 * x - 1.0) * (y);
431 derivs[3] = (4.0 * x - 3.0) * (y);
432 //Edge middle nodes
433 derivs[4] = 4.0 * (1.0 - 2.0 * x) * (1.0 - y);
434 derivs[5] = 4.0 * (1.0 - 2.0 * x) * (y);
435
436 // Derivatives in the y-direction
437 //corners
438 derivs[6] = -1.0 * (2.0 * x - 1.0) * (x - 1.0);
439 derivs[7] = -1.0 * (2.0 * x - 1.0) * (x) ;
440 derivs[8] = (2.0 * x - 1.0) * (x) ;
441 derivs[9] = (2.0 * x - 1.0) * (x - 1.0);
442 //Edge middle nodes
443 derivs[10]= -4.0 * (x) * (1.0 - x);
444 derivs[11]= 4.0 * (x) * (1.0 - x);
445
446 }
447
448 //----------------------------------------------------------------------------
449 static double vtkQLinQuadCellPCoords[18] = { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0,
450 1.0, 1.0, 0.0, 0.0, 1.0, 0.0,
451 0.5, 0.0, 0.0, 0.5, 1.0, 0.0 };
452
GetParametricCoords()453 double *vtkQuadraticLinearQuad::GetParametricCoords ()
454 {
455 return vtkQLinQuadCellPCoords;
456 }
457
458 //----------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)459 void vtkQuadraticLinearQuad::PrintSelf (ostream & os, vtkIndent indent)
460 {
461 this->Superclass::PrintSelf(os,indent);
462
463 os << indent << "Edge:\n";
464 this->Edge->PrintSelf(os,indent.GetNextIndent());
465 os << indent << "Quad:\n";
466 this->Quad->PrintSelf(os,indent.GetNextIndent());
467 os << indent << "Scalars:\n";
468 this->Scalars->PrintSelf(os,indent.GetNextIndent());
469 }
470
471