Hexa.cc 13.9 KB
Newer Older
Will Schroeder's avatar
Will Schroeder committed
1
2
/*=========================================================================

Ken Martin's avatar
Ken Martin committed
3
  Program:   Visualization Toolkit
Will Schroeder's avatar
Will Schroeder committed
4
5
6
7
8
  Module:    Hexa.cc
  Language:  C++
  Date:      $Date$
  Version:   $Revision$

Ken Martin's avatar
Ken Martin committed
9
This file is part of the Visualization Toolkit. No part of this file
Will Schroeder's avatar
Will Schroeder committed
10
11
12
13
14
15
16
or its contents may be copied, reproduced or altered in any way
without the express written consent of the authors.

Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen 1993, 1994 

=========================================================================*/
#include "Hexa.hh"
Ken Martin's avatar
Ken Martin committed
17
#include "vtkMath.hh"
Will Schroeder's avatar
Will Schroeder committed
18
19
20
#include "Line.hh"
#include "Quad.hh"
#include "CellArr.hh"
Will Schroeder's avatar
Will Schroeder committed
21

22
23
// Description:
// Deep copy of cell.
Ken Martin's avatar
Ken Martin committed
24
vtkHexahedron::vtkHexahedron(const vtkHexahedron& h)
25
26
27
28
{
  this->Points = h.Points;
  this->PointIds = h.PointIds;
}
Will Schroeder's avatar
Will Schroeder committed
29
30
31

//
//  Method to calculate parametric coordinates in an eight noded
32
//  linear hexahedron element from global coordinates.
Will Schroeder's avatar
Will Schroeder committed
33
34
35
36
//
#define MAX_ITERATION 10
#define CONVERGED 1.e-03

Ken Martin's avatar
Ken Martin committed
37
int vtkHexahedron::EvaluatePosition(float x[3], float closestPoint[3],
Will Schroeder's avatar
Will Schroeder committed
38
                                   int& subId, float pcoords[3], 
Will Schroeder's avatar
Will Schroeder committed
39
                                   float& dist2, float weights[MAX_CELL_SIZE])
