Fix FlatManifold::normal_vector for triangular faces

diff --git a/include/deal.II/grid/reference_cell.h b/include/deal.II/grid/reference_cell.h
index 9b99f6ef7624accec1d09a3a7b99e38600722ca6..28a314d2b4700d04775374dfdfb22707bb39559d 100644 (file)
--- a/include/deal.II/grid/reference_cell.h
+++ b/include/deal.II/grid/reference_cell.h
@@ -118,6 +118,71 @@ namespace ReferenceCell
     return table[dim][n_vertices];
   }
 
+
+  /**
+   * Compute the value of the $i$-th linear shape function at location $\xi$ for
+   * a given reference-cell type.
+   */
+  template <int dim>
+  inline double
+  d_linear_shape_function(const Type &       reference_cell,
+                          const Point<dim> & xi,
+                          const unsigned int i)
+  {
+    if (reference_cell == get_hypercube(dim))
+      return GeometryInfo<dim>::d_linear_shape_function(xi, i);
+
+    if (reference_cell ==
+        Type::Tri) // see also Simplex::ScalarPolynomial::compute_value
+      {
+        switch (i)
+          {
+            case 0:
+              return 1.0 - xi[0] - xi[1];
+            case 1:
+              return xi[0];
+            case 2:
+              return xi[1];
+          }
+      }
+
+    Assert(false, ExcNotImplemented());
+
+    return 0.0;
+  }
+
+  /**
+   * Compute the gradient of the $i$-th linear shape function at location $\xi$
+   * for a given reference-cell type.
+   */
+  template <int dim>
+  inline Tensor<1, dim>
+  d_linear_shape_function_gradient(const Type &       reference_cell,
+                                   const Point<dim> & xi,
+                                   const unsigned int i)
+  {
+    if (reference_cell == get_hypercube(dim))
+      return GeometryInfo<dim>::d_linear_shape_function_gradient(xi, i);
+
+    if (reference_cell ==
+        Type::Tri) // see also Simplex::ScalarPolynomial::compute_grad
+      {
+        switch (i)
+          {
+            case 0:
+              return Point<dim>(-1.0, -1.0);
+            case 1:
+              return Point<dim>(+1.0, +0.0);
+            case 2:
+              return Point<dim>(+0.0, +1.0);
+          }
+      }
+
+    Assert(false, ExcNotImplemented());
+
+    return Point<dim>(+0.0, +0.0, +0.0);
+  }
+
   namespace internal
   {
     /**

diff --git a/source/grid/manifold.cc b/source/grid/manifold.cc
index b839f903fdd3ac350ff99422aaafb0dd233f1010..b834e3f6864d1e84b817e18fdb60eff0b63b4e63 100644 (file)
--- a/source/grid/manifold.cc
+++ b/source/grid/manifold.cc
@@ -892,8 +892,19 @@ FlatManifold<dim, spacedim>::normal_vector(
   const unsigned int facedim = dim - 1;
 
   Point<facedim> xi;
-  for (unsigned int i = 0; i < facedim; ++i)
-    xi[i] = 1. / 2;
+
+  const auto face_reference_cell_type = face->reference_cell_type();
+
+  if (face_reference_cell_type == ReferenceCell::get_hypercube(facedim))
+    {
+      for (unsigned int i = 0; i < facedim; ++i)
+        xi[i] = 1. / 2;
+    }
+  else
+    {
+      for (unsigned int i = 0; i < facedim; ++i)
+        xi[i] = 1. / 3;
+    }
 
   const double        eps = 1e-12;
   Tensor<1, spacedim> grad_F[facedim];
@@ -901,17 +912,17 @@ FlatManifold<dim, spacedim>::normal_vector(
   while (true)
     {
       Point<spacedim> F;
-      for (const unsigned int v : GeometryInfo<facedim>::vertex_indices())
-        F += face->vertex(v) *
-             GeometryInfo<facedim>::d_linear_shape_function(xi, v);
+      for (const unsigned int v : face->vertex_indices())
+        F += face->vertex(v) * ReferenceCell::d_linear_shape_function(
+                                 face_reference_cell_type, xi, v);
 
       for (unsigned int i = 0; i < facedim; ++i)
         {
           grad_F[i] = 0;
-          for (const unsigned int v : GeometryInfo<facedim>::vertex_indices())
+          for (const unsigned int v : face->vertex_indices())
             grad_F[i] +=
-              face->vertex(v) *
-              GeometryInfo<facedim>::d_linear_shape_function_gradient(xi, v)[i];
+              face->vertex(v) * ReferenceCell::d_linear_shape_function_gradient(
+                                  face_reference_cell_type, xi, v)[i];
         }
 
       Tensor<1, facedim> J;