Added "normalized corner Jacobian" quality measure for triangles and quads...

Added "normalized corner Jacobian" quality measure for triangles and quads (not yet used nor tested)

Mesh/qualityMeasures.cpp

@@ -333,3 +333,313 @@ double qmQuadrangleAngles (MQuadrangle *e) {
 
 }
 
+
+namespace {
+
+
+// Compute unit vector
+inline void unitVec(const double &x0, const double &y0, const double &z0,
+                    const double &x1, const double &y1, const double &z1,
+                    double &x01n, double &y01n, double &z01n)
+{
+  const double x01 = x1-x0, y01 = y1-y0, z01 = z1-z0;
+  const double invL01 = 1./sqrt(x01*x01 + y01*y01 + z01*z01);
+  x01n = x01*invL01; y01n = y01*invL01; z01n = z01*invL01;
+}
+
+
+// Compute unit vector and gradients w.r.t. x0, y0, z0, x1, y1, z1
+inline void unitVecAndGrad(const double &x0, const double &y0, const double &z0,
+                           const double &x1, const double &y1, const double &z1,
+                           fullMatrix<double> &res)
+{
+  const double x01 = x1-x0, y01 = y1-y0, z01 = z1-z0;
+  const double x01Sq = x01*x01, y01Sq = y01*y01, z01Sq = z01*z01;
+  const double invL01Sq = 1./(x01Sq+y01Sq+z01Sq), invL01 = sqrt(invL01Sq);
+  const double invL01Cu = invL01Sq*invL01;
+
+  double &dx01ndx0 = res(0, 0), &dx01ndy0 = res(0, 1), &dx01ndz0 = res(0, 2),
+          &dx01ndx1 = res(0, 3), &dx01ndy1 = res(0, 4), &dx01ndz1 = res(0, 5), &x01n = res(0, 6);
+  double &dy01ndx0 = res(1, 0), &dy01ndy0 = res(1, 1), &dy01ndz0 = res(1, 2),
+          &dy01ndx1 = res(1, 3), &dy01ndy1 = res(1, 4), &dy01ndz1 = res(1, 5), &y01n = res(1, 6);
+  double &dz01ndx0 = res(2, 0), &dz01ndy0 = res(2, 1), &dz01ndz0 = res(2, 2),
+          &dz01ndx1 = res(2, 3), &dz01ndy1 = res(2, 4), &dz01ndz1 = res(2, 5), &z01n = res(2, 6);
+
+  dx01ndx0 = x01Sq*invL01Cu-invL01;
+  dx01ndy0 = x01*y01*invL01Cu;
+  dx01ndz0 = x01*z01*invL01Cu;
+  dx01ndx1 = -dx01ndx0;
+  dx01ndy1 = -dx01ndy0;
+  dx01ndz1 = -dx01ndz0;
+  x01n = x01*invL01;
+
+  dy01ndx0 = dx01ndy0;
+  dy01ndy0 = y01Sq*invL01Cu-invL01;
+  dy01ndz0 = y01*z01*invL01Cu;
+  dy01ndx1 = -dy01ndx0;
+  dy01ndy1 = -dy01ndy0;
+  dy01ndz1 = -dy01ndz0;
+  y01n = y01*invL01;
+
+  dz01ndx0 = dx01ndz0;
+  dz01ndy0 = dy01ndz0;
+  dz01ndz0 = z01Sq*invL01Cu-invL01;
+  dz01ndx1 = -dx01ndz0;
+  dz01ndy1 = -dz01ndy0;
+  dz01ndz1 = -dz01ndz0;
+  z01n = z01*invL01;
+}
+
+
+// Given vectors a and b, compute area of triangle defined by both vectors (times a factor)
+inline double triArea(const double fact,
+                      const double &xa, const double &ya, const double &za,
+                      const double &xb, const double &yb, const double &zb)
+{
+  const double zn = xa*yb - xb*ya, yn = za*xb - xa*zb, xn = ya*zb - za*yb;
+  return fact * sqrt(xn*xn + yn*yn + zn*zn);
+}
+
+
+// Given vectors (0, 1) and (0, 2) and their gradients, compute area of triangle
+// defined by both vectors and grad. w.r.t. x0, y0, z0, x1, y1, z1, x2, y2, z2 (times a factor)
+inline void triAreaAndGrad(const double fact, const fullMatrix<double> &vga,
+                           const fullMatrix<double> &vgb, fullVector<double> &res)
+{
+  const double &dxadx0 = vga(0, 0), &dxady0 = vga(0, 1), &dxadz0 = vga(0, 2),
+                &dxadx1 = vga(0, 3), &dxady1 = vga(0, 4), &dxadz1 = vga(0, 5), &xa = vga(0, 6);
+  const double &dyadx0 = vga(1, 0), &dyady0 = vga(1, 1), &dyadz0 = vga(1, 2),
+                &dyadx1 = vga(1, 3), &dyady1 = vga(1, 4), &dyadz1 = vga(1, 5), &ya = vga(1, 6);
+  const double &dzadx0 = vga(2, 0), &dzady0 = vga(2, 1), &dzadz0 = vga(2, 2),
+                &dzadx1 = vga(2, 3), &dzady1 = vga(2, 4), &dzadz1 = vga(2, 5), &za = vga(2, 6);
+
+  const double &dxbdx0 = vgb(0, 0), &dxbdy0 = vgb(0, 1), &dxbdz0 = vgb(0, 2),
+                &dxbdx2 = vgb(0, 3), &dxbdy2 = vgb(0, 4), &dxbdz2 = vgb(0, 5), &xb = vgb(0, 6);
+  const double &dybdx0 = vgb(1, 0), &dybdy0 = vgb(1, 1), &dybdz0 = vgb(1, 2),
+                &dybdx2 = vgb(1, 3), &dybdy2 = vgb(1, 4), &dybdz2 = vgb(1, 5), &yb = vgb(1, 6);
+  const double &dzbdx0 = vgb(2, 0), &dzbdy0 = vgb(2, 1), &dzbdz0 = vgb(2, 2),
+                &dzbdx2 = vgb(2, 3), &dzbdy2 = vgb(2, 4), &dzbdz2 = vgb(2, 5), &zb = vgb(2, 6);
+
+  const double xn = ya*zb - za*yb;
+  const double dxndx0 = ya*dzbdx0 + dyadx0*zb - yb*dzadx0 - dybdx0*za;
+  const double dxndy0 = ya*dzbdy0 + dyady0*zb - yb*dzady0 - dybdy0*za;
+  const double dxndz0 = ya*dzbdz0 + dyadz0*zb - yb*dzadz0 - dybdz0*za;
+  const double dxndx1 = dyadx1*zb - yb*dzadx1;
+  const double dxndy1 = dyady1*zb - yb*dzady1;
+  const double dxndz1 = dyadz1*zb - yb*dzadz1;
+  const double dxndx2 = ya*dzbdx2 - dybdx2*za;
+  const double dxndy2 = ya*dzbdy2 - dybdy2*za;
+  const double dxndz2 = ya*dzbdz2 - dybdz2*za;
+
+  const double yn = za*xb - xa*zb;
+  const double dyndx0 = -xa*dzbdx0 - dxadx0*zb + xb*dzadx0 + dxbdx0*za;
+  const double dyndy0 = -xa*dzbdy0 - dxady0*zb + xb*dzady0 + dxbdy0*za;
+  const double dyndz0 = -xa*dzbdz0 - dxadz0*zb + xb*dzadz0 + dxbdz0*za;
+  const double dyndx1 = xb*dzadx1 - dxadx1*zb;
+  const double dyndy1 = xb*dzady1 - dxady1*zb;
+  const double dyndz1 = xb*dzadz1 - dxadz1*zb;
+  const double dyndx2 = dxbdx2*za - xa*dzbdx2;
+  const double dyndy2 = dxbdy2*za - xa*dzbdy2;
+  const double dyndz2 = dxbdz2*za - xa*dzbdz2;
+
+  const double zn = xa*yb - xb*ya;
+  const double dzndx0 = xa*dybdx0 + dxadx0*yb - xb*dyadx0 - dxbdx0*ya;
+  const double dzndy0 = xa*dybdy0 + dxady0*yb - xb*dyady0 - dxbdy0*ya;
+  const double dzndz0 = xa*dybdz0 + dxadz0*yb - xb*dyadz0 - dxbdz0*ya;
+  const double dzndx1 = dxadx1*yb - xb*dyadx1;
+  const double dzndy1 = dxady1*yb - xb*dyady1;
+  const double dzndz1 = dxadz1*yb - xb*dyadz1;
+  const double dzndx2 = xa*dybdx2 - dxbdx2*ya;
+  const double dzndy2 = xa*dybdy2 - dxbdy2*ya;
+  const double dzndz2 = xa*dybdz2 - dxbdz2*ya;
+
+  double &dNCJdx0 = res(0), &dNCJdy0 = res(1), &dNCJdz0 = res(2),
+          &dNCJdx1 = res(3), &dNCJdy1 = res(4), &dNCJdz1 = res(5),
+          &dNCJdx2 = res(6), &dNCJdy2 = res(7), &dNCJdz2 = res(8), &NCJ = res(9);
+
+  const double n = sqrt(xn*xn + yn*yn + zn*zn), invN = 1./n, fact2 = fact*invN;
+  dNCJdx0 = (dxndx0*xn + dyndx0*yn + dzndx0*zn) * fact2;
+  dNCJdy0 = (dxndy0*xn + dyndy0*yn + dzndy0*zn) * fact2;
+  dNCJdz0 = (dxndz0*xn + dyndz0*yn + dzndz0*zn) * fact2;
+  dNCJdx1 = (dxndx1*xn + dyndx1*yn + dzndx1*zn) * fact2;
+  dNCJdy1 = (dxndy1*xn + dyndy1*yn + dzndy1*zn) * fact2;
+  dNCJdz1 = (dxndz1*xn + dyndz1*yn + dzndz1*zn) * fact2;
+  dNCJdx2 = (dxndx2*xn + dyndx2*yn + dzndx2*zn) * fact2;
+  dNCJdy2 = (dxndy2*xn + dyndy2*yn + dzndy2*zn) * fact2;
+  dNCJdz2 = (dxndz2*xn + dyndz2*yn + dzndz2*zn) * fact2;
+  NCJ = fact * n;
+}
+
+
+inline void revertVG(const fullMatrix<double> &vg, fullMatrix<double> &res)
+{
+  res(0, 0) = -vg(0, 3); res(0, 1) = -vg(0, 4); res(0, 2) = -vg(0, 5); res(0, 6) = -vg(0, 6);
+  res(1, 0) = -vg(1, 3); res(1, 1) = -vg(1, 4); res(1, 2) = -vg(1, 5); res(1, 6) = -vg(1, 6);
+  res(2, 0) = -vg(2, 3); res(2, 1) = -vg(2, 4); res(2, 2) = -vg(2, 5); res(2, 6) = -vg(2, 6);
+}
+
+
+}
+
+
+void measuresTriangle::getNCJ(const double &x0, const double &y0, const double &z0,
+                              const double &x1, const double &y1, const double &z1,
+                              const double &x2, const double &y2, const double &z2,
+                              fullVector<double> &ncj)
+{
+  // Compute unit vectors for each edge
+  double x01n, y01n, z01n, x12n, y12n, z12n, x20n, y20n, z20n;
+  unitVec(x0, y0, z0, x1, y1, z1, x01n, y01n, z01n);
+  unitVec(x1, y1, z1, x2, y2, z2, x12n, y12n, z12n);
+  unitVec(x2, y2, z2, x0, y0, z0, x20n, y20n, z20n);
+
+  // Compute NCJ at vertex from unit vectors a and b as 0.5*||a^b||/A_equilateral
+  // Factor = 2./sqrt(3.) = 0.5/A_equilateral
+  static const double fact = 1.1547005383792517;
+  ncj(0) = triArea(fact, x01n, y01n, z01n, -x20n, -y20n, -z20n);
+  ncj(1) = triArea(fact, x12n, y12n, z12n, -x01n, -y01n, -z01n);
+  ncj(2) = triArea(fact, x20n, y20n, z20n, -x12n, -y12n, -z12n);
+}
+
+
+void measuresTriangle::getNCJAndGradients(const double &x0, const double &y0, const double &z0,
+                                          const double &x1, const double &y1, const double &z1,
+                                          const double &x2, const double &y2, const double &z2,
+                                          fullMatrix<double> &res)
+{
+  // Factor = 2./sqrt(3.) = 0.5/A_equilateral
+  static const double fact = 1.1547005383792517;
+
+  // Compute unit vector and its gradients for each edge
+  fullMatrix<double> vg01(3,7);
+  unitVecAndGrad(x0, y0, z0, x1, y1, z1, vg01);
+  fullMatrix<double> vg12(3,7);
+  unitVecAndGrad(x1, y1, z1, x2, y2, z2, vg12);
+  fullMatrix<double> vg20(3,7);
+  unitVecAndGrad(x2, y2, z2, x0, y0, z0, vg20);
+  fullMatrix<double> vg10(3,7);
+  revertVG(vg01, vg10);
+  fullMatrix<double> vg21(3,7);
+  revertVG(vg12, vg21);
+  fullMatrix<double> vg02(3,7);
+  revertVG(vg20, vg02);
+
+  // Compute NCJ at vertex 0 as 0.5*||u01^u02||/A_triEqui
+  // and gradients w.r.t. x0, y0, z0, x1, y1, z1, x2, y2, z2
+  fullVector<double> res0(10);
+  triAreaAndGrad(fact, vg01, vg02, res0);
+  res(0, 0) = res0(0); res(0, 1) = res0(1); res(0, 2) = res0(2);
+  res(0, 3) = res0(3); res(0, 4) = res0(4); res(0, 5) = res0(5);
+  res(0, 6) = res0(6); res(0, 7) = res0(7); res(0, 8) = res0(8);
+  res(0, 9) = res0(9);
+
+  // Compute NCJ at vertex 1 as 0.5*||u12^u10||/A_triEqui
+  // and gradients w.r.t. x1, y1, z1, x2, y2, z2, x0, y0, z0
+  fullVector<double> res1(10);
+  triAreaAndGrad(fact, vg12, vg10, res1);
+  res(1, 0) = res1(6); res(1, 1) = res1(7); res(1, 2) = res1(8);
+  res(1, 3) = res1(0); res(1, 4) = res1(1); res(1, 5) = res1(2);
+  res(1, 6) = res1(3); res(1, 7) = res1(4); res(1, 8) = res1(5);
+  res(1, 9) = res1(9);
+
+  // Compute NCJ at vertex 2 as 0.5*||u20^u21||/A_triEqui
+  // Compute NCJ at vertex 2 and gradients w.r.t. x2, y2, z2, x0, y0, z0, x1, y1, z1
+  fullVector<double> res2(10);
+  triAreaAndGrad(fact, vg20, vg21, res2);
+  res(2, 0) = res2(3); res(2, 1) = res2(4); res(2, 2) = res2(5);
+  res(2, 3) = res2(6); res(2, 4) = res2(7); res(2, 5) = res2(8);
+  res(2, 6) = res2(0); res(2, 7) = res2(1); res(2, 8) = res2(2);
+  res(2, 9) = res2(9);
+}
+
+
+void measuresQuadrangle::getNCJ(const double &x0, const double &y0, const double &z0,
+                                const double &x1, const double &y1, const double &z1,
+                                const double &x2, const double &y2, const double &z2,
+                                const double &x3, const double &y3, const double &z3,
+                                fullVector<double> &ncj)
+{
+  // Compute unit vectors for each edge
+  double x01n, y01n, z01n, x12n, y12n, z12n, x23n, y23n, z23n, x30n, y30n, z30n;
+  unitVec(x0, y0, z0, x1, y1, z1, x01n, y01n, z01n);
+  unitVec(x1, y1, z1, x2, y2, z2, x12n, y12n, z12n);
+  unitVec(x2, y2, z2, x3, y3, z3, x23n, y23n, z23n);
+  unitVec(x3, y3, z3, x0, y0, z0, x30n, y30n, z30n);
+
+  // Compute NCJ at vertex from unit vectors a and b as 0.5*||a^b||/A_triRect
+  // Factor = 1. = 0.5/A_triRect
+  static const double fact = 1.;
+  ncj(0) = triArea(fact, x01n, y01n, z01n, -x30n, -y30n, -z30n);
+  ncj(1) = triArea(fact, x12n, y12n, z12n, -x01n, -y01n, -z01n);
+  ncj(2) = triArea(fact, x23n, y23n, z23n, -x12n, -y12n, -z12n);
+  ncj(3) = triArea(fact, x30n, y30n, z30n, -x23n, -y23n, -z23n);
+}
+
+
+void measuresQuadrangle::getNCJAndGradients(const double &x0, const double &y0, const double &z0,
+                                            const double &x1, const double &y1, const double &z1,
+                                            const double &x2, const double &y2, const double &z2,
+                                            const double &x3, const double &y3, const double &z3,
+                                            fullMatrix<double> &res)
+{
+  // Factor = 1. = 0.5/A_triRect
+  static const double fact = 1.;
+
+  // Compute unit vector and its gradients for each edge
+  fullMatrix<double> vg01(3,7);
+  unitVecAndGrad(x0, y0, z0, x1, y1, z1, vg01);
+  fullMatrix<double> vg12(3,7);
+  unitVecAndGrad(x1, y1, z1, x2, y2, z2, vg12);
+  fullMatrix<double> vg23(3,7);
+  unitVecAndGrad(x2, y2, z2, x3, y3, z3, vg23);
+  fullMatrix<double> vg30(3,7);
+  unitVecAndGrad(x3, y3, z3, x0, y0, z0, vg30);
+  fullMatrix<double> vg10(3,7);
+  revertVG(vg01, vg10);
+  fullMatrix<double> vg21(3,7);
+  revertVG(vg12, vg21);
+  fullMatrix<double> vg32(3,7);
+  revertVG(vg23, vg32);
+  fullMatrix<double> vg03(3,7);
+  revertVG(vg30, vg03);
+
+  // Compute NCJ at vertex 0 as 0.5*||u01^u03||/A_triRect
+  // and gradients w.r.t. x0, y0, z0, x1, y1, z1, x3, y3, z3
+  fullVector<double> res0(10);
+  triAreaAndGrad(fact, vg01, vg03, res0);
+  res(0, 0) = res0(0); res(0, 1) =  res0(1); res(0, 2) =  res0(2);
+  res(0, 3) = res0(3); res(0, 4) =  res0(4); res(0, 5) =  res0(5);
+  res(0, 6) =      0.; res(0, 7) =       0.; res(0, 8) =       0.;
+  res(0, 9) = res0(6); res(0, 10) = res0(7); res(0, 11) = res0(8);
+  res(0, 12) = res0(9);
+
+  // Compute NCJ at vertex 1 as 0.5*||u12^u10||/A_triRect
+  // and gradients w.r.t. x1, y1, z1, x2, y2, z2, x0, y0, z0
+  fullVector<double> res1(10);
+  triAreaAndGrad(fact, vg12, vg10, res1);
+  res(1, 0) = res1(6); res(1, 1) = res1(7); res(1, 2) = res1(8);
+  res(1, 3) = res1(0); res(1, 4) = res1(1); res(1, 5) = res1(2);
+  res(1, 6) = res1(3); res(1, 7) = res1(4); res(1, 8) = res1(5);
+  res(1, 9) =      0.; res(1, 10) =     0.; res(1, 11) =     0.;
+  res(1, 12) = res1(9);
+
+  // Compute NCJ at vertex 2 as 0.5*||u23^u21||/A_triRect
+  // and gradients w.r.t. x2, y2, z2, x3, y3, z3, x1, y1, z1
+  fullVector<double> res2(10);
+  triAreaAndGrad(fact, vg23, vg21, res2);
+  res(2, 0) =      0.; res(2, 1) =       0.; res(2, 2) =       0.;
+  res(2, 3) = res2(6); res(2, 4) =  res2(7); res(2, 5) =  res2(8);
+  res(2, 6) = res2(0); res(2, 7) =  res2(1); res(2, 8) =  res2(2);
+  res(2, 9) = res2(3); res(2, 10) = res2(4); res(2, 11) = res2(5);
+  res(2, 12) = res2(9);
+
+  // Compute NCJ at vertex 3 as 0.5*||u30^u32||/A_triRect
+  // and gradients w.r.t. x3, y3, z3, x0, y0, z0, x2, y2, z2
+  fullVector<double> res3(10);
+  triAreaAndGrad(fact, vg23, vg21, res3);
+  res(3, 0) = res3(3); res(3, 1) =  res3(4); res(3, 2) =  res3(5);
+  res(3, 3) =      0.; res(3, 4) =       0.; res(3, 5) =       0.;
+  res(3, 6) = res3(6); res(3, 7) =  res3(7); res(3, 8) =  res3(8);
+  res(3, 9) = res3(0); res(3, 10) = res3(1); res(3, 11) = res3(2);
+  res(3, 12) = res3(9);
+}

