]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Introduce ReferenceCell::unit_tangential_vectors() 11423/head
authorPeter Munch <peterrmuench@gmail.com>
Sun, 27 Dec 2020 06:34:18 +0000 (07:34 +0100)
committerPeter Munch <peterrmuench@gmail.com>
Wed, 30 Dec 2020 05:30:35 +0000 (06:30 +0100)
include/deal.II/grid/reference_cell.h
source/fe/mapping_fe.cc
tests/simplex/unit_tangential_vectors_01.cc [new file with mode: 0644]
tests/simplex/unit_tangential_vectors_01.output [new file with mode: 0644]

index 2d3d4580b1022392159845d9c2be9b6dd4349b04..2bfa0943b121fe729a3fccaa587c60a57efb49c8 100644 (file)
@@ -45,6 +45,31 @@ namespace ReferenceCell
     Invalid = static_cast<std::uint8_t>(-1)
   };
 
+  /**
+   * Return the dimension of the given reference-cell type @p type.
+   */
+  inline unsigned int
+  get_dimension(const Type &type)
+  {
+    switch (type)
+      {
+        case Type::Vertex:
+          return 0;
+        case Type::Line:
+          return 1;
+        case Type::Tri:
+        case Type::Quad:
+          return 2;
+        case Type::Tet:
+        case Type::Pyramid:
+        case Type::Wedge:
+        case Type::Hex:
+          return 3;
+        default:
+          return numbers::invalid_unsigned_int;
+      }
+  }
+
   /**
    * Return the correct simplex reference cell type for the given dimension
    * @p dim.
@@ -374,7 +399,85 @@ namespace ReferenceCell
     return temp_;
   }
 
+  /*
+   * Return i-th unit tangential vector of a face of the reference cell.
+   * The vectors are arranged such that the
+   * cross product between the two vectors returns the unit normal vector.
+   */
+  template <int dim>
+  inline Tensor<1, dim>
+  unit_tangential_vectors(const Type &       reference_cell,
+                          const unsigned int face_no,
+                          const unsigned int i)
+  {
+    AssertDimension(dim, get_dimension(reference_cell));
+    AssertIndexRange(i, dim - 1);
 
+    if (reference_cell == get_hypercube(dim))
+      {
+        AssertIndexRange(face_no, GeometryInfo<dim>::faces_per_cell);
+        return GeometryInfo<dim>::unit_tangential_vectors[face_no][i];
+      }
+    else if (reference_cell == Type::Tri)
+      {
+        AssertIndexRange(face_no, 3);
+        static const std::array<Tensor<1, dim>, 3> table = {
+          {Point<dim>(1, 0),
+           Point<dim>(-std::sqrt(0.5), +std::sqrt(0.5)),
+           Point<dim>(0, -1)}};
+
+        return table[face_no];
+      }
+    else if (reference_cell == Type::Tet)
+      {
+        AssertIndexRange(face_no, 4);
+        static const std::array<std::array<Tensor<1, dim>, 2>, 4> table = {
+          {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
+           {{Point<dim>(1, 0, 0), Point<dim>(0, 0, 1)}},
+           {{Point<dim>(0, 0, 1), Point<dim>(0, 1, 0)}},
+           {{Point<dim>(-std::pow(1.0 / 3.0, 1.0 / 4.0),
+                        +std::pow(1.0 / 3.0, 1.0 / 4.0),
+                        0),
+             Point<dim>(-std::pow(1.0 / 3.0, 1.0 / 4.0),
+                        0,
+                        +std::pow(1.0 / 3.0, 1.0 / 4.0))}}}};
+
+        return table[face_no][i];
+      }
+    else if (reference_cell == Type::Wedge)
+      {
+        AssertIndexRange(face_no, 5);
+        static const std::array<std::array<Tensor<1, dim>, 2>, 5> table = {
+          {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
+           {{Point<dim>(1, 0, 0), Point<dim>(0, 0, 1)}},
+           {{Point<dim>(-1 / std::sqrt(2.0), +1 / std::sqrt(2.0), 0),
+             Point<dim>(0, 0, 1)}},
+           {{Point<dim>(0, 0, 1), Point<dim>(0, 1, 0)}},
+           {{Point<dim>(1, 0, 0), Point<dim>(0, 0, 1)}}}};
+
+        return table[face_no][i];
+      }
+    else if (reference_cell == Type::Pyramid)
+      {
+        AssertIndexRange(face_no, 5);
+        static const std::array<std::array<Tensor<1, dim>, 2>, 5> table = {
+          {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
+           {{Point<dim>(+1.0 / sqrt(2.0), 0, +1.0 / sqrt(2.0)),
+             Point<dim>(0, 1, 0)}},
+           {{Point<dim>(+1.0 / sqrt(2.0), 0, -1.0 / sqrt(2.0)),
+             Point<dim>(0, 1, 0)}},
+           {{Point<dim>(1, 0, 0),
+             Point<dim>(0, +1.0 / sqrt(2.0), +1.0 / sqrt(2.0))}},
+           {{Point<dim>(1, 0, 0),
+             Point<dim>(0, +1.0 / sqrt(2.0), -1.0 / sqrt(2.0))}}}};
+
+        return table[face_no][i];
+      }
+
+    Assert(false, ExcNotImplemented());
+
+    return {};
+  }
 
   namespace internal
   {
index 7921ebfd074fff0fd34be600cbfbd7e0ce0ec110..c004084ca2b4996240fd75092ceea58c17c4164c 100644 (file)
@@ -151,259 +151,26 @@ MappingFE<dim, spacedim>::InternalData::initialize_face(
                  std::vector<Tensor<1, spacedim>>(n_original_q_points));
 
       // Compute tangentials to the unit cell.
-      if (this->fe.reference_cell_type() == ReferenceCell::get_hypercube(dim))
+      const auto reference_cell_type = this->fe.reference_cell_type();
+      const auto n_faces =
+        ReferenceCell::internal::Info::get_cell(reference_cell_type).n_faces();
+
+      for (unsigned int i = 0; i < n_faces; ++i)
         {
-          for (const unsigned int i : GeometryInfo<dim>::face_indices())
+          unit_tangentials[i].resize(n_original_q_points);
+          std::fill(unit_tangentials[i].begin(),
+                    unit_tangentials[i].end(),
+                    ReferenceCell::unit_tangential_vectors<dim>(
+                      reference_cell_type, i, 0));
+          if (dim > 2)
             {
-              unit_tangentials[i].resize(n_original_q_points);
-              std::fill(unit_tangentials[i].begin(),
-                        unit_tangentials[i].end(),
-                        GeometryInfo<dim>::unit_tangential_vectors[i][0]);
-              if (dim > 2)
-                {
-                  unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
-                    .resize(n_original_q_points);
-                  std::fill(
-                    unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
-                      .begin(),
-                    unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
-                      .end(),
-                    GeometryInfo<dim>::unit_tangential_vectors[i][1]);
-                }
+              unit_tangentials[n_faces + i].resize(n_original_q_points);
+              std::fill(unit_tangentials[n_faces + i].begin(),
+                        unit_tangentials[n_faces + i].end(),
+                        ReferenceCell::unit_tangential_vectors<dim>(
+                          reference_cell_type, i, 1));
             }
         }
-      else if (this->fe.reference_cell_type() == ReferenceCell::Type::Tri)
-        {
-          Tensor<1, dim> t1;
-          constexpr int  d0 = 0;
-          constexpr int  d1 = 1 % dim;
-
-          t1[d0] = 1;
-          t1[d1] = 0;
-          for (unsigned int i = 0; i < n_original_q_points; i++)
-            unit_tangentials[0].emplace_back(t1);
-          t1[d0] = -std::sqrt(0.5);
-          t1[d1] = +std::sqrt(0.5);
-          for (unsigned int i = 0; i < n_original_q_points; i++)
-            unit_tangentials[1].emplace_back(t1);
-          t1[d0] = 0;
-          t1[d1] = -1;
-          for (unsigned int i = 0; i < n_original_q_points; i++)
-            unit_tangentials[2].emplace_back(t1);
-        }
-      else if (this->fe.reference_cell_type() == ReferenceCell::Type::Tet)
-        {
-          Tensor<1, dim> t1;
-          constexpr int  d0 = 0;
-          constexpr int  d1 = 1 % dim;
-          constexpr int  d2 = 2 % dim;
-
-          t1[d0] = 0;
-          t1[d1] = 1;
-          t1[d2] = 0; // face 0
-          for (unsigned int i = 0; i < n_original_q_points; i++)
-            unit_tangentials[0].emplace_back(t1);
-
-          t1[d0] = 1;
-          t1[d1] = 0;
-          t1[d2] = 0; // face 0
-          for (unsigned int i = 0; i < n_original_q_points; i++)
-            unit_tangentials[4].emplace_back(t1);
-
-          t1[d0] = 1;
-          t1[d1] = 0;
-          t1[d2] = 0; // face 1
-          for (unsigned int i = 0; i < n_original_q_points; i++)
-            unit_tangentials[1].emplace_back(t1);
-
-          t1[d0] = 0;
-          t1[d1] = 0;
-          t1[d2] = 1; // face 1
-          for (unsigned int i = 0; i < n_original_q_points; i++)
-            unit_tangentials[5].emplace_back(t1);
-
-          t1[d0] = 0;
-          t1[d1] = 0;
-          t1[d2] = 1; // face 2
-          for (unsigned int i = 0; i < n_original_q_points; i++)
-            unit_tangentials[2].emplace_back(t1);
-
-          t1[d0] = 0;
-          t1[d1] = 1;
-          t1[d2] = 0; // face 2
-          for (unsigned int i = 0; i < n_original_q_points; i++)
-            unit_tangentials[6].emplace_back(t1);
-
-          t1[d0] = -std::pow(1.0 / 3.0, 1.0 / 4.0);
-          t1[d1] = +std::pow(1.0 / 3.0, 1.0 / 4.0);
-          t1[d2] = +0; // face 3
-          for (unsigned int i = 0; i < n_original_q_points; i++)
-            unit_tangentials[3].emplace_back(t1);
-
-          t1[d0] = -std::pow(1.0 / 3.0, 1.0 / 4.0);
-          t1[d1] = +0;
-          t1[d2] = +std::pow(1.0 / 3.0, 1.0 / 4.0); // face 3
-          for (unsigned int i = 0; i < n_original_q_points; i++)
-            unit_tangentials[7].emplace_back(t1);
-        }
-      else if (this->fe.reference_cell_type() == ReferenceCell::Type::Wedge)
-        {
-          Tensor<1, dim> t1;
-          constexpr int  d0 = 0;
-          constexpr int  d1 = 1 % dim;
-          constexpr int  d2 = 2 % dim;
-
-          // face 0
-          t1[d0] = 0;
-          t1[d1] = 1;
-          t1[d2] = 0;
-          for (unsigned int i = 0; i < n_original_q_points; i++)
-            unit_tangentials[0].emplace_back(t1);
-
-          // face 0
-          t1[d0] = 1;
-          t1[d1] = 0;
-          t1[d2] = 0;
-          for (unsigned int i = 0; i < n_original_q_points; i++)
-            unit_tangentials[5].emplace_back(t1);
-
-          // face 1
-          t1[d0] = 1;
-          t1[d1] = 0;
-          t1[d2] = 0;
-          for (unsigned int i = 0; i < n_original_q_points; i++)
-            unit_tangentials[1].emplace_back(t1);
-
-          // face 1
-          t1[d0] = 0;
-          t1[d1] = 1;
-          t1[d2] = 0;
-          for (unsigned int i = 0; i < n_original_q_points; i++)
-            unit_tangentials[6].emplace_back(t1);
-
-          // face 2
-          t1[d0] = 1;
-          t1[d1] = 0;
-          t1[d2] = 0;
-          for (unsigned int i = 0; i < n_original_q_points; i++)
-            unit_tangentials[2].emplace_back(t1);
-
-          // face 2
-          t1[d0] = 0;
-          t1[d1] = 0;
-          t1[d2] = 1;
-          for (unsigned int i = 0; i < n_original_q_points; i++)
-            unit_tangentials[7].emplace_back(t1);
-
-          // face 3
-          t1[d0] = -1 / std::sqrt(2.0);
-          t1[d1] = +1 / std::sqrt(2.0);
-          t1[d2] = 0;
-          for (unsigned int i = 0; i < n_original_q_points; i++)
-            unit_tangentials[3].emplace_back(t1);
-
-          // face 3
-          t1[d0] = 0;
-          t1[d1] = 0;
-          t1[d2] = 1;
-          for (unsigned int i = 0; i < n_original_q_points; i++)
-            unit_tangentials[8].emplace_back(t1);
-
-          // face 4
-          t1[d0] = +0;
-          t1[d1] = +0;
-          t1[d2] = +1;
-          for (unsigned int i = 0; i < n_original_q_points; i++)
-            unit_tangentials[4].emplace_back(t1);
-
-          // face 4
-          t1[d0] = +0;
-          t1[d1] = +1;
-          t1[d2] = +0;
-          for (unsigned int i = 0; i < n_original_q_points; i++)
-            unit_tangentials[9].emplace_back(t1);
-        }
-      else if (this->fe.reference_cell_type() == ReferenceCell::Type::Pyramid)
-        {
-          Tensor<1, dim> t1;
-          constexpr int  d0 = 0;
-          constexpr int  d1 = 1 % dim;
-          constexpr int  d2 = 2 % dim;
-
-          // face 0
-          t1[d0] = 0;
-          t1[d1] = 1;
-          t1[d2] = 0;
-          for (unsigned int i = 0; i < n_original_q_points; i++)
-            unit_tangentials[0].emplace_back(t1);
-
-          // face 0
-          t1[d0] = 1;
-          t1[d1] = 0;
-          t1[d2] = 0;
-          for (unsigned int i = 0; i < n_original_q_points; i++)
-            unit_tangentials[5].emplace_back(t1);
-
-          // face 1
-          t1[d0] = 1.0 / sqrt(2.0);
-          t1[d1] = 0;
-          t1[d2] = 1.0 / sqrt(2.0);
-          for (unsigned int i = 0; i < n_original_q_points; i++)
-            unit_tangentials[1].emplace_back(t1);
-
-          // face 1
-          t1[d0] = 0;
-          t1[d1] = 1;
-          t1[d2] = 0;
-          for (unsigned int i = 0; i < n_original_q_points; i++)
-            unit_tangentials[6].emplace_back(t1);
-
-          // face 2
-          t1[d0] = 1.0 / sqrt(2.0);
-          t1[d1] = 0;
-          t1[d2] = -1.0 / sqrt(2.0);
-          for (unsigned int i = 0; i < n_original_q_points; i++)
-            unit_tangentials[2].emplace_back(t1);
-
-          // face 2
-          t1[d0] = 0;
-          t1[d1] = 1;
-          t1[d2] = 0;
-          for (unsigned int i = 0; i < n_original_q_points; i++)
-            unit_tangentials[7].emplace_back(t1);
-
-          // face 3
-          t1[d0] = 1;
-          t1[d1] = 0;
-          t1[d2] = 0;
-          for (unsigned int i = 0; i < n_original_q_points; i++)
-            unit_tangentials[3].emplace_back(t1);
-
-          // face 3
-          t1[d0] = 0;
-          t1[d1] = 1.0 / sqrt(2.0);
-          t1[d2] = 1.0 / sqrt(2.0);
-          for (unsigned int i = 0; i < n_original_q_points; i++)
-            unit_tangentials[8].emplace_back(t1);
-
-          // face 4
-          t1[d0] = 1;
-          t1[d1] = 0;
-          t1[d2] = 0;
-          for (unsigned int i = 0; i < n_original_q_points; i++)
-            unit_tangentials[4].emplace_back(t1);
-
-          // face 4
-          t1[d0] = 0;
-          t1[d1] = +1.0 / sqrt(2.0);
-          t1[d2] = -1.0 / sqrt(2.0);
-          for (unsigned int i = 0; i < n_original_q_points; i++)
-            unit_tangentials[9].emplace_back(t1);
-        }
-      else
-        {
-          Assert(false, ExcNotImplemented());
-        }
     }
 }
 