Will Schroeder's avatar
Will Schroeder committed
40
41
42
43
{
  int iteration, converged;
  float  params[3];
  float  fcol[3], rcol[3], scol[3], tcol[3];
Will Schroeder's avatar
Will Schroeder committed
44
  int i, j;
Will Schroeder's avatar
Will Schroeder committed
45
  float  d, *pt;
Will Schroeder's avatar
Will Schroeder committed
46
  float derivs[24];
Ken Martin's avatar
Ken Martin committed
47
  static vtkMath math;
Will Schroeder's avatar
Will Schroeder committed
48
49
50
51
//
//  set initial position for Newton's method
//
  subId = 0;
52
  pcoords[0] = pcoords[1] = pcoords[2] = params[0] = params[1] = params[2] = 0.5;
Will Schroeder's avatar
Will Schroeder committed
53
54
55
56
57
58
59
//
//  enter iteration loop
///
  for (iteration=converged=0; !converged && (iteration < MAX_ITERATION);
  iteration++) 
    {
//
Will Schroeder's avatar
Will Schroeder committed
60
//  calculate element interpolation functions and derivatives
Will Schroeder's avatar
Will Schroeder committed
61
//
Will Schroeder's avatar
Will Schroeder committed
62
63
    this->InterpolationFunctions(pcoords, weights);
    this->InterpolationDerivs(pcoords, derivs);
Will Schroeder's avatar
Will Schroeder committed
64
65
66
67
68
69
70
71
72
73
74
75
//
//  calculate newton functions
//
    for (i=0; i<3; i++) 
      {
      fcol[i] = rcol[i] = scol[i] = tcol[i] = 0.0;
      }
    for (i=0; i<8; i++)
      {
      pt = this->Points.GetPoint(i);
      for (j=0; j<3; j++)
        {
Will Schroeder's avatar
Will Schroeder committed
76
        fcol[j] += pt[j] * weights[i];
Will Schroeder's avatar
Will Schroeder committed
77
78
79
80
81
82
83
84
        rcol[j] += pt[j] * derivs[i];
        scol[j] += pt[j] * derivs[i+8];
        tcol[j] += pt[j] * derivs[i+16];
        }
      }

    for (i=0; i<3; i++) fcol[i] -= x[i];
//
85
//  compute determinants and generate improvements
Will Schroeder's avatar
Will Schroeder committed
86
//
87
    if ( (d=math.Determinant3x3(rcol,scol,tcol)) == 0.0 ) return -1;
Will Schroeder's avatar
Will Schroeder committed
88

89
90
91
    pcoords[0] = params[0] - math.Determinant3x3 (fcol,scol,tcol) / d;
    pcoords[1] = params[1] - math.Determinant3x3 (rcol,fcol,tcol) / d;
    pcoords[2] = params[2] - math.Determinant3x3 (rcol,scol,fcol) / d;
Will Schroeder's avatar
Will Schroeder committed
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
//
//  check for convergence
//
    if ( ((fabs(pcoords[0]-params[0])) < CONVERGED) &&
    ((fabs(pcoords[1]-params[1])) < CONVERGED) &&
    ((fabs(pcoords[2]-params[2])) < CONVERGED) )
      {
      converged = 1;
      }
//
//  if not converged, repeat
//
    else 
      {
      params[0] = pcoords[0];
      params[1] = pcoords[1];
      params[2] = pcoords[2];
      }
    }
//
//  if not converged, set the parametric coordinates to arbitrary values
//  outside of element
//
115
116
117
118
119
120
121
  if ( !converged ) return -1;

  this->InterpolationFunctions(pcoords, weights);

  if ( pcoords[0] >= -0.001 && pcoords[0] <= 1.001 &&
  pcoords[1] >= -0.001 && pcoords[1] <= 1.001 &&
  pcoords[2] >= -0.001 && pcoords[2] <= 1.001 )
Will Schroeder's avatar
Will Schroeder committed
122
    {
123
124
125
    closestPoint[0] = x[0]; closestPoint[1] = x[1]; closestPoint[2] = x[2];
    dist2 = 0.0; //inside hexahedron
    return 1;
Will Schroeder's avatar
Will Schroeder committed
126
127
128
    }
  else
    {
129
130
    float pc[3], w[8];
    for (i=0; i<3; i++) //only approximate, not really true for warped hexa
Will Schroeder's avatar
Will Schroeder committed
131
      {
132
133
134
      if (pcoords[i] < 0.0) pc[i] = 0.0;
      else if (pcoords[i] > 1.0) pc[i] = 1.0;
      else pc[i] = pcoords[i];
Will Schroeder's avatar
Will Schroeder committed
135
      }
136
137
138
    this->EvaluateLocation(subId, pc, closestPoint, w);
    dist2 = math.Distance2BetweenPoints(closestPoint,x);
    return 0;
Will Schroeder's avatar
Will Schroeder committed
139
140
141
    }
}
//
Will Schroeder's avatar
Will Schroeder committed
142
// Compute iso-parametrix interpolation functions
Will Schroeder's avatar
Will Schroeder committed
143
//
Ken Martin's avatar
Ken Martin committed
144
void vtkHexahedron::InterpolationFunctions(float pcoords[3], float sf[8])
Will Schroeder's avatar
Will Schroeder committed
145
{
146
  double rm, sm, tm;
Will Schroeder's avatar
Will Schroeder committed
147
148
149
150

  rm = 1. - pcoords[0];
  sm = 1. - pcoords[1];
  tm = 1. - pcoords[2];
151
152
153
154
155
156
157
158
159

  sf[0] = rm*sm*tm;
  sf[1] = pcoords[0]*sm*tm;
  sf[2] = pcoords[0]*pcoords[1]*tm;
  sf[3] = rm*pcoords[1]*tm;
  sf[4] = rm*sm*pcoords[2];
  sf[5] = pcoords[0]*sm*pcoords[2];
  sf[6] = pcoords[0]*pcoords[1]*pcoords[2];
  sf[7] = rm*pcoords[1]*pcoords[2];
Will Schroeder's avatar
Will Schroeder committed
160
161
}

