added measure computation for distorted 2d cells

diff --git a/source/grid/tria_accessor.cc b/source/grid/tria_accessor.cc
index 290917877cd4cddc5280a3333fe5e2aa0d0ad9a9..62919885f65dfbfe9abc71402dc1e72668d9a49f 100644 (file)
--- a/source/grid/tria_accessor.cc
+++ b/source/grid/tria_accessor.cc
@@ -1286,42 +1286,35 @@ namespace
 
 
   // a 2d face in 3d space
+  template <int dim>
   double
-  measure(const dealii::TriaAccessor<2, 3, 3> &accessor)
+  measure(const dealii::TriaAccessor<2, dim, 3> &accessor)
   {
-    // 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)
-      {
-        Assert(
-          false,
-          ExcMessage(
-            "Computing the measure of a nonplanar face is not implemented!"));
-        return std::numeric_limits<double>::quiet_NaN();
-      }
+    // In general the area can be computed as
+    // 0.25*(v_0+v_1-v_2-v_3)*(v_0-v_1+v_2-v_3)
 
-    // 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();
+    const Tensor<1, 3> piece_1 = accessor.vertex(0) + accessor.vertex(1) -
+                                 accessor.vertex(2) - accessor.vertex(3);
+    const Tensor<1, 3> piece_2 = accessor.vertex(0) - accessor.vertex(1) +
+                                 accessor.vertex(2) - accessor.vertex(3);
+
+    return 0.25 * cross_product_3d(piece_1, piece_2).norm();
   }
 
+  // // a 2d cell in 3d space
+  // double
+  // measure(const dealii::TriaAccessor<2, 2, 3> &accessor)
+  // {
+  //   // In general the area can be computed as
+  //   // 0.25*(v_0+v_1-v_2-v_3)*(v_0-v_1+v_2-v_3)
+
+  //   const Tensor<1, 3> piece_1 = accessor.vertex(0) + accessor.vertex(1) -
+  //                                accessor.vertex(2) - accessor.vertex(3);
+  //   const Tensor<1, 3> piece_2 = accessor.vertex(0) - accessor.vertex(1) +
+  //                                accessor.vertex(2) - accessor.vertex(3);
+  //   return 0.25 * cross_product_3d(piece_1, piece_2).norm();
+  // }
+
 
 
   template <int structdim, int dim, int spacedim>

diff --git a/tests/codim_one/measure_of_distorted_codim_one_face.cc b/tests/codim_one/measure_of_distorted_codim_one_face.cc
new file mode 100644 (file)
index 0000000..efb6d30
--- /dev/null
+++ b/tests/codim_one/measure_of_distorted_codim_one_face.cc
@@ -0,0 +1,97 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2014 - 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// Test by Nicola Giuliani: compute the measure of a 2d distorted cell in codim
+// one
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include "../tests.h"
+
+
+
+void
+test()
+{
+  Triangulation<2, 3>      tria;
+  std::vector<Point<3>>    vertices;
+  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});
+  vertices.push_back(Point<3>{1, 1, 1});
+
+  cells.resize(1);
+
+  cells[0].vertices[0] = 0;
+  cells[0].vertices[1] = 1;
+  cells[0].vertices[2] = 2;
+  cells[0].vertices[3] = 3;
+
+  tria.create_triangulation(vertices, cells, subcelldata);
+
+  FE_Q<2, 3>       fe(1);
+  DoFHandler<2, 3> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+
+  QGauss<2>       quadrature_formula(1);
+  MappingQ1<2, 3> mapping;
+  FEValues<2, 3>  fe_values(mapping, fe, quadrature_formula, update_JxW_values);
+
+  double sum_1    = 0;
+  double sum_2    = 0;
+  auto   cell     = dof_handler.begin_active();
+  auto   accessor = tria.begin_active();
+
+  sum_1 = cell->measure();
+  fe_values.reinit(cell);
+  for (unsigned int q = 0; q < quadrature_formula.size(); ++q)
+    {
+      sum_2 += fe_values.JxW(q);
+    }
+  if (std::abs(sum_1 - sum_2) < tol)
+    deallog << "Test passed" << std::endl;
+  else
+    deallog << sum_1 << " " << sum_2 << std::endl;
+  deallog << std::endl;
+}
+
+
+
+int
+main()
+{
+  initlog();
+  deallog << std::setprecision(5);
+
+
+  test();
+}

diff --git a/tests/codim_one/measure_of_distorted_codim_one_face.output b/tests/codim_one/measure_of_distorted_codim_one_face.output
new file mode 100644 (file)
index 0000000..7ca18d1
--- /dev/null
+++ b/tests/codim_one/measure_of_distorted_codim_one_face.output
@@ -0,0 +1,3 @@
+
+DEAL::Test passed
+DEAL::