diff --git a/tests/simplex/unit_tangential_vectors_01.cc b/tests/simplex/unit_tangential_vectors_01.cc
new file mode 100644 (file)
index 0000000..660cb58
--- /dev/null
@@ -0,0 +1,55 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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 ReferenceCell::unit_tangential_vectors() for all reference-cell types.
+
+
+#include <deal.II/grid/reference_cell.h>
+
+#include <deal.II/simplex/quadrature_lib.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim>
+void
+test(const ReferenceCell::Type &reference_cell)
+{
+  for (const auto face_no :
+       ReferenceCell::internal::Info::get_cell(reference_cell).face_indices())
+    for (unsigned int i = 0; i < dim - 1; ++i)
+      deallog << ReferenceCell::unit_tangential_vectors<dim>(reference_cell,
+                                                             face_no,
+                                                             i)
+              << std::endl;
+  deallog << std::endl;
+}
+
+
+
+int
+main()
+{
+  initlog();
+
+  test<2>(ReferenceCell::Type::Tri);
+  test<2>(ReferenceCell::Type::Quad);
+  test<3>(ReferenceCell::Type::Tet);
+  test<3>(ReferenceCell::Type::Pyramid);
+  test<3>(ReferenceCell::Type::Wedge);
+  test<3>(ReferenceCell::Type::Hex);
+}
diff --git a/tests/simplex/unit_tangential_vectors_01.output b/tests/simplex/unit_tangential_vectors_01.output
new file mode 100644 (file)
index 0000000..46a477e
--- /dev/null
@@ -0,0 +1,54 @@
+
+DEAL::1.00000 0.00000
+DEAL::-0.707107 0.707107
+DEAL::0.00000 -1.00000
+DEAL::
+DEAL::0.00000 -1.00000
+DEAL::0.00000 1.00000
+DEAL::1.00000 0.00000
+DEAL::-1.00000 0.00000
+DEAL::
+DEAL::0.00000 1.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::0.00000 0.00000 1.00000
+DEAL::0.00000 0.00000 1.00000
+DEAL::0.00000 1.00000 0.00000
+DEAL::-0.759836 0.759836 0.00000
+DEAL::-0.759836 0.00000 0.759836
+DEAL::
+DEAL::0.00000 1.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::0.707107 0.00000 0.707107
+DEAL::0.00000 1.00000 0.00000
+DEAL::0.707107 0.00000 -0.707107
+DEAL::0.00000 1.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::0.00000 0.707107 0.707107
+DEAL::1.00000 0.00000 0.00000
+DEAL::0.00000 0.707107 -0.707107
+DEAL::
+DEAL::0.00000 1.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::0.00000 0.00000 1.00000
+DEAL::-0.707107 0.707107 0.00000
+DEAL::0.00000 0.00000 1.00000
+DEAL::0.00000 0.00000 1.00000
+DEAL::0.00000 1.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::0.00000 0.00000 1.00000
+DEAL::
+DEAL::0.00000 -1.00000 0.00000
+DEAL::0.00000 0.00000 1.00000
+DEAL::0.00000 1.00000 0.00000
+DEAL::0.00000 0.00000 1.00000
+DEAL::0.00000 0.00000 -1.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::0.00000 0.00000 1.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::-1.00000 0.00000 0.00000
+DEAL::0.00000 1.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::0.00000 1.00000 0.00000
+DEAL::

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.