Overture  Version 25
MeshQuality.h
Go to the documentation of this file.
1 #ifndef __MESH_QUALITY__
2 #define __MESH_QUALITY__
3 
4 #include "GenericDataBase.h"
5 #include "OvertureTypes.h"
6 #include "UnstructuredMapping.h"
7 #include "CompositeGrid.h"
8 #include "CompositeGridFunction.h"
9 #include "GL_GraphicsInterface.h"
10 #include "ArraySimple.h"
11 #include "InterpolatePoints.h"
12 
13 class MetricEvaluator {
14 public:
15  virtual ArraySimpleFixed<real,3,3,1,1> computeMetric(const ArraySimpleFixed<real,3,1,1,1> &x)=0;
16  virtual ArraySimpleFixed<real,3,3,1,1> computeMetric(int id)=0;
17  virtual ArraySimpleFixed<real,3,3,1,1> computeMetric(void *entity)=0;
18 };
19 
20 class IdentityMetricEvaluator : public MetricEvaluator {
21 
22 public:
23  IdentityMetricEvaluator() { }
24  virtual ~IdentityMetricEvaluator() { }
25 
26  virtual ArraySimpleFixed<real,3,3,1,1> computeMetric(const ArraySimpleFixed<real,3,1,1,1> &x);
27  virtual ArraySimpleFixed<real,3,3,1,1> computeMetric(int id);
28  virtual ArraySimpleFixed<real,3,3,1,1> computeMetric(void *entity);
29 
30 };
31 
32 class MetricCGFunctionEvaluator : public MetricEvaluator {
33 
34 public:
35  MetricCGFunctionEvaluator(RealCompositeGridFunction *cgf_=0) : cgf(cgf_){ }
36  virtual ~MetricCGFunctionEvaluator() { }
37 
38  virtual ArraySimpleFixed<real,3,3,1,1> computeMetric(const ArraySimpleFixed<real,3,1,1,1> &x);
39  virtual ArraySimpleFixed<real,3,3,1,1> computeMetric(int id);
40  virtual ArraySimpleFixed<real,3,3,1,1> computeMetric(void *entity);
41 
42 private:
43  RealCompositeGridFunction *cgf;
44 
45 };
46 
47 
48 enum MeshQualityMetric {
49  volumeMetric=0,
50  shapeMetric,
51  volumeShapeMetric,
52  numberOfQualityMetrics
53 };
54 
55 class MeshQualityMetrics
56 {
57 
58 public:
59 
60  enum JacobianProperties {
61  normSquared=0,
62  determinant,
63  conditionNum,
64  numberOfProperties
65  };
66 
67  MeshQualityMetrics();
68  MeshQualityMetrics(UnstructuredMapping &umap_);
69 
70  const realArray &computeMetric(MeshQualityMetric metric);
71 
72  void setReferenceTransformation(MetricEvaluator *rt);
73 
74  void plot(GL_GraphicsInterface &gi);
75 
76  void outputHistogram(aString fileName="meshQuality.dat");
77 
78  const realArray & getJacobianProperties() { return jacobianProperties; }
79  const realArray & getJacobians() { return jacobians; }
80 
81  // 2d triangle
82  inline ArraySimpleFixed<real,2,2,1,1> computeJacobian(const ArraySimpleFixed<real, 2,1,1,1> &x0,
83  const ArraySimpleFixed<real,2,1,1,1> &xm,
84  const ArraySimpleFixed<real,2,1,1,1> &xmp1,
85  const ArraySimpleFixed<real,2,2,1,1> &T);
86 
87  // 2d quadrilateral
88  inline ArraySimpleFixed<real,2,2,1,1> computeJacobian(const ArraySimpleFixed<real,2,1,1,1> &x0,
89  const ArraySimpleFixed<real,2,1,1,1> &x1,
90  const ArraySimpleFixed<real,2,1,1,1> &x2,
91  const ArraySimpleFixed<real,2,1,1,1> &x3,
92  const ArraySimpleFixed<real,2,2,1,1> &T);
93 
94  // 3d tetrahedron
95  inline ArraySimpleFixed<real,3,3,1,1> computeJacobian(const ArraySimpleFixed<real,3,1,1,1> &x0,
96  const ArraySimpleFixed<real,3,1,1,1> &x1,
97  const ArraySimpleFixed<real,3,1,1,1> &x2,
98  const ArraySimpleFixed<real,3,1,1,1> &x3,
99  const ArraySimpleFixed<real,3,3,1,1> &T);
100 
101  // 3d pyramid
102  inline ArraySimpleFixed<real,3,3,1,1> computeJacobian(const ArraySimpleFixed<real,3,1,1,1> &x0,
103  const ArraySimpleFixed<real,3,1,1,1> &x1,
104  const ArraySimpleFixed<real,3,1,1,1> &x2,
105  const ArraySimpleFixed<real,3,1,1,1> &x3,
106  const ArraySimpleFixed<real,3,1,1,1> &x4,
107  const ArraySimpleFixed<real,3,3,1,1> &T);
108 
109  // 3d hexahedron
110  inline ArraySimpleFixed<real,3,3,1,1> computeJacobian(const ArraySimpleFixed<real,3,1,1,1> &x0,
111  const ArraySimpleFixed<real,3,1,1,1> &x1,
112  const ArraySimpleFixed<real,3,1,1,1> &x2,
113  const ArraySimpleFixed<real,3,1,1,1> &x3,
114  const ArraySimpleFixed<real,3,1,1,1> &x4,
115  const ArraySimpleFixed<real,3,1,1,1> &x5,
116  const ArraySimpleFixed<real,3,1,1,1> &x6,
117  const ArraySimpleFixed<real,3,1,1,1> &x7,
118  const ArraySimpleFixed<real,3,3,1,1> &T);
119 
120  inline ArraySimpleFixed<real,2,2,1,1> computeWeight(const ArraySimpleFixed<real,2,1,1,1> &xc,
121  UnstructuredMapping::ElementType et);
122 
123  inline ArraySimpleFixed<real,3,3,1,1> computeWeight(const ArraySimpleFixed<real,3,1,1,1> &xc,
124  UnstructuredMapping::ElementType et);
125 
126  inline void computeJacobianProperties(real &N2,
127  real &det,
128  real &K, const ArraySimpleFixed<real,2,2,1,1> &J) const;
129 
130  inline void computeJacobianProperties(real &N2,
131  real &det,
132  real &K, const ArraySimpleFixed<real,3,3,1,1> &J) const;
133 
134  inline real jacobianNodeDerivative(UnstructuredMapping::ElementType, int n, int c) const;
135 
136 
137 protected:
138  void computeJacobianProperties();
139 
140 private:
141  UnstructuredMapping *umap;
142  MetricEvaluator *referenceTransformation;
143 
144  realArray jacobianProperties;
145  realArray jacobians;
146  realArray metrics[int(numberOfQualityMetrics)];
147 
148 };
149 
150 
151 #if 0
152 // stupid sun compiler can't understand function templates with non-type parameters
153 template<int DIM>
154 void
155 interpolateFromControlFunction( const ArraySimpleFixed<real,DIM,1,1,1> &midPt,
156  ArraySimpleFixed<real,DIM,DIM,1,1> &T,
157  RealCompositeGridFunction &controlFunction )
158 {
159 
160  for ( int a1=0; a1<DIM; a1++ )
161  for ( int a2=0; a2<DIM; a2++ )
162  T(a1,a2) = 0.;
163 
164  Index Iv[3], &I1=Iv[0], &I2=Iv[1], &I3=Iv[2];
165  CompositeGrid &controlGrid = *controlFunction.getCompositeGrid();
166  getIndex(controlGrid[0].gridIndexRange(),I1,I2,I3);
167 
168  RealArray controlBounds(2,3);
169  realArray & vertices = controlGrid[0].vertex();
170 
171  real dx[3];
172  int ii[3];
173  ii[2] = I3.getBase();
174  real dxa[3];
175  dx[2] = dxa[2] = 0.0;
176  int axis;
177 
178  // compute indices into the control grid
179  for ( axis=0; axis<DIM; axis++ )
180  {
181  controlBounds(0,axis) = vertices(I1.getBase(), I2.getBase(), I3.getBase(),axis);
182  controlBounds(1,axis) = vertices(I1.getBound(), I2.getBound(), I3.getBound(), axis);
183  dx[axis] = (controlBounds(1,axis) - controlBounds(0,axis))/real(Iv[axis].getLength()-1);
184  ii[axis] = int( (midPt(axis)-controlBounds(0,axis))/dx[axis] );
185  }
186 
187  // compute linear interpolation coefficients
188  for ( axis=0; axis<DIM; axis++ )
189  dxa[axis] = (midPt(axis) - vertices(ii[0],ii[1],ii[2],axis))/dx[axis];
190 
191  // compute the interpolation
192  if ( DIM==2 )
193  {
194  for ( int ti=0; ti<DIM; ti++ )
195  for ( int tj=0; tj<DIM; tj++ )
196  T(ti,tj) = (( 1.0 -dxa[1] ) * ( ( 1.0-dxa[0] ) * controlFunction[0](ii[0],ii[1],ii[2],ti,tj) +
197  ( dxa[0] ) * controlFunction[0](ii[0]+1,ii[1],ii[2],ti,tj) ) +
198  ( dxa[1] ) * ( ( 1.0-dxa[0] ) * controlFunction[0](ii[0],ii[1]+1,ii[2],ti,tj) +
199  ( dxa[0] ) * controlFunction[0](ii[0]+1, ii[1]+1,ii[2],ti,tj) ) );
200  }
201  else
202  {
203  for ( int ti=0; ti<DIM; ti++ )
204  for ( int tj=0; tj<DIM; tj++ )
205  T(ti,tj) = ( (1.0-dxa[2])*( (( 1.0 -dxa[1] ) * ( ( 1.0-dxa[0] ) * controlFunction[0](ii[0],ii[1],ii[2],ti,tj) +
206  ( dxa[0] ) * controlFunction[0](ii[0]+1,ii[1],ii[2],ti,tj) ) +
207  ( dxa[1] ) * ( ( 1.0-dxa[0] ) * controlFunction[0](ii[0],ii[1]+1,ii[2],ti,tj) +
208  ( dxa[0] ) * controlFunction[0](ii[0]+1, ii[1]+1,ii[2],ti,tj) ) ) ) +
209  ( dxa[2])*( (( 1.0 -dxa[1] ) * ( ( 1.0-dxa[0] ) * controlFunction[0](ii[0],ii[1],ii[2]+1,ti,tj) +
210  ( dxa[0] ) * controlFunction[0](ii[0]+1,ii[1],ii[2]+1,ti,tj) ) +
211  ( dxa[1] ) * ( ( 1.0-dxa[0] ) * controlFunction[0](ii[0],ii[1]+1,ii[2]+1,ti,tj) +
212  ( dxa[0] ) * controlFunction[0](ii[0]+1, ii[1]+1,ii[2]+1,ti,tj) ) ) ));
213 
214  }
215 
216 }
217 #else
218 inline void
219 interpolateFromControlFunction( const ArraySimpleFixed<real,2,1,1,1> &midPt,
220  ArraySimpleFixed<real,2,2,1,1> &T,
221  RealCompositeGridFunction &controlFunction )
222 {
223 
224  int DIM=2;
225  for ( int a1=0; a1<DIM; a1++ )
226  for ( int a2=0; a2<DIM; a2++ )
227  T(a1,a2) = 0.;
228 
229  Index Iv[3], &I1=Iv[0], &I2=Iv[1], &I3=Iv[2];
230  CompositeGrid &controlGrid = *controlFunction.getCompositeGrid();
231  getIndex(controlGrid[0].gridIndexRange(),I1,I2,I3);
232 
233  RealArray controlBounds(2,3);
234  realArray & vertices = controlGrid[0].vertex();
235 
236  real dx[3];
237  int ii[3];
238  ii[2] = I3.getBase();
239  real dxa[3];
240  dx[2] = dxa[2] = 0.0;
241  int axis;
242 
243  // compute indices into the control grid
244  for ( axis=0; axis<DIM; axis++ )
245  {
246  controlBounds(0,axis) = vertices(I1.getBase(), I2.getBase(), I3.getBase(),axis);
247  controlBounds(1,axis) = vertices(I1.getBound(), I2.getBound(), I3.getBound(), axis);
248  dx[axis] = (controlBounds(1,axis) - controlBounds(0,axis))/real(Iv[axis].getLength()-1);
249  ii[axis] = min(Iv[axis].getBound(),max(Iv[axis].getBase(),int( (midPt(axis)-controlBounds(0,axis))/dx[axis] ) ));
250  }
251 
252  // compute linear interpolation coefficients
253  for ( axis=0; axis<DIM; axis++ )
254  dxa[axis] = (midPt(axis) - vertices(ii[0],ii[1],ii[2],axis))/dx[axis];
255 
256  for ( int ti=0; ti<DIM; ti++ )
257  for ( int tj=0; tj<DIM; tj++ )
258  T(ti,tj) = (( 1.0 -dxa[1] ) * ( ( 1.0-dxa[0] ) * controlFunction[0](ii[0],ii[1],ii[2],ti,tj) +
259  ( dxa[0] ) * controlFunction[0](ii[0]+1,ii[1],ii[2],ti,tj) ) +
260  ( dxa[1] ) * ( ( 1.0-dxa[0] ) * controlFunction[0](ii[0],ii[1]+1,ii[2],ti,tj) +
261  ( dxa[0] ) * controlFunction[0](ii[0]+1, ii[1]+1,ii[2],ti,tj) ) );
262 }
263 
264 inline void
265 interpolateFromControlFunction( const ArraySimpleFixed<real,3,1,1,1> &midPt,
266  ArraySimpleFixed<real,3,3,1,1> &T,
267  RealCompositeGridFunction &controlFunction )
268 {
269 
270  int DIM=3;
271  for ( int a1=0; a1<DIM; a1++ )
272  for ( int a2=0; a2<DIM; a2++ )
273  T(a1,a2) = 0.;
274 
275  Index Iv[3], &I1=Iv[0], &I2=Iv[1], &I3=Iv[2];
276  CompositeGrid &controlGrid = *controlFunction.getCompositeGrid();
277  getIndex(controlGrid[0].gridIndexRange(),I1,I2,I3);
278 
279  RealArray controlBounds(2,3);
280  realArray & vertices = controlGrid[0].vertex();
281 
282  real dx[3];
283  int ii[3];
284  ii[2] = I3.getBase();
285  real dxa[3];
286  dx[2] = dxa[2] = 0.0;
287  int axis;
288 
289  // compute indices into the control grid
290  for ( axis=0; axis<DIM; axis++ )
291  {
292  controlBounds(0,axis) = vertices(I1.getBase(), I2.getBase(), I3.getBase(),axis);
293  controlBounds(1,axis) = vertices(I1.getBound(), I2.getBound(), I3.getBound(), axis);
294  dx[axis] = (controlBounds(1,axis) - controlBounds(0,axis))/real(Iv[axis].getLength()-1);
295  // ii[axis] = int( (midPt(axis)-controlBounds(0,axis))/dx[axis] );
296  ii[axis] = min(Iv[axis].getBound(),max(Iv[axis].getBase(),int( (midPt(axis)-controlBounds(0,axis))/dx[axis] ) ));
297  }
298 
299  // compute linear interpolation coefficients
300  for ( axis=0; axis<DIM; axis++ )
301  dxa[axis] = (midPt(axis) - vertices(ii[0],ii[1],ii[2],axis))/dx[axis];
302 
303 
304  for ( int ti=0; ti<DIM; ti++ )
305  for ( int tj=0; tj<DIM; tj++ )
306  T(ti,tj) = ( (1.0-dxa[2])*( (( 1.0 -dxa[1] ) * ( ( 1.0-dxa[0] ) * controlFunction[0](ii[0],ii[1],ii[2],ti,tj) +
307  ( dxa[0] ) * controlFunction[0](ii[0]+1,ii[1],ii[2],ti,tj) ) +
308  ( dxa[1] ) * ( ( 1.0-dxa[0] ) * controlFunction[0](ii[0],ii[1]+1,ii[2],ti,tj) +
309  ( dxa[0] ) * controlFunction[0](ii[0]+1, ii[1]+1,ii[2],ti,tj) ) ) ) +
310  ( dxa[2])*( (( 1.0 -dxa[1] ) * ( ( 1.0-dxa[0] ) * controlFunction[0](ii[0],ii[1],ii[2]+1,ti,tj) +
311  ( dxa[0] ) * controlFunction[0](ii[0]+1,ii[1],ii[2]+1,ti,tj) ) +
312  ( dxa[1] ) * ( ( 1.0-dxa[0] ) * controlFunction[0](ii[0],ii[1]+1,ii[2]+1,ti,tj) +
313  ( dxa[0] ) * controlFunction[0](ii[0]+1, ii[1]+1,ii[2]+1,ti,tj) ) ) ));
314 
315 }
316 #endif
317 
318 
319 inline ArraySimpleFixed<real,2,2,1,1>
320 MeshQualityMetrics::
321 computeJacobian(const ArraySimpleFixed<real, 2,1,1,1> &x0,
322  const ArraySimpleFixed<real,2,1,1,1> &xm,
323  const ArraySimpleFixed<real,2,1,1,1> &xmp1,
324  const ArraySimpleFixed<real,2,2,1,1> &T)
325 {
326 
327  ArraySimpleFixed<real,2,2,1,1> J,Jt;
328 
329  Jt(0,0) = xm[0]-x0[0];
330  Jt(0,1) = xmp1[0]-x0[0];
331 
332  Jt(1,0) = xm[1]-x0[1];
333  Jt(1,1) = xmp1[1]-x0[1];
334 
335  int r,c,cc;
336  for ( r=0; r<2; r++ )
337  for ( c=0; c<2; c++ )
338  J(r,c) = 0;
339 
340  for ( r=0; r<2; r++ )
341  for ( c=0; c<2; c++ )
342  for ( cc=0; cc<2; cc++ )
343  J(r,c) += Jt(r,cc)*T(cc,c);
344 
345  return J;
346 }
347 
348 inline ArraySimpleFixed<real,2,2,1,1>
349 MeshQualityMetrics::
350 computeJacobian(const ArraySimpleFixed<real,2,1,1,1> &x0,
351  const ArraySimpleFixed<real,2,1,1,1> &x1,
352  const ArraySimpleFixed<real,2,1,1,1> &x2,
353  const ArraySimpleFixed<real,2,1,1,1> &x3,
354  const ArraySimpleFixed<real,2,2,1,1> &T)
355 {
356 
357  ArraySimpleFixed<real,2,2,1,1> J,Jt;
358 
359  Jt(0,0) = 0.5*(x1[0]+x2[0]-x3[0]-x0[0]);
360  Jt(0,1) = 0.5*(x3[0]+x2[0]-x1[0]-x0[0]);
361  Jt(1,0) = 0.5*(x1[1]+x2[1]-x3[1]-x0[1]);
362  Jt(1,1) = 0.5*(x3[1]+x2[1]-x1[1]-x0[1]);
363 
364  int r,c,cc;
365  for ( r=0; r<2; r++ )
366  for ( c=0; c<2; c++ )
367  J(r,c) = 0;
368 
369  for ( r=0; r<2; r++ )
370  for ( c=0; c<2; c++ )
371  for ( cc=0; cc<2; cc++ )
372  J(r,c) += Jt(r,cc)*T(cc,c);
373 
374  return J;
375 }
376 
377 inline ArraySimpleFixed<real,3,3,1,1>
378 MeshQualityMetrics::
379 computeJacobian(const ArraySimpleFixed<real,3,1,1,1> &x0,
380  const ArraySimpleFixed<real,3,1,1,1> &x1,
381  const ArraySimpleFixed<real,3,1,1,1> &x2,
382  const ArraySimpleFixed<real,3,1,1,1> &x3,
383  const ArraySimpleFixed<real,3,3,1,1> &T)
384 {
385 
386  ArraySimpleFixed<real,3,3,1,1> J,Jt;
387 
388  Jt(0,0) = x1[0]-x0[0];
389  Jt(0,1) = x2[0]-x0[0];
390  Jt(0,2) = x3[0]-x0[0];
391 
392  Jt(1,0) = x1[1]-x0[1];
393  Jt(1,1) = x2[1]-x0[1];
394  Jt(1,2) = x3[1]-x0[1];
395 
396  Jt(2,0) = x1[2]-x0[2];
397  Jt(2,1) = x2[2]-x0[2];
398  Jt(2,2) = x3[2]-x0[2];
399 
400  int r,c,cc;
401  for ( r=0; r<3; r++ )
402  for ( c=0; c<3; c++ )
403  J(r,c) = 0;
404 
405  for ( r=0; r<3; r++ )
406  for ( c=0; c<3; c++ )
407  for ( cc=0; cc<3; cc++ )
408  J(r,c) += Jt(r,cc)*T(cc,c);
409 
410  return J;
411 }
412 
413 inline ArraySimpleFixed<real,3,3,1,1>
414 MeshQualityMetrics::
415 computeJacobian(const ArraySimpleFixed<real,3,1,1,1> &x0,
416  const ArraySimpleFixed<real,3,1,1,1> &x1,
417  const ArraySimpleFixed<real,3,1,1,1> &x2,
418  const ArraySimpleFixed<real,3,1,1,1> &x3,
419  const ArraySimpleFixed<real,3,1,1,1> &x4,
420  const ArraySimpleFixed<real,3,3,1,1> &T)
421 {
422  ArraySimpleFixed<real,3,3,1,1> J,Jt;
423 
424  Jt(0,0) = -( x0(0)+x3(0)+x4(0)-
425  (x1(0)+x2(0)+x4(0)))/3.;
426  Jt(1,0) = -( x0(1)+x3(1)+x4(1)-
427  (x1(1)+x2(1)+x4(1)))/3.;
428  Jt(2,0) = -( x0(2)+x3(2)+x4(2)-
429  (x1(2)+x2(2)+x4(2)))/3.;
430 
431  Jt(0,1) = -( x0(0)+x1(0)+x4(0)-
432  (x2(0)+x3(0)+x4(0)))/3.;
433  Jt(1,1) = -( x0(1)+x1(1)+x4(1)-
434  (x2(1)+x3(1)+x4(1)))/3.;
435  Jt(2,1) = -( x0(2)+x1(2)+x4(2)-
436  (x2(2)+x3(2)+x4(2)))/3.;
437 
438  Jt(0,2) = -2*( (x0(0)+x1(0)+
439  x2(0)+x3(0))/4.-
440  (x4(0)));
441  Jt(1,2) = -2*( (x0(1)+x1(1)+x2(1)+
442  x3(1))/4.-
443  (x4(1)));
444  Jt(2,2) = -2*( (x0(2)+x1(2)+x2(2)+
445  x3(2))/4.-
446  (x4(2)));
447 
448  int r,c,cc;
449  for ( r=0; r<3; r++ )
450  for ( c=0; c<3; c++ )
451  J(r,c) = 0;
452 
453  for ( r=0; r<3; r++ )
454  for ( c=0; c<3; c++ )
455  for ( cc=0; cc<3; cc++ )
456  J(r,c) += Jt(r,cc)*T(cc,c);
457 
458  return J;
459 }
460 
461 inline ArraySimpleFixed<real,3,3,1,1>
462 MeshQualityMetrics::
463 computeJacobian(const ArraySimpleFixed<real,3,1,1,1> &x0,
464  const ArraySimpleFixed<real,3,1,1,1> &x1,
465  const ArraySimpleFixed<real,3,1,1,1> &x2,
466  const ArraySimpleFixed<real,3,1,1,1> &x3,
467  const ArraySimpleFixed<real,3,1,1,1> &x4,
468  const ArraySimpleFixed<real,3,1,1,1> &x5,
469  const ArraySimpleFixed<real,3,1,1,1> &x6,
470  const ArraySimpleFixed<real,3,1,1,1> &x7,
471  const ArraySimpleFixed<real,3,3,1,1> &T)
472 {
473  ArraySimpleFixed<real,3,3,1,1> J,Jt;
474 
475  Jt(0,0) = -( x0(0)+x3(0)+x4(0)+x7(0)-
476  (x1(0)+x2(0)+x5(0)+x6(0)))/4.;
477  Jt(1,0) = -( x0(1)+x3(1)+x4(1)+x7(1)-
478  (x1(1)+x2(1)+x5(1)+x6(1)))/4.;
479  Jt(2,0) = -( x0(2)+x3(2)+x4(2)+x7(2)-
480  (x1(2)+x2(2)+x5(2)+x6(2)))/4.;
481 
482  Jt(0,1) = -( x0(0)+x1(0)+x4(0)+x5(0)-
483  (x2(0)+x3(0)+x6(0)+x7(0)))/4.;
484  Jt(1,1) = -( x0(1)+x1(1)+x4(1)+x5(1)-
485  (x2(1)+x3(1)+x6(1)+x7(1)))/4.;
486  Jt(2,1) = -( x0(2)+x1(2)+x4(2)+x5(2)-
487  (x2(2)+x3(2)+x6(2)+x7(2)))/4.;
488 
489  Jt(0,2) = -( x0(0)+x1(0)+x2(0)+x3(0)-
490  (x4(0)+x5(0)+x6(0)+x7(0)))/4.;
491  Jt(1,2) = -( x0(1)+x1(1)+x2(1)+x3(1)-
492  (x4(1)+x5(1)+x6(1)+x7(1)))/4.;
493  Jt(2,2) = -( x0(2)+x1(2)+x2(2)+x3(2)-
494  (x4(2)+x5(2)+x6(2)+x7(2)))/4.;
495 
496  int r,c,cc;
497  for ( r=0; r<3; r++ )
498  for ( c=0; c<3; c++ )
499  J(r,c) = 0;
500 
501  for ( r=0; r<3; r++ )
502  for ( c=0; c<3; c++ )
503  for ( cc=0; cc<3; cc++ )
504  J(r,c) += Jt(r,cc)*T(cc,c);
505 
506  return J;
507 }
508 
509 
510 inline void
511 MeshQualityMetrics::
512 computeJacobianProperties(real &N2,
513  real &det,
514  real &K, const ArraySimpleFixed<real,2,2,1,1> &J) const
515 {
516  int r,c;
517  N2 = det = K = 0.;
518 
519  for ( r=0; r<2; r++ )
520  for ( c=0; c<2; c++ )
521  N2 += J(r,c)*J(r,c);
522 
523  det = J(0,0)*J(1,1) - J(1,0)*J(0,1);
524 
525  if ( det > 0.0 )
526  {
527  real normOfInverse = 0.;
528 
529  normOfInverse += J(1,1)*J(1,1);
530  normOfInverse += J(1,0)*J(1,0);
531  normOfInverse += J(0,1)*J(0,1);
532  normOfInverse += J(0,0)*J(0,0);
533 
534  normOfInverse = sqrt(normOfInverse)/det;
535 
536  K = sqrt(N2)*normOfInverse;
537 
538  }
539  else
540  K = REAL_MAX;
541 
542 }
543 
544 inline void
545 MeshQualityMetrics::
546 computeJacobianProperties(real &N2,
547  real &det,
548  real &K, const ArraySimpleFixed<real,3,3,1,1> &J) const
549 {
550  int r,c;
551  N2 = det = K = 0.;
552 
553  for ( r=0; r<3; r++ )
554  for ( c=0; c<3; c++ )
555  N2 += J(r,c)*J(r,c);
556 
557  det =
558  J(0,0)*(J(1,1)*J(2,2)-J(1,2)*J(2,1)) -
559  J(0,1)*(J(1,0)*J(2,2)-J(1,2)*J(2,0)) +
560  J(0,2)*(J(1,0)*J(2,1)-J(1,1)*J(2,0));
561 
562  if ( det>0.0 )
563  {
564  real normOfInverse = 0.;
565 
566  normOfInverse += (J(1,1)*J(2,2)-J(1,2)*J(2,1))*
567  (J(1,1)*J(2,2)-J(1,2)*J(2,1));
568  normOfInverse += (J(1,0)*J(2,2)-J(1,2)*J(2,0))*
569  (J(1,0)*J(2,2)-J(1,2)*J(2,0));
570  normOfInverse += (J(1,0)*J(2,1)-J(1,1)*J(2,0))*
571  (J(1,0)*J(2,1)-J(1,1)*J(2,0));
572  normOfInverse += (J(0,1)*J(2,2)-J(0,2)*J(2,1))*
573  (J(0,1)*J(2,2)-J(0,2)*J(2,1));
574  normOfInverse += (J(0,0)*J(2,2)-J(0,2)*J(2,0))*
575  (J(0,0)*J(2,2)-J(0,2)*J(2,0));
576  normOfInverse += (J(0,0)*J(2,1)-J(0,1)*J(2,0))*
577  (J(0,0)*J(2,1)-J(0,1)*J(2,0));
578  normOfInverse += (J(0,1)*J(1,2)-J(0,2)*J(1,1))*
579  (J(0,1)*J(1,2)-J(0,2)*J(1,1));
580  normOfInverse += (J(0,0)*J(1,2)-J(0,2)*J(1,0))*
581  (J(0,0)*J(1,2)-J(0,2)*J(1,0));
582  normOfInverse += (J(0,0)*J(1,1)-J(0,1)*J(1,0))*
583  (J(0,0)*J(1,1)-J(0,1)*J(1,0));
584 
585  normOfInverse = sqrt(normOfInverse)/det;
586 
587  K = sqrt(N2)*normOfInverse;
588  }
589  else
590  K = REAL_MAX;
591 
592 }
593 
594 inline ArraySimpleFixed<real,2,2,1,1>
595 MeshQualityMetrics::
596 computeWeight(const ArraySimpleFixed<real, 2,1,1,1> &xc_,
597  UnstructuredMapping::ElementType et)
598 {
599  //realArray xc(1,2), T(2,2);
600  ArraySimpleFixed<real,2,2,1,1> T;
601 
602  T(0,0) = T(1,1) = 1.;
603  T(0,1) = T(1,0) = 0.;
604 
605  if ( referenceTransformation!=NULL )
606  {
607  ArraySimpleFixed<real,3,1,1,1> x3d;
608  ArraySimpleFixed<real,3,3,1,1> T3d;
609  x3d=0;
610  T3d=0;
611  x3d[0] = xc_[0];
612  x3d[1] = xc_[1];
613  T3d = referenceTransformation->computeMetric(x3d);
614  T(0,0) = T3d(0,0);
615  T(1,0) = T3d(1,0);
616  T(1,1) = T3d(1,1);
617  T(0,1) = T3d(0,1);
618  }
619 
620  //interpolateFromControlFunction(xc_, T, *referenceTransformation);
621 
622  ArraySimpleFixed<real,2,2,1,1> Winv;
623  if ( et==UnstructuredMapping::triangle )
624  {
625  Winv(0,0) = T(0,0) - T(1,0)/sqrt(3.);
626  Winv(0,1) = T(0,1) - T(1,1)/sqrt(3.);
627  Winv(1,0) = 2*T(1,0)/sqrt(3.);
628  Winv(1,1) = 2*T(1,1)/sqrt(3.);
629  }
630  else
631  for ( int r=0; r<2; r++ )
632  for ( int c=0; c<2; c++ )
633  Winv(r,c) = T(r,c);
634 
635  return Winv;
636 }
637 
638 inline ArraySimpleFixed<real,3,3,1,1>
639 MeshQualityMetrics::
640 computeWeight(const ArraySimpleFixed<real, 3,1,1,1> &xc_,
641  UnstructuredMapping::ElementType et)
642 {
643  ArraySimpleFixed<real,3,3,1,1> T;
644 
645  T(0,0) = T(1,1) = T(2,2) = 1.;
646  T(0,1) = T(1,0) = T(0,2) = T(1,2) = T(2,1) = T(2,0) = 0.;
647 
648  if ( referenceTransformation!=NULL )
649  T = referenceTransformation->computeMetric(xc_);
650  // interpolateFromControlFunction(xc_, T, *referenceTransformation);
651 
652  ArraySimpleFixed<real,3,3,1,1> Winv;
653  if ( et==UnstructuredMapping::tetrahedron )
654  {
655  Winv(0,0) = T(0,0) - T(1,0)/sqrt(3.) - T(2,0)/sqrt(3.)/sqrt(2.);
656  Winv(0,1) = T(0,1) - T(1,1)/sqrt(3.) - T(2,1)/sqrt(3.)/sqrt(2.);
657  Winv(0,2) = T(0,2) - T(1,2)/sqrt(3.) - T(2,2)/sqrt(3.)/sqrt(2.);
658 
659  Winv(1,0) = 2*T(1,0)/sqrt(3.) - T(2,0)/sqrt(3.)/sqrt(2.);
660  Winv(1,1) = 2*T(1,1)/sqrt(3.) - T(2,1)/sqrt(3.)/sqrt(2.);
661  Winv(1,2) = 2*T(1,2)/sqrt(3.) - T(2,2)/sqrt(3.)/sqrt(2.);
662 
663  Winv(2,0) = sqrt(3.)*T(2,0)/sqrt(2.);
664  Winv(2,1) = sqrt(3.)*T(2,1)/sqrt(2.);
665  Winv(2,2) = sqrt(3.)*T(2,2)/sqrt(2.);
666  }
667  else if ( et==UnstructuredMapping::pyramid )
668  {
669  Winv(0,0) = T(0,0)*3./2;
670  Winv(0,1) = T(0,1)*3./2;
671  Winv(0,2) = T(0,2);
672 
673  Winv(1,0) = T(1,0)*3./2;
674  Winv(1,1) = T(1,1)*3./2;
675  Winv(1,2) = T(1,2)/2.;
676 
677  Winv(2,0) = T(2,0)*3./2;
678  Winv(2,1) = T(2,1)*3./2;
679  Winv(2,2) = T(2,2);
680  }
681  else
682  for ( int r=0; r<3; r++ )
683  for ( int c=0; c<3; c++ )
684  Winv(r,c) = T(r,c);
685 
686  return Winv;
687 }
688 
689 inline
690 real
691 MeshQualityMetrics::
692 jacobianNodeDerivative(UnstructuredMapping::ElementType e, int n, int c) const
693 {
694  real res = -1;
695 
696  switch(e) {
697 
698  case UnstructuredMapping::triangle : case UnstructuredMapping::tetrahedron :
699  if ( n==0 )
700  res = -1;
701  else if ( n==(c+1) )
702  res = 1;
703  else
704  res = 0;
705  break;
706 
707  case UnstructuredMapping::quadrilateral :
708  switch(c) {
709  case(0):
710  if ( n==0 || n==3 )
711  res=-.5;
712  else
713  res= .5;
714  break;
715  case(1):
716  if ( n==0 || n==1 )
717  res=-.5;
718  else
719  res= .5;
720  break;
721  }
722  break;
723 
724  case UnstructuredMapping::hexahedron :
725  switch(c) {
726  case(0):
727  if ( n==0 || n==3 || n==4 || n==7 )
728  res = -.25;
729  else
730  res = .25;
731  break;
732  case(1):
733  if ( n==0 || n==1 || n==4 || n==5 )
734  res = -.25;
735  else
736  res = .25;
737  break;
738  case(2):
739  if ( n==0 || n==1 || n==2 || n==3 )
740  res = -.25;
741  else
742  res = .25;
743  break;
744  default:
745  break;
746  }
747  break;
748  case UnstructuredMapping::pyramid :
749  switch(c) {
750  case(0):
751  if ( n==0 || n==3 || n==4 )
752  res = -1./3;
753  else
754  res = 1./3.;
755  break;
756  case(1):
757  if ( n==0 || n==1 || n==4 )
758  res = -1./3;
759  else
760  res = 1./3.;
761  break;
762  case(2):
763  if ( n==4 )
764  res = 2.;
765  else
766  res = -.5;
767  break;
768  }
769  break;
770 
771  default:
772  res = -1;
773  break;
774  }
775 
776  return res;
777 }
778 
779 #endif
Generated on Fri Jan 4 2013 10:17:56 for Overture by   doxygen 1.8.3