Ken Martin's avatar
Ken Martin committed
162
void vtkHexahedron::InterpolationDerivs(float pcoords[3], float derivs[24])
Will Schroeder's avatar
Will Schroeder committed
163
{
164
  double rm, sm, tm;
Will Schroeder's avatar
Will Schroeder committed
165
166
167
168

  rm = 1. - pcoords[0];
  sm = 1. - pcoords[1];
  tm = 1. - pcoords[2];
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193

  derivs[0] = -sm*tm;
  derivs[1] = sm*tm;
  derivs[2] = pcoords[1]*tm;
  derivs[3] = -pcoords[1]*tm;
  derivs[4] = -sm*pcoords[2];
  derivs[5] = sm*pcoords[2];
  derivs[6] = pcoords[1]*pcoords[2];
  derivs[7] = -pcoords[1]*pcoords[2];
  derivs[8] = -rm*tm;
  derivs[9] = -pcoords[0]*tm;
  derivs[10] = pcoords[0]*tm;
  derivs[11] = rm*tm;
  derivs[12] = -rm*pcoords[2];
  derivs[13] = -pcoords[0]*pcoords[2];
  derivs[14] = pcoords[0]*pcoords[2];
  derivs[15] = rm*pcoords[2];
  derivs[16] = -rm*sm;
  derivs[17] = -pcoords[0]*sm;
  derivs[18] = -pcoords[0]*pcoords[1];
  derivs[19] = -rm*pcoords[1];
  derivs[20] = rm*sm;
  derivs[21] = pcoords[0]*sm;
  derivs[22] = pcoords[0]*pcoords[1];
  derivs[23] = rm*pcoords[1];
Will Schroeder's avatar
Will Schroeder committed
194
195
}

Ken Martin's avatar
Ken Martin committed
196
void vtkHexahedron::EvaluateLocation(int& subId, float pcoords[3], float x[3],
Will Schroeder's avatar
Will Schroeder committed
197
                                    float weights[MAX_CELL_SIZE])
Will Schroeder's avatar
Will Schroeder committed
198
199
{
  int i, j;
Will Schroeder's avatar
Will Schroeder committed
200
  float *pt;
Will Schroeder's avatar
Will Schroeder committed
201

Will Schroeder's avatar
Will Schroeder committed
202
  this->InterpolationFunctions(pcoords, weights);
Will Schroeder's avatar
Will Schroeder committed
203
204
205
206
207
208
209

  x[0] = x[1] = x[2] = 0.0;
  for (i=0; i<8; i++)
    {
    pt = this->Points.GetPoint(i);
    for (j=0; j<3; j++)
      {
Will Schroeder's avatar
Will Schroeder committed
210
      x[j] += pt[j] * weights[i];
Will Schroeder's avatar
Will Schroeder committed
211
212
213
      }
    }
}
214

