1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkQuad.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 #include "vtkQuad.h"
16
17 #include "vtkObjectFactory.h"
18 #include "vtkCellArray.h"
19 #include "vtkCellData.h"
20 #include "vtkLine.h"
21 #include "vtkTriangle.h"
22 #include "vtkMath.h"
23 #include "vtkPlane.h"
24 #include "vtkPointData.h"
25 #include "vtkIncrementalPointLocator.h"
26 #include "vtkPoints.h"
27
28 vtkStandardNewMacro(vtkQuad);
29
30 static const double VTK_DIVERGED = 1.e6;
31
32 //----------------------------------------------------------------------------
33 // Construct the quad with four points.
vtkQuad()34 vtkQuad::vtkQuad()
35 {
36 this->Points->SetNumberOfPoints(4);
37 this->PointIds->SetNumberOfIds(4);
38 for (int i = 0; i < 4; i++)
39 {
40 this->Points->SetPoint(i, 0.0, 0.0, 0.0);
41 this->PointIds->SetId(i,0);
42 }
43 this->Line = vtkLine::New();
44 this->Triangle = vtkTriangle::New();
45 }
46
47 //----------------------------------------------------------------------------
~vtkQuad()48 vtkQuad::~vtkQuad()
49 {
50 this->Line->Delete();
51 this->Triangle->Delete();
52 }
53
54 //----------------------------------------------------------------------------
55 static const int VTK_QUAD_MAX_ITERATION=20;
56 static const double VTK_QUAD_CONVERGED=1.e-04;
57
ComputeNormal(vtkQuad * self,double pt1[3],double pt2[3],double pt3[3],double n[3])58 inline static void ComputeNormal(vtkQuad *self, double pt1[3], double pt2[3],
59 double pt3[3], double n[3])
60 {
61 vtkTriangle::ComputeNormal (pt1, pt2, pt3, n);
62
63 // If first three points are co-linear, then use fourth point
64 //
65 double pt4[3];
66 if ( n[0] == 0.0 && n[1] == 0.0 && n[2] == 0.0 )
67 {
68 self->Points->GetPoint(3,pt4);
69 vtkTriangle::ComputeNormal (pt2, pt3, pt4, n);
70 }
71 }
72
73 //----------------------------------------------------------------------------
EvaluatePosition(const double x[3],double closestPoint[3],int & subId,double pcoords[3],double & dist2,double weights[])74 int vtkQuad::EvaluatePosition(const double x[3], double closestPoint[3],
75 int& subId, double pcoords[3],
76 double& dist2, double weights[])
77 {
78 int i, j;
79 double pt1[3], pt2[3], pt3[3], pt[3], n[3];
80 double det;
81 double maxComponent;
82 int idx=0, indices[2];
83 int iteration, converged;
84 double params[2];
85 double fcol[2], rcol[2], scol[2], cp[3];
86 double derivs[8];
87
88 subId = 0;
89 pcoords[0] = pcoords[1] = params[0] = params[1] = 0.5;
90 pcoords[2] = 0.0;
91
92 // Get normal for quadrilateral
93 //
94 this->Points->GetPoint(0, pt1);
95 this->Points->GetPoint(1, pt2);
96 this->Points->GetPoint(2, pt3);
97 ComputeNormal (this, pt1, pt2, pt3, n);
98
99 // Project point to plane
100 //
101 vtkPlane::ProjectPoint(x,pt1,n,cp);
102
103 // Construct matrices. Since we have over determined system, need to find
104 // which 2 out of 3 equations to use to develop equations. (Any 2 should
105 // work since we've projected point to plane.)
106 //
107 for (maxComponent=0.0, i=0; i<3; i++)
108 {
109 if (fabs(n[i]) > maxComponent)
110 {
111 maxComponent = fabs(n[i]);
112 idx = i;
113 }
114 }
115 for (j=0, i=0; i<3; i++)
116 {
117 if ( i != idx )
118 {
119 indices[j++] = i;
120 }
121 }
122
123 // Use Newton's method to solve for parametric coordinates
124 //
125 for (iteration=converged=0; !converged
126 && (iteration < VTK_QUAD_MAX_ITERATION);
127 iteration++)
128 {
129 // calculate element interpolation functions and derivatives
130 //
131 this->InterpolationFunctions(pcoords, weights);
132 this->InterpolationDerivs(pcoords, derivs);
133
134 // calculate newton functions
135 //
136 for (i=0; i<2; i++)
137 {
138 fcol[i] = rcol[i] = scol[i] = 0.0;
139 }
140 for (i=0; i<4; i++)
141 {
142 this->Points->GetPoint(i, pt);
143 for (j=0; j<2; j++)
144 {
145 fcol[j] += pt[indices[j]] * weights[i];
146 rcol[j] += pt[indices[j]] * derivs[i];
147 scol[j] += pt[indices[j]] * derivs[i+4];
148 }
149 }
150
151 for (j=0; j<2; j++)
152 {
153 fcol[j] -= cp[indices[j]];
154 }
155
156 // compute determinants and generate improvements
157 //
158 if ( (det=vtkMath::Determinant2x2(rcol,scol)) == 0.0 )
159 {
160 return -1;
161 }
162
163 pcoords[0] = params[0] - vtkMath::Determinant2x2 (fcol,scol) / det;
164 pcoords[1] = params[1] - vtkMath::Determinant2x2 (rcol,fcol) / det;
165
166 // check for convergence
167 //
168 if ( ((fabs(pcoords[0]-params[0])) < VTK_QUAD_CONVERGED) &&
169 ((fabs(pcoords[1]-params[1])) < VTK_QUAD_CONVERGED) )
170 {
171 converged = 1;
172 }
173 // Test for bad divergence (S.Hirschberg 11.12.2001)
174 else if ((fabs(pcoords[0]) > VTK_DIVERGED) ||
175 (fabs(pcoords[1]) > VTK_DIVERGED))
176 {
177 return -1;
178 }
179
180 // if not converged, repeat
181 //
182 else
183 {
184 params[0] = pcoords[0];
185 params[1] = pcoords[1];
186 }
187 }
188
189 // if not converged, set the parametric coordinates to arbitrary values
190 // outside of element
191 //
192 if ( !converged )
193 {
194 return -1;
195 }
196
197 this->InterpolationFunctions(pcoords, weights);
198
199 if ( pcoords[0] >= -0.001 && pcoords[0] <= 1.001 &&
200 pcoords[1] >= -0.001 && pcoords[1] <= 1.001 )
201 {
202 if (closestPoint)
203 {
204 dist2 =
205 vtkMath::Distance2BetweenPoints(cp,x); //projection distance
206 closestPoint[0] = cp[0];
207 closestPoint[1] = cp[1];
208 closestPoint[2] = cp[2];
209 }
210 return 1;
211 }
212 else
213 {
214 double t;
215 double pt4[3];
216
217 if (closestPoint)
218 {
219 this->Points->GetPoint(3, pt4);
220
221 if ( pcoords[0] < 0.0 && pcoords[1] < 0.0 )
222 {
223 dist2 = vtkMath::Distance2BetweenPoints(x,pt1);
224 for (i=0; i<3; i++)
225 {
226 closestPoint[i] = pt1[i];
227 }
228 }
229 else if ( pcoords[0] > 1.0 && pcoords[1] < 0.0 )
230 {
231 dist2 = vtkMath::Distance2BetweenPoints(x,pt2);
232 for (i=0; i<3; i++)
233 {
234 closestPoint[i] = pt2[i];
235 }
236 }
237 else if ( pcoords[0] > 1.0 && pcoords[1] > 1.0 )
238 {
239 dist2 = vtkMath::Distance2BetweenPoints(x,pt3);
240 for (i=0; i<3; i++)
241 {
242 closestPoint[i] = pt3[i];
243 }
244 }
245 else if ( pcoords[0] < 0.0 && pcoords[1] > 1.0 )
246 {
247 dist2 = vtkMath::Distance2BetweenPoints(x,pt4);
248 for (i=0; i<3; i++)
249 {
250 closestPoint[i] = pt4[i];
251 }
252 }
253 else if ( pcoords[0] < 0.0 )
254 {
255 dist2 = vtkLine::DistanceToLine(x,pt1,pt4,t,closestPoint);
256 }
257 else if ( pcoords[0] > 1.0 )
258 {
259 dist2 = vtkLine::DistanceToLine(x,pt2,pt3,t,closestPoint);
260 }
261 else if ( pcoords[1] < 0.0 )
262 {
263 dist2 = vtkLine::DistanceToLine(x,pt1,pt2,t,closestPoint);
264 }
265 else if ( pcoords[1] > 1.0 )
266 {
267 dist2 = vtkLine::DistanceToLine(x,pt3,pt4,t,closestPoint);
268 }
269 }
270 return 0;
271 }
272 }
273
274 //----------------------------------------------------------------------------
EvaluateLocation(int & vtkNotUsed (subId),const double pcoords[3],double x[3],double * weights)275 void vtkQuad::EvaluateLocation(int& vtkNotUsed(subId), const double pcoords[3],
276 double x[3], double *weights)
277 {
278 int i, j;
279 double pt[3];
280
281 this->InterpolationFunctions(pcoords, weights);
282
283 x[0] = x[1] = x[2] = 0.0;
284 for (i=0; i<4; i++)
285 {
286 this->Points->GetPoint(i, pt);
287 for (j=0; j<3; j++)
288 {
289 x[j] += pt[j] * weights[i];
290 }
291 }
292 }
293
294 //----------------------------------------------------------------------------
295 // Compute iso-parametric interpolation functions
296 //
InterpolationFunctions(const double pcoords[3],double sf[4])297 void vtkQuad::InterpolationFunctions(const double pcoords[3], double sf[4])
298 {
299 double rm, sm;
300
301 rm = 1. - pcoords[0];
302 sm = 1. - pcoords[1];
303
304 sf[0] = rm * sm;
305 sf[1] = pcoords[0] * sm;
306 sf[2] = pcoords[0] * pcoords[1];
307 sf[3] = rm * pcoords[1];
308 }
309
310 //----------------------------------------------------------------------------
InterpolationDerivs(const double pcoords[3],double derivs[8])311 void vtkQuad::InterpolationDerivs(const double pcoords[3], double derivs[8])
312 {
313 double rm, sm;
314
315 rm = 1. - pcoords[0];
316 sm = 1. - pcoords[1];
317
318 derivs[0] = -sm;
319 derivs[1] = sm;
320 derivs[2] = pcoords[1];
321 derivs[3] = -pcoords[1];
322 derivs[4] = -rm;
323 derivs[5] = -pcoords[0];
324 derivs[6] = pcoords[0];
325 derivs[7] = rm;
326 }
327
328 //----------------------------------------------------------------------------
CellBoundary(int vtkNotUsed (subId),const double pcoords[3],vtkIdList * pts)329 int vtkQuad::CellBoundary(int vtkNotUsed(subId), const double pcoords[3],
330 vtkIdList *pts)
331 {
332 double t1=pcoords[0]-pcoords[1];
333 double t2=1.0-pcoords[0]-pcoords[1];
334
335 pts->SetNumberOfIds(2);
336
337 // compare against two lines in parametric space that divide element
338 // into four pieces.
339 if ( t1 >= 0.0 && t2 >= 0.0 )
340 {
341 pts->SetId(0,this->PointIds->GetId(0));
342 pts->SetId(1,this->PointIds->GetId(1));
343 }
344
345 else if ( t1 >= 0.0 && t2 < 0.0 )
346 {
347 pts->SetId(0,this->PointIds->GetId(1));
348 pts->SetId(1,this->PointIds->GetId(2));
349 }
350
351 else if ( t1 < 0.0 && t2 < 0.0 )
352 {
353 pts->SetId(0,this->PointIds->GetId(2));
354 pts->SetId(1,this->PointIds->GetId(3));
355 }
356
357 else //( t1 < 0.0 && t2 >= 0.0 )
358 {
359 pts->SetId(0,this->PointIds->GetId(3));
360 pts->SetId(1,this->PointIds->GetId(0));
361 }
362
363 if ( pcoords[0] < 0.0 || pcoords[0] > 1.0 ||
364 pcoords[1] < 0.0 || pcoords[1] > 1.0 )
365 {
366 return 0;
367 }
368 else
369 {
370 return 1;
371 }
372 }
373
374 //----------------------------------------------------------------------------
375 // Marching (convex) quadrilaterals
376 //
377 namespace { //required so we don't violate ODR
378 static int edges[4][2] = { {0,1}, {1,2}, {3,2}, {0,3} };
379
380 typedef int EDGE_LIST;
381 typedef struct {
382 EDGE_LIST edges[5];
383 } LINE_CASES;
384
385 static LINE_CASES lineCases[] = {
386 {{-1, -1, -1, -1, -1}},
387 {{0, 3, -1, -1, -1}},
388 {{1, 0, -1, -1, -1}},
389 {{1, 3, -1, -1, -1}},
390 {{2, 1, -1, -1, -1}},
391 {{0, 3, 2, 1, -1}},
392 {{2, 0, -1, -1, -1}},
393 {{2, 3, -1, -1, -1}},
394 {{3, 2, -1, -1, -1}},
395 {{0, 2, -1, -1, -1}},
396 {{1, 0, 3, 2, -1}},
397 {{1, 2, -1, -1, -1}},
398 {{3, 1, -1, -1, -1}},
399 {{0, 1, -1, -1, -1}},
400 {{3, 0, -1, -1, -1}},
401 {{-1, -1, -1, -1, -1}}
402 };
403 }
404
405 //----------------------------------------------------------------------------
GetEdgeArray(int edgeId)406 int *vtkQuad::GetEdgeArray(int edgeId)
407 {
408 return edges[edgeId];
409 }
410
411 //----------------------------------------------------------------------------
Contour(double value,vtkDataArray * cellScalars,vtkIncrementalPointLocator * locator,vtkCellArray * verts,vtkCellArray * lines,vtkCellArray * vtkNotUsed (polys),vtkPointData * inPd,vtkPointData * outPd,vtkCellData * inCd,vtkIdType cellId,vtkCellData * outCd)412 void vtkQuad::Contour(double value, vtkDataArray *cellScalars,
413 vtkIncrementalPointLocator *locator,
414 vtkCellArray *verts,
415 vtkCellArray *lines,
416 vtkCellArray *vtkNotUsed(polys),
417 vtkPointData *inPd, vtkPointData *outPd,
418 vtkCellData *inCd, vtkIdType cellId, vtkCellData *outCd)
419 {
420 static const int CASE_MASK[4] = {1,2,4,8};
421 LINE_CASES *lineCase;
422 EDGE_LIST *edge;
423 int i, j, index, *vert;
424 int newCellId;
425 vtkIdType pts[2];
426 int e1, e2;
427 double t, x1[3], x2[3], x[3], deltaScalar;
428 vtkIdType offset = verts->GetNumberOfCells();
429
430 // Build the case table
431 for ( i=0, index = 0; i < 4; i++)
432 {
433 if (cellScalars->GetComponent(i,0) >= value)
434 {
435 index |= CASE_MASK[i];
436 }
437 }
438
439 lineCase = lineCases + index;
440 edge = lineCase->edges;
441
442 for ( ; edge[0] > -1; edge += 2 )
443 {
444 for (i=0; i<2; i++) // insert line
445 {
446 vert = edges[edge[i]];
447 // calculate a preferred interpolation direction
448 deltaScalar = (cellScalars->GetComponent(vert[1],0)
449 - cellScalars->GetComponent(vert[0],0));
450 if (deltaScalar > 0)
451 {
452 e1 = vert[0]; e2 = vert[1];
453 }
454 else
455 {
456 e1 = vert[1]; e2 = vert[0];
457 deltaScalar = -deltaScalar;
458 }
459
460 // linear interpolation
461 if (deltaScalar == 0.0)
462 {
463 t = 0.0;
464 }
465 else
466 {
467 t = (value - cellScalars->GetComponent(e1,0)) / deltaScalar;
468 }
469
470 this->Points->GetPoint(e1, x1);
471 this->Points->GetPoint(e2, x2);
472
473 for (j=0; j<3; j++)
474 {
475 x[j] = x1[j] + t * (x2[j] - x1[j]);
476 }
477 if ( locator->InsertUniquePoint(x, pts[i]) )
478 {
479 if ( outPd )
480 {
481 vtkIdType p1 = this->PointIds->GetId(e1);
482 vtkIdType p2 = this->PointIds->GetId(e2);
483 outPd->InterpolateEdge(inPd,pts[i],p1,p2,t);
484 }
485 }
486 }
487 // check for degenerate line
488 if ( pts[0] != pts[1] )
489 {
490 newCellId = offset + lines->InsertNextCell(2,pts);
491 if (outCd)
492 {
493 outCd->CopyData(inCd, cellId, newCellId);
494 }
495 }
496 }
497 }
498
499 //----------------------------------------------------------------------------
GetEdge(int edgeId)500 vtkCell *vtkQuad::GetEdge(int edgeId)
501 {
502 int edgeIdPlus1 = edgeId + 1;
503
504 if (edgeIdPlus1 > 3)
505 {
506 edgeIdPlus1 = 0;
507 }
508
509 // load point id's
510 this->Line->PointIds->SetId(0,this->PointIds->GetId(edgeId));
511 this->Line->PointIds->SetId(1,this->PointIds->GetId(edgeIdPlus1));
512
513 // load coordinates
514 this->Line->Points->SetPoint(0,this->Points->GetPoint(edgeId));
515 this->Line->Points->SetPoint(1,this->Points->GetPoint(edgeIdPlus1));
516
517 return this->Line;
518 }
519
520
521 //----------------------------------------------------------------------------
522 // Intersect plane; see whether point is in quadrilateral. This code
523 // splits the quad into two triangles and intersects them (because the
524 // quad may be non-planar).
525 //
IntersectWithLine(const double p1[3],const double p2[3],double tol,double & t,double x[3],double pcoords[3],int & subId)526 int vtkQuad::IntersectWithLine(const double p1[3], const double p2[3], double tol, double& t,
527 double x[3], double pcoords[3], int& subId)
528 {
529 int diagonalCase;
530 double d1 = vtkMath::Distance2BetweenPoints(this->Points->GetPoint(0),
531 this->Points->GetPoint(2));
532 double d2 = vtkMath::Distance2BetweenPoints(this->Points->GetPoint(1),
533 this->Points->GetPoint(3));
534 subId = 0;
535
536 // Figure out how to uniquely tessellate the quad. Watch out for
537 // equivalent triangulations (i.e., the triangulation is equivalent
538 // no matter where the diagonal). In this case use the point ids as
539 // a tie breaker to insure unique triangulation across the quad.
540 //
541 if ( d1 == d2 ) //rare case; discriminate based on point id
542 {
543 int i, id, maxId=0, maxIdx=0;
544 for (i=0; i<4; i++) //find the maximum id
545 {
546 if ( (id=this->PointIds->GetId(i)) > maxId )
547 {
548 maxId = id;
549 maxIdx = i;
550 }
551 }
552 if ( maxIdx == 0 || maxIdx == 2) diagonalCase = 0;
553 else diagonalCase = 1;
554 }
555 else if ( d1 < d2 )
556 {
557 diagonalCase = 0;
558 }
559 else //d2 < d1
560 {
561 diagonalCase = 1;
562 }
563
564 // Note: in the following code the parametric coords must be adjusted to
565 // reflect the use of the triangle parametric coordinate system.
566 switch (diagonalCase)
567 {
568 case 0:
569 this->Triangle->Points->SetPoint(0,this->Points->GetPoint(0));
570 this->Triangle->Points->SetPoint(1,this->Points->GetPoint(1));
571 this->Triangle->Points->SetPoint(2,this->Points->GetPoint(2));
572 if (this->Triangle->IntersectWithLine(p1, p2, tol, t, x, pcoords, subId) )
573 {
574 pcoords[0] = pcoords[0] + pcoords[1];
575 return 1;
576 }
577 this->Triangle->Points->SetPoint(0,this->Points->GetPoint(2));
578 this->Triangle->Points->SetPoint(1,this->Points->GetPoint(3));
579 this->Triangle->Points->SetPoint(2,this->Points->GetPoint(0));
580 if (this->Triangle->IntersectWithLine(p1, p2, tol, t, x, pcoords, subId) )
581 {
582 pcoords[0] = 1.0 - (pcoords[0]+pcoords[1]);
583 pcoords[1] = 1.0 - pcoords[1];
584 return 1;
585 }
586 return 0;
587
588 case 1:
589 this->Triangle->Points->SetPoint(0,this->Points->GetPoint(0));
590 this->Triangle->Points->SetPoint(1,this->Points->GetPoint(1));
591 this->Triangle->Points->SetPoint(2,this->Points->GetPoint(3));
592 if (this->Triangle->IntersectWithLine(p1, p2, tol, t, x, pcoords, subId) )
593 {
594 return 1;
595 }
596 this->Triangle->Points->SetPoint(0,this->Points->GetPoint(2));
597 this->Triangle->Points->SetPoint(1,this->Points->GetPoint(3));
598 this->Triangle->Points->SetPoint(2,this->Points->GetPoint(1));
599 if (this->Triangle->IntersectWithLine(p1, p2, tol, t, x, pcoords, subId) )
600 {
601 pcoords[0] = 1.0 - pcoords[0];
602 pcoords[1] = 1.0 - pcoords[1];
603 return 1;
604 }
605
606 return 0;
607 }
608
609 return 0;
610 }
611
612 //----------------------------------------------------------------------------
Triangulate(int vtkNotUsed (index),vtkIdList * ptIds,vtkPoints * pts)613 int vtkQuad::Triangulate(int vtkNotUsed(index), vtkIdList *ptIds,
614 vtkPoints *pts)
615 {
616 double d1, d2;
617
618 pts->Reset();
619 ptIds->Reset();
620
621 // use minimum diagonal (Delaunay triangles) - assumed convex
622 d1 = vtkMath::Distance2BetweenPoints(this->Points->GetPoint(0),
623 this->Points->GetPoint(2));
624 d2 = vtkMath::Distance2BetweenPoints(this->Points->GetPoint(1),
625 this->Points->GetPoint(3));
626
627 if ( d1 <= d2 )
628 {
629 ptIds->InsertId(0,this->PointIds->GetId(0));
630 pts->InsertPoint(0,this->Points->GetPoint(0));
631 ptIds->InsertId(1,this->PointIds->GetId(1));
632 pts->InsertPoint(1,this->Points->GetPoint(1));
633 ptIds->InsertId(2,this->PointIds->GetId(2));
634 pts->InsertPoint(2,this->Points->GetPoint(2));
635
636 ptIds->InsertId(3,this->PointIds->GetId(0));
637 pts->InsertPoint(3,this->Points->GetPoint(0));
638 ptIds->InsertId(4,this->PointIds->GetId(2));
639 pts->InsertPoint(4,this->Points->GetPoint(2));
640 ptIds->InsertId(5,this->PointIds->GetId(3));
641 pts->InsertPoint(5,this->Points->GetPoint(3));
642 }
643 else
644 {
645 ptIds->InsertId(0,this->PointIds->GetId(0));
646 pts->InsertPoint(0,this->Points->GetPoint(0));
647 ptIds->InsertId(1,this->PointIds->GetId(1));
648 pts->InsertPoint(1,this->Points->GetPoint(1));
649 ptIds->InsertId(2,this->PointIds->GetId(3));
650 pts->InsertPoint(2,this->Points->GetPoint(3));
651
652 ptIds->InsertId(3,this->PointIds->GetId(1));
653 pts->InsertPoint(3,this->Points->GetPoint(1));
654 ptIds->InsertId(4,this->PointIds->GetId(2));
655 pts->InsertPoint(4,this->Points->GetPoint(2));
656 ptIds->InsertId(5,this->PointIds->GetId(3));
657 pts->InsertPoint(5,this->Points->GetPoint(3));
658 }
659
660 return 1;
661 }
662
663 //----------------------------------------------------------------------------
Derivatives(int vtkNotUsed (subId),const double pcoords[3],const double * values,int dim,double * derivs)664 void vtkQuad::Derivatives(int vtkNotUsed(subId), const double pcoords[3],
665 const double *values, int dim, double *derivs)
666 {
667 double v0[2], v1[2], v2[2], v3[2], v10[3], v20[3], lenX;
668 double x0[3], x1[3], x2[3], x3[3], n[3], vec20[3], vec30[3];
669 double *J[2], J0[2], J1[2];
670 double *JI[2], JI0[2], JI1[2];
671 double funcDerivs[8], sum[2], dBydx, dBydy;
672 int i, j;
673
674 // Project points of quad into 2D system
675 this->Points->GetPoint(0, x0);
676 this->Points->GetPoint(1, x1);
677 this->Points->GetPoint(2, x2);
678 ComputeNormal (this,x0, x1, x2, n);
679 this->Points->GetPoint(3, x3);
680
681 for (i=0; i < 3; i++)
682 {
683 v10[i] = x1[i] - x0[i];
684 vec20[i] = x2[i] - x0[i];
685 vec30[i] = x3[i] - x0[i];
686 }
687
688 vtkMath::Cross(n,v10,v20); //creates local y' axis
689
690 if ( (lenX=vtkMath::Normalize(v10)) <= 0.0
691 || vtkMath::Normalize(v20) <= 0.0 ) //degenerate
692 {
693 for ( j=0; j < dim; j++ )
694 {
695 for ( i=0; i < 3; i++ )
696 {
697 derivs[j*dim + i] = 0.0;
698 }
699 }
700 return;
701 }
702
703 v0[0] = v0[1] = 0.0; //convert points to 2D (i.e., local system)
704 v1[0] = lenX; v1[1] = 0.0;
705 v2[0] = vtkMath::Dot(vec20,v10);
706 v2[1] = vtkMath::Dot(vec20,v20);
707 v3[0] = vtkMath::Dot(vec30,v10);
708 v3[1] = vtkMath::Dot(vec30,v20);
709
710 this->InterpolationDerivs(pcoords, funcDerivs);
711
712 // Compute Jacobian and inverse Jacobian
713 J[0] = J0; J[1] = J1;
714 JI[0] = JI0; JI[1] = JI1;
715
716 J[0][0] = v0[0]*funcDerivs[0] + v1[0]*funcDerivs[1] +
717 v2[0]*funcDerivs[2] + v3[0]*funcDerivs[3];
718 J[0][1] = v0[1]*funcDerivs[0] + v1[1]*funcDerivs[1] +
719 v2[1]*funcDerivs[2] + v3[1]*funcDerivs[3];
720 J[1][0] = v0[0]*funcDerivs[4] + v1[0]*funcDerivs[5] +
721 v2[0]*funcDerivs[6] + v3[0]*funcDerivs[7];
722 J[1][1] = v0[1]*funcDerivs[4] + v1[1]*funcDerivs[5] +
723 v2[1]*funcDerivs[6] + v3[1]*funcDerivs[7];
724
725 // Compute inverse Jacobian, return if Jacobian is singular
726 if (!vtkMath::InvertMatrix(J,JI,2))
727 {
728 for ( j=0; j < dim; j++ )
729 {
730 for ( i=0; i < 3; i++ )
731 {
732 derivs[j*dim + i] = 0.0;
733 }
734 }
735 return;
736 }
737
738 // Loop over "dim" derivative values. For each set of values,
739 // compute derivatives
740 // in local system and then transform into modelling system.
741 // First compute derivatives in local x'-y' coordinate system
742 for ( j=0; j < dim; j++ )
743 {
744 sum[0] = sum[1] = 0.0;
745 for ( i=0; i < 4; i++) //loop over interp. function derivatives
746 {
747 sum[0] += funcDerivs[i] * values[dim*i + j];
748 sum[1] += funcDerivs[4 + i] * values[dim*i + j];
749 }
750 dBydx = sum[0]*JI[0][0] + sum[1]*JI[0][1];
751 dBydy = sum[0]*JI[1][0] + sum[1]*JI[1][1];
752
753 // Transform into global system (dot product with global axes)
754 derivs[3*j] = dBydx * v10[0] + dBydy * v20[0];
755 derivs[3*j + 1] = dBydx * v10[1] + dBydy * v20[1];
756 derivs[3*j + 2] = dBydx * v10[2] + dBydy * v20[2];
757 }
758 }
759
760 //----------------------------------------------------------------------------
761 // support quad clipping
762 typedef int QUAD_EDGE_LIST;
763 typedef struct {
764 QUAD_EDGE_LIST edges[14];
765 } QUAD_CASES;
766
767 static QUAD_CASES quadCases[] = {
768 {{ -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 0
769 {{ 3, 100, 0, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 1
770 {{ 3, 101, 1, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 2
771 {{ 4, 100, 101, 1, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 3
772 {{ 3, 102, 2, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 4
773 {{ 3, 100, 0, 3, 3, 102, 2, 1, 4, 0, 1, 2, 3, -1}}, // 5
774 {{ 4, 101, 102, 2, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 6
775 {{ 3, 100, 101, 3, 3, 101, 2, 3, 3, 101, 102, 2, -1, -1}}, // 7
776 {{ 3, 103, 3, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 8
777 {{ 4, 100, 0, 2, 103, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 9
778 {{ 3, 101, 1, 0, 3, 103, 3, 2, 4, 0, 1, 2, 3, -1}}, // 10
779 {{ 3, 100, 101, 1, 3, 100, 1, 2, 3, 100, 2, 103, -1, -1}}, // 11
780 {{ 4, 102, 103, 3, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 12
781 {{ 3, 100, 0, 103, 3, 0, 1, 103, 3, 1, 102, 103, -1, -1}}, // 13
782 {{ 3, 0, 101, 102, 3, 0, 102, 3, 3, 102, 103, 3, -1, -1}}, // 14
783 {{ 4, 100, 101, 102, 103, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 15
784 };
785
786 static QUAD_CASES quadCasesComplement[] = {
787 {{ -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 0
788 {{ 3, 100, 0, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 1
789 {{ 3, 101, 1, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 2
790 {{ 4, 100, 101, 1, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 3
791 {{ 3, 102, 2, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 4
792 {{ 3, 100, 0, 3, 3, 102, 2, 1, -1, -1, -1, -1, -1, -1}}, // 5
793 {{ 4, 101, 102, 2, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 6
794 {{ 3, 100, 101, 3, 3, 101, 2, 3, 3, 101, 102, 2, -1, -1}}, // 7
795 {{ 3, 103, 3, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 8
796 {{ 4, 100, 0, 2, 103, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 9
797 {{ 3, 101, 1, 0, 3, 103, 3, 2, -1, -1, -1, -1, -1, -1}}, // 10
798 {{ 3, 100, 101, 1, 3, 100, 1, 2, 3, 100, 2, 103, -1, -1}}, // 11
799 {{ 4, 102, 103, 3, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 12
800 {{ 3, 100, 0, 103, 3, 0, 1, 103, 3, 1, 102, 103, -1, -1}}, // 13
801 {{ 3, 0, 101, 102, 3, 0, 102, 3, 3, 102, 103, 3, -1, -1}}, // 14
802 {{ 4, 100, 101, 102, 103, -1, -1, -1, -1, -1, -1, -1, -1, -1}}, // 15
803 };
804
805 //----------------------------------------------------------------------------
806 // Clip this quad using scalar value provided. Like contouring, except
807 // that it cuts the quad to produce other quads and/or triangles.
Clip(double value,vtkDataArray * cellScalars,vtkIncrementalPointLocator * locator,vtkCellArray * polys,vtkPointData * inPd,vtkPointData * outPd,vtkCellData * inCd,vtkIdType cellId,vtkCellData * outCd,int insideOut)808 void vtkQuad::Clip(double value, vtkDataArray *cellScalars,
809 vtkIncrementalPointLocator *locator, vtkCellArray *polys,
810 vtkPointData *inPd, vtkPointData *outPd,
811 vtkCellData *inCd, vtkIdType cellId, vtkCellData *outCd,
812 int insideOut)
813 {
814 static const int CASE_MASK[4] = {1,2,4,8};
815 QUAD_CASES *quadCase;
816 QUAD_EDGE_LIST *edge;
817 int i, j, index, *vert;
818 int e1, e2;
819 int newCellId;
820 vtkIdType pts[4];
821 int vertexId;
822 double t, x1[3], x2[3], x[3], deltaScalar;
823 double scalar0, scalar1, e1Scalar;
824
825 // Build the index into the case table
826 if ( insideOut )
827 {
828 for ( i=0, index = 0; i < 4; i++)
829 {
830 if (cellScalars->GetComponent(i,0) <= value)
831 {
832 index |= CASE_MASK[i];
833 }
834 }
835 // Select case based on the index and get the list of edges for this case
836 quadCase = quadCases + index;
837 }
838 else
839 {
840 for ( i=0, index = 0; i < 4; i++)
841 {
842 if (cellScalars->GetComponent(i,0) > value)
843 {
844 index |= CASE_MASK[i];
845 }
846 }
847 // Select case based on the index and get the list of edges for this case
848 quadCase = quadCasesComplement + index;
849 }
850
851 edge = quadCase->edges;
852
853 // generate each quad
854 for ( ; edge[0] > -1; edge += edge[0]+1 )
855 {
856 for (i=0; i < edge[0]; i++) // insert quad or triangle
857 {
858 // vertex exists, and need not be interpolated
859 if (edge[i+1] >= 100)
860 {
861 vertexId = edge[i+1] - 100;
862 this->Points->GetPoint(vertexId, x);
863 if ( locator->InsertUniquePoint(x, pts[i]) )
864 {
865 outPd->CopyData(inPd,this->PointIds->GetId(vertexId),pts[i]);
866 }
867 }
868
869 else //new vertex, interpolate
870 {
871 vert = edges[edge[i+1]];
872
873 // calculate a preferred interpolation direction
874 scalar0 = cellScalars->GetComponent(vert[0],0);
875 scalar1 = cellScalars->GetComponent(vert[1],0);
876 deltaScalar = scalar1 - scalar0;
877
878 if (deltaScalar > 0)
879 {
880 e1 = vert[0]; e2 = vert[1];
881 e1Scalar = scalar0;
882 }
883 else
884 {
885 e1 = vert[1]; e2 = vert[0];
886 e1Scalar = scalar1;
887 deltaScalar = -deltaScalar;
888 }
889
890 // linear interpolation
891 if (deltaScalar == 0.0)
892 {
893 t = 0.0;
894 }
895 else
896 {
897 t = (value - e1Scalar) / deltaScalar;
898 }
899
900 this->Points->GetPoint(e1, x1);
901 this->Points->GetPoint(e2, x2);
902
903 for (j=0; j<3; j++)
904 {
905 x[j] = x1[j] + t * (x2[j] - x1[j]);
906 }
907
908 if ( locator->InsertUniquePoint(x, pts[i]) )
909 {
910 vtkIdType p1 = this->PointIds->GetId(e1);
911 vtkIdType p2 = this->PointIds->GetId(e2);
912 outPd->InterpolateEdge(inPd,pts[i],p1,p2,t);
913 }
914 }
915 }
916 // check for degenerate output
917 if ( edge[0] == 3 ) //i.e., a triangle
918 {
919 if (pts[0] == pts[1] || pts[0] == pts[2] || pts[1] == pts[2] )
920 {
921 continue;
922 }
923 }
924 else // a quad
925 {
926 if ((pts[0] == pts[3] && pts[1] == pts[2]) ||
927 (pts[0] == pts[1] && pts[3] == pts[2]) )
928 {
929 continue;
930 }
931 }
932
933 newCellId = polys->InsertNextCell(edge[0],pts);
934 outCd->CopyData(inCd,cellId,newCellId);
935 }
936 }
937
938 //----------------------------------------------------------------------------
939 static double vtkQuadCellPCoords[12] = {0.0,0.0,0.0, 1.0,0.0,0.0,
940 1.0,1.0,0.0, 0.0,1.0,0.0};
GetParametricCoords()941 double *vtkQuad::GetParametricCoords()
942 {
943 return vtkQuadCellPCoords;
944 }
945
946 //----------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)947 void vtkQuad::PrintSelf(ostream& os, vtkIndent indent)
948 {
949 this->Superclass::PrintSelf(os,indent);
950
951 os << indent << "Line:\n";
952 this->Line->PrintSelf(os,indent.GetNextIndent());
953 os << indent << "Triangle:\n";
954 this->Triangle->PrintSelf(os,indent.GetNextIndent());
955 }
956