added numerical quadrature for non planar cells

diff --git a/source/grid/tria_accessor.cc b/source/grid/tria_accessor.cc
index a7edc8283da75f8743dfe76689d13667311d35e1..09015c006a71ab184a89a279f278f8623f8d73a8 100644 (file)
--- a/source/grid/tria_accessor.cc
+++ b/source/grid/tria_accessor.cc
@@ -1285,41 +1285,125 @@ namespace
   }
 
 
+  // a 2d face in 3d space
+
   // a 2d face in 3d space
   template <int dim>
   double
-  measure(const dealii::TriaAccessor<2, dim, 3> &accessor)
+  measure(const TriaAccessor<2, dim, 3> &accessor)
   {
-    // In general the area can be computed as
-    // the integral of the cross product of the two tangential vectors
-
-    // If we assume a bilinear patch parametrized in u and v we get that
-    // t_u = (v_1 - v_0) + v (v_3 - v_2 - v_1 + v_0)
-    // t_v = (v_2 - v_0) + u (v_3 - v_2 - v_1 + v_0)
-    // So t_u x t_v = (v_1 - v_0) x (v_2 - v_0) + u (v_1 - v_0) x (v_3 - v_2 -
-    // v_1 + v_0) + v (v_3 - v_2 - v_1 + v_0) x (v_2 - v_0) t_u x t_v = w_1 + u
-    // w_2 + v w_3 we can integrate the square norm (t_u x t_v) * (t_u x t_v) =
-    // w_1*w_1 + u^2 w_2*w_2 + v^2 w_3*w_3 + 2u w_1*w_2 + 2v w_1*w_3 + 2uv
-    // w_2*w_3 in u and v getting (between zero and one) w_1*w_1 + 1/3 w_2*w_2 +
-    // 1/3 w_3*w_3 + w_1*w_2 + w_1*w_3 + 1/2 w_2*w_3
-
-    const Tensor<1, 3> w_1 =
-      cross_product_3d(accessor.vertex(1) - accessor.vertex(0),
-                       accessor.vertex(2) - accessor.vertex(0));
-    const Tensor<1, 3> w_2 =
-      cross_product_3d(accessor.vertex(1) - accessor.vertex(0),
-                       accessor.vertex(3) - accessor.vertex(2) -
-                         accessor.vertex(1) + accessor.vertex(0));
-    const Tensor<1, 3> w_3 =
-      cross_product_3d(accessor.vertex(3) - accessor.vertex(2) -
-                         accessor.vertex(1) + accessor.vertex(0),
-                       accessor.vertex(2) - accessor.vertex(0));
-
-
-    return std::sqrt(scalar_product(w_1, w_1) + scalar_product(w_1, w_2) +
-                     scalar_product(w_1, w_3) + 0.5 * scalar_product(w_2, w_3) +
-                     1. / 3 * scalar_product(w_2, w_2) +
-                     1. / 3 * scalar_product(w_3, w_3));
+    // If the face is planar, the diagonal from vertex 0 to vertex 3,
+    // v_03, should be in the plane P_012 of vertices 0, 1 and 2.  Get
+    // the normal vector of P_012 and test if v_03 is orthogonal to
+    // that. If so, the face is planar and computing its area is simple.
+    const Tensor<1, 3> v01 = accessor.vertex(1) - accessor.vertex(0);
+    const Tensor<1, 3> v02 = accessor.vertex(2) - accessor.vertex(0);
+
+    const Tensor<1, 3> normal = cross_product_3d(v01, v02);
+
+    const Tensor<1, 3> v03 = accessor.vertex(3) - accessor.vertex(0);
+
+    // check whether v03 does not lie in the plane of v01 and v02
+    // (i.e., whether the face is not planar). we do so by checking
+    // whether the triple product (v01 x v02) * v03 forms a positive
+    // volume relative to |v01|*|v02|*|v03|. the test checks the
+    // squares of these to avoid taking norms/square roots:
+    if (std::abs((v03 * normal) * (v03 * normal) /
+                 ((v03 * v03) * (v01 * v01) * (v02 * v02))) >= 1e-24)
+      {
+        // If the vectors are non planar we integrate the norm of the normal
+        // vector using a numerical Gauss scheme of order 4.
+        const Tensor<1, 3> w_1 =
+          cross_product_3d(accessor.vertex(1) - accessor.vertex(0),
+                           accessor.vertex(2) - accessor.vertex(0));
+        const Tensor<1, 3> w_2 =
+          cross_product_3d(accessor.vertex(1) - accessor.vertex(0),
+                           accessor.vertex(3) - accessor.vertex(2) -
+                             accessor.vertex(1) + accessor.vertex(0));
+        const Tensor<1, 3> w_3 =
+          cross_product_3d(accessor.vertex(3) - accessor.vertex(2) -
+                             accessor.vertex(1) + accessor.vertex(0),
+                           accessor.vertex(2) - accessor.vertex(0));
+
+        double a = scalar_product(w_1, w_1);
+        double b = scalar_product(w_2, w_2);
+        double c = scalar_product(w_3, w_3);
+        double d = scalar_product(w_1, w_2);
+        double e = scalar_product(w_1, w_3);
+        double f = scalar_product(w_2, w_3);
+
+        return 0.03025074832140047 *
+                 std::sqrt(a + 0.0048207809894260144 * b +
+                           0.0048207809894260144 * c + 0.13886368840594743 * d +
+                           0.13886368840594743 * e +
+                           0.0096415619788520288 * f) +
+               0.056712962962962937 *
+                 std::sqrt(a + 0.0048207809894260144 * b +
+                           0.10890625570683385 * c + 0.13886368840594743 * d +
+                           0.66001895641514374 * e + 0.045826333352825557 * f) +
+               0.056712962962962937 *
+                 std::sqrt(a + 0.0048207809894260144 * b +
+                           0.44888729929169013 * c + 0.13886368840594743 * d +
+                           1.3399810435848563 * e + 0.09303735505312187 * f) +
+               0.03025074832140047 *
+                 std::sqrt(a + 0.0048207809894260144 * b +
+                           0.86595709258347853 * c + 0.13886368840594743 * d +
+                           1.8611363115940525 * e + 0.12922212642709538 * f) +
+               0.056712962962962937 *
+                 std::sqrt(a + 0.10890625570683385 * b +
+                           0.0048207809894260144 * c + 0.66001895641514374 * d +
+                           0.13886368840594743 * e + 0.045826333352825557 * f) +
+               0.10632332575267359 *
+                 std::sqrt(a + 0.10890625570683385 * b +
+                           0.10890625570683385 * c + 0.66001895641514374 * d +
+                           0.66001895641514374 * e + 0.2178125114136677 * f) +
+               0.10632332575267359 *
+                 std::sqrt(a + 0.10890625570683385 * b +
+                           0.44888729929169013 * c + 0.66001895641514374 * d +
+                           1.3399810435848563 * e + 0.44220644500147605 * f) +
+               0.056712962962962937 *
+                 std::sqrt(a + 0.10890625570683385 * b +
+                           0.86595709258347853 * c + 0.66001895641514374 * d +
+                           1.8611363115940525 * e + 0.61419262306231814 * f) +
+               0.056712962962962937 *
+                 std::sqrt(a + 0.44888729929169013 * b +
+                           0.0048207809894260144 * c + 1.3399810435848563 * d +
+                           0.13886368840594743 * e + 0.09303735505312187 * f) +
+               0.10632332575267359 *
+                 std::sqrt(a + 0.44888729929169013 * b +
+                           0.10890625570683385 * c + 1.3399810435848563 * d +
+                           0.66001895641514374 * e + 0.44220644500147605 * f) +
+               0.10632332575267359 *
+                 std::sqrt(a + 0.44888729929169013 * b +
+                           0.44888729929169013 * c + 1.3399810435848563 * d +
+                           1.3399810435848563 * e + 0.89777459858338027 * f) +
+               0.056712962962962937 *
+                 std::sqrt(a + 0.44888729929169013 * b +
+                           0.86595709258347853 * c + 1.3399810435848563 * d +
+                           1.8611363115940525 * e + 1.2469436885317342 * f) +
+               0.03025074832140047 *
+                 std::sqrt(a + 0.86595709258347853 * b +
+                           0.0048207809894260144 * c + 1.8611363115940525 * d +
+                           0.13886368840594743 * e + 0.12922212642709538 * f) +
+               0.056712962962962937 *
+                 std::sqrt(a + 0.86595709258347853 * b +
+                           0.10890625570683385 * c + 1.8611363115940525 * d +
+                           0.66001895641514374 * e + 0.61419262306231814 * f) +
+               0.056712962962962937 *
+                 std::sqrt(a + 0.86595709258347853 * b +
+                           0.44888729929169013 * c + 1.8611363115940525 * d +
+                           1.3399810435848563 * e + 1.2469436885317342 * f) +
+               0.03025074832140047 *
+                 std::sqrt(a + 0.86595709258347853 * b +
+                           0.86595709258347853 * c + 1.8611363115940525 * d +
+                           1.8611363115940525 * e + 1.7319141851669571 * f);
+      }
+
+    // the face is planar. then its area is 1/2 of the norm of the
+    // cross product of the two diagonals
+    const Tensor<1, 3> v12        = accessor.vertex(2) - accessor.vertex(1);
+    const Tensor<1, 3> twice_area = cross_product_3d(v03, v12);
+    return 0.5 * twice_area.norm();
   }
 
 