Ken Martin's avatar
Ken Martin committed
215
int vtkHexahedron::CellBoundary(int subId, float pcoords[3], vtkIdList& pts)
Will Schroeder's avatar
Will Schroeder committed
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
{
  float t1=pcoords[0]-pcoords[1];
  float t2=1.0-pcoords[0]-pcoords[1];
  float t3=pcoords[1]-pcoords[2];
  float t4=1.0-pcoords[1]-pcoords[2];
  float t5=pcoords[2]-pcoords[0];
  float t6=1.0-pcoords[2]-pcoords[0];

  pts.Reset();

  // compare against six planes in parametric space that divide element
  // into six pieces.
  if ( t3 >= 0.0 && t4 >= 0.0 && t5 < 0.0 && t6 >= 0.0 )
    {
    pts.SetId(0,this->PointIds.GetId(0));
    pts.SetId(1,this->PointIds.GetId(1));
    pts.SetId(2,this->PointIds.GetId(2));
    pts.SetId(3,this->PointIds.GetId(3));
    }

  else if ( t1 >= 0.0 && t2 < 0.0 && t5 < 0.0 && t6 < 0.0 )
    {
    pts.SetId(0,this->PointIds.GetId(1));
    pts.SetId(1,this->PointIds.GetId(2));
    pts.SetId(2,this->PointIds.GetId(6));
    pts.SetId(3,this->PointIds.GetId(5));
    }

  else if ( t1 >= 0.0 && t2 >= 0.0 && t3 < 0.0 && t4 >= 0.0 )
    {
    pts.SetId(0,this->PointIds.GetId(0));
    pts.SetId(1,this->PointIds.GetId(1));
    pts.SetId(2,this->PointIds.GetId(5));
    pts.SetId(3,this->PointIds.GetId(4));
    }

  else if ( t3 < 0.0 && t4 < 0.0 && t5 >= 0.0 && t6 < 0.0 )
    {
    pts.SetId(0,this->PointIds.GetId(4));
    pts.SetId(1,this->PointIds.GetId(5));
    pts.SetId(2,this->PointIds.GetId(6));
    pts.SetId(3,this->PointIds.GetId(7));
    }

  else if ( t1 < 0.0 && t2 >= 0.0 && t5 >= 0.0 && t6 >= 0.0 )
    {
    pts.SetId(0,this->PointIds.GetId(0));
    pts.SetId(1,this->PointIds.GetId(4));
    pts.SetId(2,this->PointIds.GetId(7));
    pts.SetId(3,this->PointIds.GetId(3));
    }

  else // if ( t1 < 0.0 && t2 < 0.0 && t3 >= 0.0 && t6 < 0.0 )
    {
    pts.SetId(0,this->PointIds.GetId(2));
    pts.SetId(1,this->PointIds.GetId(3));
    pts.SetId(2,this->PointIds.GetId(7));
    pts.SetId(3,this->PointIds.GetId(6));
    }


  if ( pcoords[0] < 0.0 || pcoords[0] > 1.0 ||
  pcoords[1] < 0.0 || pcoords[1] > 1.0 || pcoords[2] < 0.0 || pcoords[2] > 1.0 )
    return 0;
  else
    return 1;
}

Will Schroeder's avatar
Will Schroeder committed
284
285
286
287
288
289
static int edges[12][2] = { {0,1}, {1,2}, {2,3}, {3,0},
                            {4,5}, {5,6}, {6,7}, {7,4},
                            {0,4}, {1,5}, {3,7}, {2,6}};
static int faces[6][4] = { {0,4,7,3}, {1,2,6,5},
                           {0,1,5,4}, {3,7,6,2},
                           {0,3,2,1}, {4,5,6,7} };
290
291
292
293
294
//
// Marching cubes case table
//
#include "MC_Cases.h"

Ken Martin's avatar
Ken Martin committed
295
296
297
298
void vtkHexahedron::Contour(float value, vtkFloatScalars *cellScalars, 
                      vtkFloatPoints *points,
                      vtkCellArray *verts, vtkCellArray *lines, 
                      vtkCellArray *polys, vtkFloatScalars *scalars)
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
{
  static int CASE_MASK[8] = {1,2,4,8,16,32,64,128};
  TRIANGLE_CASES *triCase;
  EDGE_LIST  *edge;
  int i, j, index, *vert;
  int pts[3];
  float t, *x1, *x2, x[3];

  // Build the case table
  for ( i=0, index = 0; i < 8; i++)
      if (cellScalars->GetScalar(i) >= value)
          index |= CASE_MASK[i];

  triCase = triCases + index;
  edge = triCase->edges;

  for ( ; edge[0] > -1; edge += 3 )
    {
    for (i=0; i<3; i++) // insert triangle
      {
      vert = edges[edge[i]];
      t = (value - cellScalars->GetScalar(vert[0])) /
          (cellScalars->GetScalar(vert[1]) - cellScalars->GetScalar(vert[0]));
      x1 = this->Points.GetPoint(vert[0]);
      x2 = this->Points.GetPoint(vert[1]);
      for (j=0; j<3; j++) x[j] = x1[j] + t * (x2[j] - x1[j]);
      pts[i] = points->InsertNextPoint(x);
      scalars->InsertNextScalar(value);
      }
    polys->InsertNextCell(3,pts);
    }
}
Will Schroeder's avatar
Will Schroeder committed
331

