PRAgMaTIc  master
ElementProperty.h
Go to the documentation of this file.
1 /* Copyright (C) 2010 Imperial College London and others.
2  *
3  * Please see the AUTHORS file in the main source directory for a
4  * full list of copyright holders.
5  *
6  * Gerard Gorman
7  * Applied Modelling and Computation Group
8  * Department of Earth Science and Engineering
9  * Imperial College London
10  *
11  * g.gorman@imperial.ac.uk
12  *
13  * Redistribution and use in source and binary forms, with or without
14  * modification, are permitted provided that the following conditions
15  * are met:
16  * 1. Redistributions of source code must retain the above copyright
17  * notice, this list of conditions and the following disclaimer.
18  * 2. Redistributions in binary form must reproduce the above
19  * copyright notice, this list of conditions and the following
20  * disclaimer in the documentation and/or other materials provided
21  * with the distribution.
22  *
23  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
24  * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
25  * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
26  * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
27  * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS
28  * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
29  * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
30  * TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
31  * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
32  * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
33  * TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF
34  * THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
35  * SUCH DAMAGE.
36  */
37 
38 #ifndef ELEMENTPROPERTY_H
39 #define ELEMENTPROPERTY_H
40 
41 #include <algorithm>
42 #include <cassert>
43 #include <cstdlib>
44 #include <iostream>
45 #include <map>
46 #include <vector>
47 #include <set>
48 #include <cmath>
49 
50 #include <cfloat>
51 
64 template<typename real_t>
65 class ElementProperty{
66  public:
72  ElementProperty(const real_t *x0, const real_t *x1, const real_t *x2): inv2(0.5), inv3(1.0/3.0), inv4(0.25), inv6(1.0/6.0), lipnikov_const2d(20.784609690826528), lipnikov_const3d(1832.8207768355312), condition_const2d(20.784609690826528), condition_const3d(940.60406122874042){
73  orientation = 1;
74 
75  double A = area(x0, x1, x2);
76  if(A<0)
77  orientation = -1;
78  }
79 
86  ElementProperty(const real_t *x0, const real_t *x1, const real_t *x2, const real_t *x3): inv2(0.5), inv3(1.0/3.0), inv4(0.25), inv6(1.0/6.0), lipnikov_const2d(20.784609690826528), lipnikov_const3d(1832.8207768355312), condition_const2d(20.784609690826528), condition_const3d(940.60406122874042){
87  orientation = 1;
88 
89  double V = volume(x0, x1, x2, x3);
90  if(V<0)
91  orientation = -1;
92  }
93 
99  real_t area(const real_t *x0, const real_t *x1, const real_t *x2) const{
100  real_t x01 = (x0[0] - x1[0]);
101  real_t y01 = (x0[1] - x1[1]);
102 
103  real_t x02 = (x0[0] - x2[0]);
104  real_t y02 = (x0[1] - x2[1]);
105 
106  return orientation*inv2*(y02*x01 - y01*x02);
107  }
108 
115  real_t volume(const real_t *x0, const real_t *x1, const real_t *x2, const real_t *x3) const{
116 
117  real_t x01 = (x0[0] - x1[0]);
118  real_t x02 = (x0[0] - x2[0]);
119  real_t x03 = (x0[0] - x3[0]);
120 
121  real_t y01 = (x0[1] - x1[1]);
122  real_t y02 = (x0[1] - x2[1]);
123  real_t y03 = (x0[1] - x3[1]);
124 
125  real_t z01 = (x0[2] - x1[2]);
126  real_t z02 = (x0[2] - x2[2]);
127  real_t z03 = (x0[2] - x3[2]);
128 
129  return orientation*inv6*(-x03*(z02*y01 - z01*y02) + x02*(z03*y01 - z01*y03) - x01*(z03*y02 - z02*y03));
130  }
131 
138  template<int dim>
139  double length(const real_t x0[], const real_t x1[], const double m[]) const{
140  if(dim==2){
141  return length2d(x0, x1, m);
142  }else{ //if(dim==3)
143  return length3d(x0, x1, m);
144  }
145  }
146 
153  static double length2d(const real_t x0[], const real_t x1[], const double m[]){
154  double x=x0[0] - x1[0];
155  double y=x0[1] - x1[1];
156 
157  assert((m[1]*x + m[2]*y)*y + (m[0]*x + m[1]*y)*x >= 0.0);
158 
159  return sqrt(((m[1]*x + m[2]*y)*y + (m[0]*x + m[1]*y)*x));
160  }
161 
168  static double length3d(const real_t x0[], const real_t x1[], const double m[]){
169  double x=x0[0] - x1[0];
170  double y=x0[1] - x1[1];
171  double z=x0[2] - x1[2];
172 
173  assert(z*(z*m[5] + y*m[4] + x*m[2]) +
174  y*(z*m[4] + y*m[3] + x*m[1]) +
175  x*(z*m[2] + y*m[1] + x*m[0]) >= 0.0);
176 
177  return sqrt(z*(z*m[5] + y*m[4] + x*m[2]) +
178  y*(z*m[4] + y*m[3] + x*m[1]) +
179  x*(z*m[2] + y*m[1] + x*m[0]));
180  }
181 
195  double lipnikov(const real_t *x0, const real_t *x1, const real_t *x2,
196  const double *m0, const double *m1, const double *m2){
197  // Metric tensor averaged over the element
198  double m00 = (m0[0] + m1[0] + m2[0])*inv3;
199  double m01 = (m0[1] + m1[1] + m2[1])*inv3;
200  double m11 = (m0[2] + m1[2] + m2[2])*inv3;
201 
202  return lipnikov(x0, x1, x2, m00, m01, m11);
203  }
204 
218  double lipnikov(const real_t *x0, const real_t *x1, const real_t *x2,
219  double m00, double m01, double m11){
220  // l is the length of the perimeter, measured in metric space
221  double x01 = x0[0] - x1[0];
222  double y01 = x0[1] - x1[1];
223  double x02 = x0[0] - x2[0];
224  double y02 = x0[1] - x2[1];
225  double x21 = x2[0] - x1[0];
226  double y21 = x2[1] - x1[1];
227 
228  double l =
229  sqrt(y01*(y01*m11 + x01*m01) +
230  x01*(y01*m01 + x01*m00))+
231  sqrt(y02*(y02*m11 + x02*m01) +
232  x02*(y02*m01 + x02*m00))+
233  sqrt(y21*(y21*m11 + x21*m01) +
234  x21*(y21*m01 + x21*m00));
235 
236  double invl = 1.0/l;
237 
238  // Area in physical space
239  double a=orientation*inv2*(y02*x01 - y01*x02);
240 
241  // Area in metric space
242  double a_m = a*sqrt(m00*m11 - m01*m01);
243 
244  // Function
245  double f = std::min(l*inv3, 3.0*invl);
246  double tf = f * (2.0 - f);
247  double F = tf*tf*tf;
248  double quality = lipnikov_const2d*a_m*F*invl*invl;
249 
250  return quality;
251  }
252 
253  // Gradient of lipnikov functional n0 using a central difference approximation.
254  void lipnikov_grad(int moving,
255  const double *x0, const double *x1, const double *x2,
256  const double *m0,
257  double *grad){
258  const double sqrt_eps = sqrt(DBL_EPSILON);
259 
260  // df/dx, df/dy
261  for(size_t i=0;i<2;i++){
262  double h = std::max(fabs(sqrt_eps*x0[i]), sqrt_eps);
263 
264  volatile double xnh = x0[i]-h;
265  volatile double xph = x0[i]+h;
266 
267  double Xn[] = {x0[0], x0[1]};
268  Xn[i] = xnh;
269  double Fxnh = lipnikov(Xn, x1, x2, m0[0], m0[1], m0[2]);
270 
271  double Xp[] = {x0[0], x0[1]};
272  Xp[i] = xph;
273  double Fxph = lipnikov(Xp, x1, x2, m0[0], m0[1], m0[2]);
274 
275  double two_dx = xph - xnh;
276  grad[i] = (Fxph - Fxnh)/two_dx;
277  }
278 
279  return;
280  }
281 
297  double lipnikov(const real_t *x0, const real_t *x1, const real_t *x2, const real_t *x3,
298  const double *m0, const double *m1, const double *m2, const double *m3){
299  // Metric tensor
300  double m00 = (m0[0] + m1[0] + m2[0] + m3[0])*inv4;
301  double m01 = (m0[1] + m1[1] + m2[1] + m3[1])*inv4;
302  double m02 = (m0[2] + m1[2] + m2[2] + m3[2])*inv4;
303  double m11 = (m0[3] + m1[3] + m2[3] + m3[3])*inv4;
304  double m12 = (m0[4] + m1[4] + m2[4] + m3[4])*inv4;
305  double m22 = (m0[5] + m1[5] + m2[5] + m3[5])*inv4;
306 
307  // l is the length of the edges of the tet, in metric space
308  double z01 = (x0[2] - x1[2]);
309  double y01 = (x0[1] - x1[1]);
310  double x01 = (x0[0] - x1[0]);
311 
312  double z12 = (x1[2] - x2[2]);
313  double y12 = (x1[1] - x2[1]);
314  double x12 = (x1[0] - x2[0]);
315 
316  double z02 = (x0[2] - x2[2]);
317  double y02 = (x0[1] - x2[1]);
318  double x02 = (x0[0] - x2[0]);
319 
320  double z03 = (x0[2] - x3[2]);
321  double y03 = (x0[1] - x3[1]);
322  double x03 = (x0[0] - x3[0]);
323 
324  double z13 = (x1[2] - x3[2]);
325  double y13 = (x1[1] - x3[1]);
326  double x13 = (x1[0] - x3[0]);
327 
328  double z23 = (x2[2] - x3[2]);
329  double y23 = (x2[1] - x3[1]);
330  double x23 = (x2[0] - x3[0]);
331 
332  double dl0 = (z01*(z01*m22 + y01*m12 + x01*m02) + y01*(z01*m12 + y01*m11 + x01*m01) + x01*(z01*m02 + y01*m01 + x01*m00));
333  double dl1 = (z12*(z12*m22 + y12*m12 + x12*m02) + y12*(z12*m12 + y12*m11 + x12*m01) + x12*(z12*m02 + y12*m01 + x12*m00));
334  double dl2 = (z02*(z02*m22 + y02*m12 + x02*m02) + y02*(z02*m12 + y02*m11 + x02*m01) + x02*(z02*m02 + y02*m01 + x02*m00));
335  double dl3 = (z03*(z03*m22 + y03*m12 + x03*m02) + y03*(z03*m12 + y03*m11 + x03*m01) + x03*(z03*m02 + y03*m01 + x03*m00));
336  double dl4 = (z13*(z13*m22 + y13*m12 + x13*m02) + y13*(z13*m12 + y13*m11 + x13*m01) + x13*(z13*m02 + y13*m01 + x13*m00));
337  double dl5 = (z23*(z23*m22 + y23*m12 + x23*m02) + y23*(z23*m12 + y23*m11 + x23*m01) + x23*(z23*m02 + y23*m01 + x23*m00));
338 
339  double l = sqrt(dl0)+sqrt(dl1)+sqrt(dl2)+sqrt(dl3)+sqrt(dl4)+sqrt(dl5);
340  double invl = 1.0/l;
341 
342  // Volume
343  double v=orientation*inv6*(-x03*(z02*y01 - z01*y02) + x02*(z03*y01 - z01*y03) - x01*(z03*y02 - z02*y03));
344 
345  // Volume in metric space
346  double v_m = v*sqrt(((m11*m22 - m12*m12)*m00 - (m01*m22 - m02*m12)*m01 + (m01*m12 - m02*m11)*m02));
347 
348  // Function
349  double f = std::min(l*inv6, 6*invl);
350  double tf = f * (2.0 - f);
351  double F = tf*tf*tf;
352  double quality = lipnikov_const3d * v_m * F *invl*invl*invl;
353 
354  return quality;
355  }
356 
369  double lipnikov(const real_t *x0, const real_t *x1, const real_t *x2, const real_t *x3,
370  const double *m0){
371  // Metric tensor
372  double m00 = m0[0];
373  double m01 = m0[1];
374  double m02 = m0[2];
375  double m11 = m0[3];
376  double m12 = m0[4];
377  double m22 = m0[5];
378 
379  // l is the length of the edges of the tet, in metric space
380  double z01 = (x0[2] - x1[2]);
381  double y01 = (x0[1] - x1[1]);
382  double x01 = (x0[0] - x1[0]);
383 
384  double z12 = (x1[2] - x2[2]);
385  double y12 = (x1[1] - x2[1]);
386  double x12 = (x1[0] - x2[0]);
387 
388  double z02 = (x0[2] - x2[2]);
389  double y02 = (x0[1] - x2[1]);
390  double x02 = (x0[0] - x2[0]);
391 
392  double z03 = (x0[2] - x3[2]);
393  double y03 = (x0[1] - x3[1]);
394  double x03 = (x0[0] - x3[0]);
395 
396  double z13 = (x1[2] - x3[2]);
397  double y13 = (x1[1] - x3[1]);
398  double x13 = (x1[0] - x3[0]);
399 
400  double z23 = (x2[2] - x3[2]);
401  double y23 = (x2[1] - x3[1]);
402  double x23 = (x2[0] - x3[0]);
403 
404  double dl0 = (z01*(z01*m22 + y01*m12 + x01*m02) + y01*(z01*m12 + y01*m11 + x01*m01) + x01*(z01*m02 + y01*m01 + x01*m00));
405  double dl1 = (z12*(z12*m22 + y12*m12 + x12*m02) + y12*(z12*m12 + y12*m11 + x12*m01) + x12*(z12*m02 + y12*m01 + x12*m00));
406  double dl2 = (z02*(z02*m22 + y02*m12 + x02*m02) + y02*(z02*m12 + y02*m11 + x02*m01) + x02*(z02*m02 + y02*m01 + x02*m00));
407  double dl3 = (z03*(z03*m22 + y03*m12 + x03*m02) + y03*(z03*m12 + y03*m11 + x03*m01) + x03*(z03*m02 + y03*m01 + x03*m00));
408  double dl4 = (z13*(z13*m22 + y13*m12 + x13*m02) + y13*(z13*m12 + y13*m11 + x13*m01) + x13*(z13*m02 + y13*m01 + x13*m00));
409  double dl5 = (z23*(z23*m22 + y23*m12 + x23*m02) + y23*(z23*m12 + y23*m11 + x23*m01) + x23*(z23*m02 + y23*m01 + x23*m00));
410 
411  double l = sqrt(dl0)+sqrt(dl1)+sqrt(dl2)+sqrt(dl3)+sqrt(dl4)+sqrt(dl5);
412  double invl = 1.0/l;
413 
414  // Volume
415  double v=orientation*inv6*(-x03*(z02*y01 - z01*y02) + x02*(z03*y01 - z01*y03) - x01*(z03*y02 - z02*y03));
416 
417  // Volume in metric space
418  double v_m = v*sqrt(((m11*m22 - m12*m12)*m00 - (m01*m22 - m02*m12)*m01 + (m01*m12 - m02*m11)*m02));
419 
420  // Function
421  double f = std::min(l*inv6, 6*invl);
422  double tf = f * (2.0 - f);
423  double F = tf*tf*tf;
424  double quality = lipnikov_const3d * v_m * F *invl*invl*invl;
425 
426  return quality;
427  }
428 
429 
430  // Gradient of lipnikov functional n0 using a central difference approximation.
431  void lipnikov_grad(int moving,
432  const double *x0, const double *x1, const double *x2, const double *x3,
433  const double *m0,
434  double *grad){
435  const double sqrt_eps = sqrt(DBL_EPSILON);
436 
437  // df/dx, df/dy, df/dz
438  for(size_t i=0;i<3;i++){
439  double h = std::max(fabs(sqrt_eps*x0[i]), sqrt_eps);
440 
441  volatile double xnh = x0[i]-h;
442  volatile double xph = x0[i]+h;
443 
444  double Xn[] = {x0[0], x0[1], x0[2]};
445  Xn[i] = xnh;
446  double Fxnh = lipnikov(Xn, x1, x2, x3, m0);
447 
448  double Xp[] = {x0[0], x0[1], x0[2]};
449  Xp[i] = xph;
450  double Fxph = lipnikov(Xp, x1, x2, x3, m0);
451 
452  double two_dx = xph - xnh;
453  grad[i] = (Fxph - Fxnh)/two_dx;
454  }
455 
456  return;
457  }
458 
472  real_t sliver(const real_t *x0, const real_t *x1, const real_t *x2, const real_t *x3,
473  const double *m0, const double *m1, const double *m2, const double *m3){
474  // Metric tensor
475  double m00 = (m0[0] + m1[0] + m2[0] + m3[0])*inv4;
476  double m01 = (m0[1] + m1[1] + m2[1] + m3[1])*inv4;
477  double m02 = (m0[2] + m1[2] + m2[2] + m3[2])*inv4;
478  double m11 = (m0[3] + m1[3] + m2[3] + m3[3])*inv4;
479  double m12 = (m0[4] + m1[4] + m2[4] + m3[4])*inv4;
480  double m22 = (m0[5] + m1[5] + m2[5] + m3[5])*inv4;
481 
482  double z01 = (x0[2] - x1[2]);
483  double y01 = (x0[1] - x1[1]);
484  double x01 = (x0[0] - x1[0]);
485 
486  double z12 = (x1[2] - x2[2]);
487  double y12 = (x1[1] - x2[1]);
488  double x12 = (x1[0] - x2[0]);
489 
490  double z02 = (x0[2] - x2[2]);
491  double y02 = (x0[1] - x2[1]);
492  double x02 = (x0[0] - x2[0]);
493 
494  double z03 = (x0[2] - x3[2]);
495  double y03 = (x0[1] - x3[1]);
496  double x03 = (x0[0] - x3[0]);
497 
498  double z13 = (x1[2] - x3[2]);
499  double y13 = (x1[1] - x3[1]);
500  double x13 = (x1[0] - x3[0]);
501 
502  double z23 = (x2[2] - x3[2]);
503  double y23 = (x2[1] - x3[1]);
504  double x23 = (x2[0] - x3[0]);
505 
506  // l is the length of the edges of the tet, in metric space
507  double dl0 = (z01*(z01*m22 + y01*m12 + x01*m02) + y01*(z01*m12 + y01*m11 + x01*m01) + x01*(z01*m02 + y01*m01 + x01*m00));
508  double dl1 = (z12*(z12*m22 + y12*m12 + x12*m02) + y12*(z12*m12 + y12*m11 + x12*m01) + x12*(z12*m02 + y12*m01 + x12*m00));
509  double dl2 = (z02*(z02*m22 + y02*m12 + x02*m02) + y02*(z02*m12 + y02*m11 + x02*m01) + x02*(z02*m02 + y02*m01 + x02*m00));
510  double dl3 = (z03*(z03*m22 + y03*m12 + x03*m02) + y03*(z03*m12 + y03*m11 + x03*m01) + x03*(z03*m02 + y03*m01 + x03*m00));
511  double dl4 = (z13*(z13*m22 + y13*m12 + x13*m02) + y13*(z13*m12 + y13*m11 + x13*m01) + x13*(z13*m02 + y13*m01 + x13*m00));
512  double dl5 = (z23*(z23*m22 + y23*m12 + x23*m02) + y23*(z23*m12 + y23*m11 + x23*m01) + x23*(z23*m02 + y23*m01 + x23*m00));
513 
514  // Volume
515  double v=orientation*inv6*(-x03*(z02*y01 - z01*y02) + x02*(z03*y01 - z01*y03) - (x0[0] - x1[0])*(z03*y02 - z02*y03));
516 
517  // Volume in metric space
518  double v_m = v*sqrt(((m11*m22 - m12*m12)*m00 - (m01*m22 - m02*m12)*m01 + (m01*m12 - m02*m11)*m02));
519 
520  // Sliver functional
521  double ts = dl0+dl1+dl2+dl3+dl4+dl5;
522  double sliver = 15552.0*(v_m*v_m)/(ts*ts*ts);
523 
524  return sliver;
525  }
526 
540  double condition(const real_t *x0, const real_t *x1, const real_t *x2,
541  const double *m0, const double *m1, const double *m2){
542  // Metric tensor averaged over the element
543  double m00 = (m0[0] + m1[0] + m2[0])*inv3;
544  double m01 = (m0[1] + m1[1] + m2[1])*inv3;
545  double m11 = (m0[2] + m1[2] + m2[2])*inv3;
546 
547  return condition(x0, x1, x2, m00, m01, m11);
548  }
549 
550 
564  double condition(const real_t *x0, const real_t *x1, const real_t *x2,
565  double m00, double m01, double m11){
566  // l is the length of the perimeter, measured in metric space
567  double x01 = x0[0] - x1[0];
568  double y01 = x0[1] - x1[1];
569  double x02 = x0[0] - x2[0];
570  double y02 = x0[1] - x2[1];
571  double x21 = x2[0] - x1[0];
572  double y21 = x2[1] - x1[1];
573 
574  double l = y01*(y01*m11 + x01*m01) +
575  x01*(y01*m01 + x01*m00) +
576  y02*(y02*m11 + x02*m01) +
577  x02*(y02*m01 + x02*m00) +
578  y21*(y21*m11 + x21*m01) +
579  x21*(y21*m01 + x21*m00);
580 
581  double invl = 1.0/l;
582 
583  // Area in physical space
584  double a=orientation*inv2*(y02*x01 - y01*x02);
585 
586  // Area in metric space
587  double a_m = a*sqrt(m00*m11 - m01*m01);
588 
589  // Function
590  double quality = condition_const2d*a_m*invl;
591 
592  return quality;
593  }
594 
610  double condition(const real_t *x0, const real_t *x1, const real_t *x2, const real_t *x3,
611  const double *m0, const double *m1, const double *m2, const double *m3){
612  // Metric tensor
613  double m00 = (m0[0] + m1[0] + m2[0] + m3[0])*inv4;
614  double m01 = (m0[1] + m1[1] + m2[1] + m3[1])*inv4;
615  double m02 = (m0[2] + m1[2] + m2[2] + m3[2])*inv4;
616  double m11 = (m0[3] + m1[3] + m2[3] + m3[3])*inv4;
617  double m12 = (m0[4] + m1[4] + m2[4] + m3[4])*inv4;
618  double m22 = (m0[5] + m1[5] + m2[5] + m3[5])*inv4;
619 
620  // l is the length of the edges of the tet, in metric space
621  double z01 = (x0[2] - x1[2]);
622  double y01 = (x0[1] - x1[1]);
623  double x01 = (x0[0] - x1[0]);
624 
625  double z12 = (x1[2] - x2[2]);
626  double y12 = (x1[1] - x2[1]);
627  double x12 = (x1[0] - x2[0]);
628 
629  double z02 = (x0[2] - x2[2]);
630  double y02 = (x0[1] - x2[1]);
631  double x02 = (x0[0] - x2[0]);
632 
633  double z03 = (x0[2] - x3[2]);
634  double y03 = (x0[1] - x3[1]);
635  double x03 = (x0[0] - x3[0]);
636 
637  double z13 = (x1[2] - x3[2]);
638  double y13 = (x1[1] - x3[1]);
639  double x13 = (x1[0] - x3[0]);
640 
641  double z23 = (x2[2] - x3[2]);
642  double y23 = (x2[1] - x3[1]);
643  double x23 = (x2[0] - x3[0]);
644 
645  double dl0 = (z01*(z01*m22 + y01*m12 + x01*m02) + y01*(z01*m12 + y01*m11 + x01*m01) + x01*(z01*m02 + y01*m01 + x01*m00));
646  double dl1 = (z12*(z12*m22 + y12*m12 + x12*m02) + y12*(z12*m12 + y12*m11 + x12*m01) + x12*(z12*m02 + y12*m01 + x12*m00));
647  double dl2 = (z02*(z02*m22 + y02*m12 + x02*m02) + y02*(z02*m12 + y02*m11 + x02*m01) + x02*(z02*m02 + y02*m01 + x02*m00));
648  double dl3 = (z03*(z03*m22 + y03*m12 + x03*m02) + y03*(z03*m12 + y03*m11 + x03*m01) + x03*(z03*m02 + y03*m01 + x03*m00));
649  double dl4 = (z13*(z13*m22 + y13*m12 + x13*m02) + y13*(z13*m12 + y13*m11 + x13*m01) + x13*(z13*m02 + y13*m01 + x13*m00));
650  double dl5 = (z23*(z23*m22 + y23*m12 + x23*m02) + y23*(z23*m12 + y23*m11 + x23*m01) + x23*(z23*m02 + y23*m01 + x23*m00));
651 
652  double l2 = dl0+dl1+dl2+dl3+dl4+dl5;
653  double invl2 = 1.0/l2;
654 
655  // areas of the faces of the tet, in metric space (Sympy)
656  double tmp00 = (2.0*m00*m12 - 2.0*m01*m02);
657  double tmp11 = (-2.0*m01*m12 + 2.0*m02*m11);
658  double tmp22 = (-2.0*m01*m22 + 2.0*m02*m12);
659  double tmp01 = (m00*m11 - m01*m01);
660  double tmp02 = (m00*m22 - m02*m02);
661  double tmp12 = (m11*m22 - m12*m12);
662 
663  double sq_area012 = tmp00*x01*x01*y02*z02 - tmp00*x01*x02*y01*z02 - tmp00*x01*x02*y02*z01 + tmp00*x02*x02*y01*z01 + tmp01*x01*x01*y02*y02 - 2.0*tmp01*x01*x02*y01*y02 + tmp01*x02*x02*y01*y01 + tmp02*x01*x01*z02*z02 - 2.0*tmp02*x01*x02*z01*z02 + tmp02*x02*x02*z01*z01 - tmp11*x01*y01*y02*z02 + tmp11*x01*y02*y02*z01 + tmp11*x02*y01*y01*z02 - tmp11*x02*y01*y02*z01 + tmp12*y01*y01*z02*z02 - 2.0*tmp12*y01*y02*z01*z02 + tmp12*y02*y02*z01*z01 - tmp22*x01*y01*z02*z02 + tmp22*x01*y02*z01*z02 + tmp22*x02*y01*z01*z02 - tmp22*x02*y02*z01*z01;
664 
665  double sq_area013 = tmp00*x01*x01*y13*z13 - tmp00*x01*x13*y01*z13 - tmp00*x01*x13*y13*z01 + tmp00*x13*x13*y01*z01 + tmp01*x01*x01*y13*y13 - 2.0*tmp01*x01*x13*y01*y13 + tmp01*x13*x13*y01*y01 + tmp02*x01*x01*z13*z13 - 2.0*tmp02*x01*x13*z01*z13 + tmp02*x13*x13*z01*z01 - tmp11*x01*y01*y13*z13 + tmp11*x01*y13*y13*z01 + tmp11*x13*y01*y01*z13 - tmp11*x13*y01*y13*z01 + tmp12*y01*y01*z13*z13 - 2.0*tmp12*y01*y13*z01*z13 + tmp12*y13*y13*z01*z01 - tmp22*x01*y01*z13*z13 + tmp22*x01*y13*z01*z13 + tmp22*x13*y01*z01*z13 - tmp22*x13*y13*z01*z01;
666 
667  double sq_area123 = tmp00*x12*x12*y23*z23 - tmp00*x12*x23*y12*z23 - tmp00*x12*x23*y23*z12 + tmp00*x23*x23*y12*z12 + tmp01*x12*x12*y23*y23 - 2.0*tmp01*x12*x23*y12*y23 + tmp01*x23*x23*y12*y12 + tmp02*x12*x12*z23*z23 - 2.0*tmp02*x12*x23*z12*z23 + tmp02*x23*x23*z12*z12 - tmp11*x12*y12*y23*z23 + tmp11*x12*y23*y23*z12 + tmp11*x23*y12*y12*z23 - tmp11*x23*y12*y23*z12 + tmp12*y12*y12*z23*z23 - 2.0*tmp12*y12*y23*z12*z23 + tmp12*y23*y23*z12*z12 - tmp22*x12*y12*z23*z23 + tmp22*x12*y23*z12*z23 + tmp22*x23*y12*z12*z23 - tmp22*x23*y23*z12*z12;
668 
669  double sq_area023 = tmp00*x01*x01*y13*z13 - tmp00*x01*x13*y01*z13 - tmp00*x01*x13*y13*z01 + tmp00*x13*x13*y01*z01 + tmp01*x01*x01*y13*y13 - 2.0*tmp01*x01*x13*y01*y13 + tmp01*x13*x13*y01*y01 + tmp02*x01*x01*z13*z13 - 2.0*tmp02*x01*x13*z01*z13 + tmp02*x13*x13*z01*z01 - tmp11*x01*y01*y13*z13 + tmp11*x01*y13*y13*z01 + tmp11*x13*y01*y01*z13 - tmp11*x13*y01*y13*z01 + tmp12*y01*y01*z13*z13 - 2.0*tmp12*y01*y13*z01*z13 + tmp12*y13*y13*z01*z01 - tmp22*x01*y01*z13*z13 + tmp22*x01*y13*z01*z13 + tmp22*x13*y01*z01*z13 - tmp22*x13*y13*z01*z01;
670 
671  double area2 = sq_area012 + sq_area013 + sq_area123 + sq_area023;
672  double invarea2 = 1.0/area2;
673 
674  // Volume
675  double v=orientation*inv6*(-x03*(z02*y01 - z01*y02) + x02*(z03*y01 - z01*y03) - x01*(z03*y02 - z02*y03));
676 
677  // Volume in metric space
678  double v_m = v*sqrt(((m11*m22 - m12*m12)*m00 - (m01*m22 - m02*m12)*m01 + (m01*m12 - m02*m11)*m02));
679 
680  // Function
681  double quality = condition_const3d * v_m * sqrt(invarea2*invl2);
682 
683  return quality;
684  }
685 
686  int getOrientation() {
687  return orientation;
688  }
689 
690  private:
691  const double inv2;
692  const double inv3;
693  const double inv4;
694  const double inv6;
695 
696  const double lipnikov_const2d; // 12.0*sqrt(3.0);
697  const double lipnikov_const3d; // pow(6.0, 4)*sqrt(2.0);
698  const double condition_const2d; // (4/sqrt(3))*3^2 = 12.0*sqrt(3.0)
699  const double condition_const3d; // (12/sqrt(2))*6*(4*(4/sqrt(3))) = 1152*sqrt(6)
700 
701  int orientation;
702 };
703 #endif
generate_Steiner_ellipse_3d.x2
tuple x2
Definition: generate_Steiner_ellipse_3d.py:9
ElementProperty::ElementProperty
ElementProperty(const real_t *x0, const real_t *x1, const real_t *x2)
Definition: ElementProperty.h:72
ElementProperty::lipnikov_grad
void lipnikov_grad(int moving, const double *x0, const double *x1, const double *x2, const double *m0, double *grad)
Definition: ElementProperty.h:254
ElementProperty::length2d
static double length2d(const real_t x0[], const real_t x1[], const double m[])
Definition: ElementProperty.h:153
fit_ellipsoid_3d.A
tuple A
Definition: fit_ellipsoid_3d.py:16
ElementProperty::length3d
static double length3d(const real_t x0[], const real_t x1[], const double m[])
Definition: ElementProperty.h:168
ElementProperty
Calculates a number of element properties.
Definition: ElementProperty.h:65
ElementProperty::condition
double condition(const real_t *x0, const real_t *x1, const real_t *x2, double m00, double m01, double m11)
Definition: ElementProperty.h:564
ElementProperty::lipnikov_grad
void lipnikov_grad(int moving, const double *x0, const double *x1, const double *x2, const double *x3, const double *m0, double *grad)
Definition: ElementProperty.h:431
ElementProperty::condition
double condition(const real_t *x0, const real_t *x1, const real_t *x2, const real_t *x3, const double *m0, const double *m1, const double *m2, const double *m3)
Definition: ElementProperty.h:610
ElementProperty::condition
double condition(const real_t *x0, const real_t *x1, const real_t *x2, const double *m0, const double *m1, const double *m2)
Definition: ElementProperty.h:540
ElementProperty::lipnikov
double lipnikov(const real_t *x0, const real_t *x1, const real_t *x2, const real_t *x3, const double *m0, const double *m1, const double *m2, const double *m3)
Definition: ElementProperty.h:297
ElementProperty::volume
real_t volume(const real_t *x0, const real_t *x1, const real_t *x2, const real_t *x3) const
Definition: ElementProperty.h:115
ElementProperty::area
real_t area(const real_t *x0, const real_t *x1, const real_t *x2) const
Definition: ElementProperty.h:99
ElementProperty::lipnikov
double lipnikov(const real_t *x0, const real_t *x1, const real_t *x2, const real_t *x3, const double *m0)
Definition: ElementProperty.h:369
ElementProperty::lipnikov
double lipnikov(const real_t *x0, const real_t *x1, const real_t *x2, const double *m0, const double *m1, const double *m2)
Definition: ElementProperty.h:195
ElementProperty::lipnikov
double lipnikov(const real_t *x0, const real_t *x1, const real_t *x2, double m00, double m01, double m11)
Definition: ElementProperty.h:218
ElementProperty::length
double length(const real_t x0[], const real_t x1[], const double m[]) const
Definition: ElementProperty.h:139
generate_Steiner_ellipse_3d.x1
tuple x1
Definition: generate_Steiner_ellipse_3d.py:8
ellipse_convergence.a
tuple a
Definition: ellipse_convergence.py:114
ElementProperty::sliver
real_t sliver(const real_t *x0, const real_t *x1, const real_t *x2, const real_t *x3, const double *m0, const double *m1, const double *m2, const double *m3)
Definition: ElementProperty.h:472
generate_Steiner_ellipse_3d.x3
tuple x3
Definition: generate_Steiner_ellipse_3d.py:10
ellipse_convergence.V
tuple V
Definition: ellipse_convergence.py:111
ElementProperty::ElementProperty
ElementProperty(const real_t *x0, const real_t *x1, const real_t *x2, const real_t *x3)
Definition: ElementProperty.h:86
Generated on Mon Feb 23 2015 22:50:44 for PRAgMaTIc by   doxygen 1.8.7