diff --git a/tests/codim_one/measure_of_distorted_codim_one_face.cc b/tests/codim_one/measure_of_distorted_codim_one_face.cc
index efb6d30c990ce345bc154c161eaae33daf30de06..342a64f10026cd710a4b1f2d587d37112bafb96d 100644 (file)
--- a/tests/codim_one/measure_of_distorted_codim_one_face.cc
+++ b/tests/codim_one/measure_of_distorted_codim_one_face.cc
@@ -43,10 +43,14 @@ test()
   std::vector<CellData<2>> cells;
   SubCellData              subcelldata;
 
-  double tol = 1e-15;
-  vertices.push_back(Point<3>{10, 56, 0});
-  vertices.push_back(Point<3>{22, 1, 0});
-  vertices.push_back(Point<3>{15, 44, 0});
+  double tol = 1e-12;
+  // vertices.push_back(Point<3>{10, 56, 0});
+  // vertices.push_back(Point<3>{22, 1, 0});
+  // vertices.push_back(Point<3>{15, 44, 0});
+  // vertices.push_back(Point<3>{1, 1, 1});
+  vertices.push_back(Point<3>{0, 0, 1});
+  vertices.push_back(Point<3>{1, 0, -10});
+  vertices.push_back(Point<3>{0, 1, -1});
   vertices.push_back(Point<3>{1, 1, 1});
 
   cells.resize(1);
@@ -62,7 +66,7 @@ test()
   DoFHandler<2, 3> dof_handler(tria);
   dof_handler.distribute_dofs(fe);
 
-  QGauss<2>       quadrature_formula(1);
+  QGauss<2>       quadrature_formula(4);
   MappingQ1<2, 3> mapping;
   FEValues<2, 3>  fe_values(mapping, fe, quadrature_formula, update_JxW_values);
 
@@ -90,7 +94,7 @@ int
 main()
 {
   initlog();
-  deallog << std::setprecision(5);
+  deallog << std::setprecision(8);
 
 
   test();