Ken Martin's avatar
Ken Martin committed
332
vtkCell *vtkHexahedron::GetEdge(int edgeId)
Will Schroeder's avatar
Will Schroeder committed
333
334
{
  int *verts;
Ken Martin's avatar
Ken Martin committed
335
  static vtkLine line;
Will Schroeder's avatar
Will Schroeder committed
336
337
338
339
340
341
342
343
344
345
346
347
348
349

  verts = edges[edgeId];

  // load point id's
  line.PointIds.SetId(0,this->PointIds.GetId(verts[0]));
  line.PointIds.SetId(1,this->PointIds.GetId(verts[1]));

  // load coordinates
  line.Points.SetPoint(0,this->Points.GetPoint(verts[0]));
  line.Points.SetPoint(1,this->Points.GetPoint(verts[1]));

  return &line;
}

Ken Martin's avatar
Ken Martin committed
350
vtkCell *vtkHexahedron::GetFace(int faceId)
Will Schroeder's avatar
Will Schroeder committed
351
352
{
  int *verts, i;
Ken Martin's avatar
Ken Martin committed
353
  static vtkQuad theQuad; // using "quad" bothers IBM xlc compiler!
Will Schroeder's avatar
Will Schroeder committed
354
355
356
357
358

  verts = faces[faceId];

  for (i=0; i<4; i++)
    {
359
360
    theQuad.PointIds.SetId(i,this->PointIds.GetId(verts[i]));
    theQuad.Points.SetPoint(i,this->Points.GetPoint(verts[i]));
Will Schroeder's avatar
Will Schroeder committed
361
362
    }

363
  return &theQuad;
Will Schroeder's avatar
Will Schroeder committed
364
}
365
366
367
// 
// Intersect hexa faces against line. Each hexa face is a quadrilateral.
//
Ken Martin's avatar
Ken Martin committed
368
int vtkHexahedron::IntersectWithLine(float p1[3], float p2[3], float tol,
369
370
371
372
373
374
375
376
                                    float &t, float x[3], float pcoords[3],
                                    int& subId)
{
  int intersection=0;
  float *pt1, *pt2, *pt3, *pt4;
  float tTemp;
  float pc[3], xTemp[3];
  int faceNum;
Ken Martin's avatar
Ken Martin committed
377
  static vtkQuad theQuad; // using "quad" bothers IBM xlc compiler!
378
379
380
381
382
383
384
385
386

  t = LARGE_FLOAT;
  for (faceNum=0; faceNum<6; faceNum++)
    {
    pt1 = this->Points.GetPoint(faces[faceNum][0]);
    pt2 = this->Points.GetPoint(faces[faceNum][1]);
    pt3 = this->Points.GetPoint(faces[faceNum][2]);
    pt4 = this->Points.GetPoint(faces[faceNum][3]);

387
388
389
390
    theQuad.Points.SetPoint(0,pt1);
    theQuad.Points.SetPoint(1,pt2);
    theQuad.Points.SetPoint(2,pt3);
    theQuad.Points.SetPoint(3,pt4);
391

392
    if ( theQuad.IntersectWithLine(p1, p2, tol, tTemp, xTemp, pc, subId) )
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
      {
      intersection = 1;
      if ( tTemp < t )
        {
        t = tTemp;
        x[0] = xTemp[0]; x[1] = xTemp[1]; x[2] = xTemp[2]; 
        switch (faceNum)
          {
          case 0:
            pcoords[0] = 0.0; pcoords[0] = pc[0]; pcoords[1] = 0.0;
            break;

          case 1:
            pcoords[0] = 1.0; pcoords[0] = pc[0]; pcoords[1] = 0.0;
            break;

          case 2:
            pcoords[0] = pc[0]; pcoords[1] = 0.0; pcoords[2] = pc[1];
            break;

          case 3:
            pcoords[0] = pc[0]; pcoords[1] = 1.0; pcoords[2] = pc[1];
            break;

          case 4:
            pcoords[0] = pc[0]; pcoords[1] = pc[1]; pcoords[2] = 0.0;
            break;

          case 5:
            pcoords[0] = pc[0]; pcoords[1] = pc[1]; pcoords[2] = 1.0;
            break;
          }
        }
      }
    }
  return intersection;
}
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501