Mesh/qualityMeasures.h

@@ -6,6 +6,8 @@
 #ifndef _QUALITY_MEASURES_H_
 #define _QUALITY_MEASURES_H_
 
+#include "fullMatrix.h"
+
 class BDS_Point;
 class BDS_Face;
 class MVertex;
@@ -41,4 +43,32 @@ double qmTet(const double &x1, const double &y1, const double &z1,
              const double &x4, const double &y4, const double &z4, 
              const qualityMeasure4Tet &cr, double *volume = 0);
 
+class measuresTriangle
+{
+public:
+  static void getNCJ(const double &x0, const double &y0, const double &z0,
+                     const double &x1, const double &y1, const double &z1,
+                     const double &x2, const double &y2, const double &z2,
+                     fullVector<double> &ncj);
+  static void getNCJAndGradients(const double &x0, const double &y0, const double &z0,
+                                 const double &x1, const double &y1, const double &z1,
+                                 const double &x2, const double &y2, const double &z2,
+                                 fullMatrix<double> &ncj);
+};
+
+class measuresQuadrangle
+{
+public:
+  static void getNCJ(const double &x0, const double &y0, const double &z0,
+                     const double &x1, const double &y1, const double &z1,
+                     const double &x2, const double &y2, const double &z2,
+                     const double &x3, const double &y3, const double &z3,
+                     fullVector<double> &ncj);
+  static void getNCJAndGradients(const double &x0, const double &y0, const double &z0,
+                                 const double &x1, const double &y1, const double &z1,
+                                 const double &x2, const double &y2, const double &z2,
+                                 const double &x3, const double &y3, const double &z3,
+                                 fullMatrix<double> &ncj);
+};
+
 #endif