int vtkHexahedron::Triangulate(int index, vtkFloatPoints &pts)
{
  pts.Reset();
//
// Create five tetrahedron. Triangulation varies depending upon index. This
// is necessary to insure compatible voxel triangulations.
//
  if ( index % 2 )
    {
    pts.InsertNextPoint(this->Points.GetPoint(0));
    pts.InsertNextPoint(this->Points.GetPoint(1));
    pts.InsertNextPoint(this->Points.GetPoint(4));
    pts.InsertNextPoint(this->Points.GetPoint(3));

    pts.InsertNextPoint(this->Points.GetPoint(1));
    pts.InsertNextPoint(this->Points.GetPoint(4));
    pts.InsertNextPoint(this->Points.GetPoint(7));
    pts.InsertNextPoint(this->Points.GetPoint(5));

    pts.InsertNextPoint(this->Points.GetPoint(1));
    pts.InsertNextPoint(this->Points.GetPoint(4));
    pts.InsertNextPoint(this->Points.GetPoint(3));
    pts.InsertNextPoint(this->Points.GetPoint(6));

    pts.InsertNextPoint(this->Points.GetPoint(1));
    pts.InsertNextPoint(this->Points.GetPoint(3));
    pts.InsertNextPoint(this->Points.GetPoint(2));
    pts.InsertNextPoint(this->Points.GetPoint(6));

    pts.InsertNextPoint(this->Points.GetPoint(3));
    pts.InsertNextPoint(this->Points.GetPoint(6));
    pts.InsertNextPoint(this->Points.GetPoint(4));
    pts.InsertNextPoint(this->Points.GetPoint(7));
    }
  else
    {
    pts.InsertNextPoint(this->Points.GetPoint(2));
    pts.InsertNextPoint(this->Points.GetPoint(1));
    pts.InsertNextPoint(this->Points.GetPoint(0));
    pts.InsertNextPoint(this->Points.GetPoint(5));

    pts.InsertNextPoint(this->Points.GetPoint(0));
    pts.InsertNextPoint(this->Points.GetPoint(2));
    pts.InsertNextPoint(this->Points.GetPoint(7));
    pts.InsertNextPoint(this->Points.GetPoint(3));

    pts.InsertNextPoint(this->Points.GetPoint(2));
    pts.InsertNextPoint(this->Points.GetPoint(5));
    pts.InsertNextPoint(this->Points.GetPoint(7));
    pts.InsertNextPoint(this->Points.GetPoint(6));

    pts.InsertNextPoint(this->Points.GetPoint(0));
    pts.InsertNextPoint(this->Points.GetPoint(7));
    pts.InsertNextPoint(this->Points.GetPoint(5));
    pts.InsertNextPoint(this->Points.GetPoint(4));

    pts.InsertNextPoint(this->Points.GetPoint(1));
    pts.InsertNextPoint(this->Points.GetPoint(2));
    pts.InsertNextPoint(this->Points.GetPoint(5));
    pts.InsertNextPoint(this->Points.GetPoint(7));
    }

  return 1;
}

void vtkHexahedron::Derivatives(int subId, float pcoords[3], float *values, 
                                int dim, float *derivs)
{

}