From: Peter Munch <peterrmuench@gmail.com>
Date: Sun, 17 Jan 2021 10:22:17 +0000 (+0100)
Subject: Generalize VectorTools::get_position_vector for simplex meshes
X-Git-Tag: v9.3.0-rc1~597^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=589854d4f5bffc952fd132d36931b919b54c8611;p=dealii.git

Generalize VectorTools::get_position_vector for simplex meshes
---

diff --git a/include/deal.II/grid/reference_cell.h b/include/deal.II/grid/reference_cell.h
index ca56d20f89..fe79ef35ef 100644
--- a/include/deal.II/grid/reference_cell.h
+++ b/include/deal.II/grid/reference_cell.h
@@ -507,6 +507,18 @@ namespace ReferenceCell
   make_triangulation(const Type &                  reference_cell,
                      Triangulation<dim, spacedim> &tria);
 
+  /**
+   * Return a default mapping of degree @ degree matching the given reference
+   * cell. If this reference cell is a hypercube, then the returned mapping is a
+   * MappingQGeneric; otherwise, it is an object of type MappingFE initialized
+   * with Simplex::FE_P (if the reference cell is a triangle and tetrahedron),
+   * with Simplex::FE_PyramidP (if the reference cell is a pyramid), or with
+   * Simplex::FE_WedgeP (if the reference cell is a wedge).
+   */
+  template <int dim, int spacedim>
+  std::unique_ptr<Mapping<dim, spacedim>>
+  get_default_mapping(const Type &reference_cell, const unsigned int degree);
+
   /**
    * Return a default linear mapping matching the given reference cell.
    * If this reference cell is a hypercube, then the returned mapping
diff --git a/include/deal.II/numerics/vector_tools_interpolate.templates.h b/include/deal.II/numerics/vector_tools_interpolate.templates.h
index 7edd39de41..80138b7c19 100644
--- a/include/deal.II/numerics/vector_tools_interpolate.templates.h
+++ b/include/deal.II/numerics/vector_tools_interpolate.templates.h
@@ -653,8 +653,12 @@ namespace VectorTools
       {
         const Quadrature<dim> quad(fe.get_unit_support_points());
 
-        MappingQGeneric<dim, spacedim> map_q(fe.degree);
-        FEValues<dim, spacedim> fe_v(map_q, fe, quad, update_quadrature_points);
+        const auto map_q = ReferenceCell::get_default_mapping<dim, spacedim>(
+          fe.reference_cell_type(), fe.degree);
+        FEValues<dim, spacedim>              fe_v(*map_q,
+                                     fe,
+                                     quad,
+                                     update_quadrature_points);
         std::vector<types::global_dof_index> dofs(fe.n_dofs_per_cell());
 
         AssertDimension(fe.n_dofs_per_cell(),
diff --git a/source/grid/reference_cell.cc b/source/grid/reference_cell.cc
index f670dde1f8..68a599be18 100644
--- a/source/grid/reference_cell.cc
+++ b/source/grid/reference_cell.cc
@@ -102,6 +102,32 @@ namespace ReferenceCell
   }
 
 
+  template <int dim, int spacedim>
+  std::unique_ptr<Mapping<dim, spacedim>>
+  get_default_mapping(const Type &reference_cell, const unsigned int degree)
+  {
+    AssertDimension(dim, get_dimension(reference_cell));
+
+    if (reference_cell == get_hypercube(dim))
+      return std::make_unique<MappingQGeneric<dim, spacedim>>(degree);
+    else if (reference_cell == Type::Tri || reference_cell == Type::Tet)
+      return std::make_unique<MappingFE<dim, spacedim>>(
+        Simplex::FE_P<dim, spacedim>(degree));
+    else if (reference_cell == Type::Pyramid)
+      return std::make_unique<MappingFE<dim, spacedim>>(
+        Simplex::FE_PyramidP<dim, spacedim>(degree));
+    else if (reference_cell == Type::Wedge)
+      return std::make_unique<MappingFE<dim, spacedim>>(
+        Simplex::FE_WedgeP<dim, spacedim>(degree));
+    else
+      {
+        Assert(false, ExcNotImplemented());
+      }
+
+    return std::make_unique<MappingQGeneric<dim, spacedim>>(degree);
+  }
+
+
   template <int dim, int spacedim>
   const Mapping<dim, spacedim> &
   get_default_linear_mapping(const Type &reference_cell)
diff --git a/source/grid/reference_cell.inst.in b/source/grid/reference_cell.inst.in
index 37f9309819..29eec7c46d 100644
--- a/source/grid/reference_cell.inst.in
+++ b/source/grid/reference_cell.inst.in
@@ -20,6 +20,10 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
       const Type &                                               reference_cell,
       Triangulation<deal_II_dimension, deal_II_space_dimension> &tria);
 
+    template std::unique_ptr<
+      Mapping<deal_II_dimension, deal_II_space_dimension>>
+    get_default_mapping(const Type &reference_cell, const unsigned int degree);
+
     template const Mapping<deal_II_dimension, deal_II_space_dimension>
       &get_default_linear_mapping(const Type &reference_cell);
 
diff --git a/tests/simplex/mapping_fe_fields_01.cc b/tests/simplex/mapping_fe_fields_01.cc
index 4f96c4a0fc..9168411635 100644
--- a/tests/simplex/mapping_fe_fields_01.cc
+++ b/tests/simplex/mapping_fe_fields_01.cc
@@ -14,7 +14,8 @@
 // ---------------------------------------------------------------------
 
 
-// Test MappingFEField for simplex meshes.
+// Test MappingFEField and VectorTools::get_position_vector() for simplex
+// meshes.
 
 #include <deal.II/dofs/dof_handler.h>
 
@@ -36,21 +37,6 @@
 
 using namespace dealii;
 
-template <int dim>
-class Solution : public Function<dim>
-{
-public:
-  Solution()
-    : Function<dim>(dim)
-  {}
-
-  double
-  value(const Point<dim> &point, const unsigned int compontent) const
-  {
-    return point[compontent];
-  }
-};
-
 void
 test()
 {
@@ -59,7 +45,7 @@ test()
   Triangulation<dim> tria;
   GridGenerator::subdivided_hyper_cube_with_simplices(tria, 1);
 
-  Simplex::FE_P<dim> fe(2);
+  Simplex::FE_P<dim> fe(1);
   FESystem<dim>      euler_fe(fe, dim);
 
   DoFHandler<dim> dof_handler(tria);
@@ -70,14 +56,7 @@ test()
 
   Vector<double> euler_vector(euler_dof_handler.n_dofs());
 
-  // TODO: not working (missing mapping)
-  // VectorTools::get_position_vector(euler_dof_handler, euler_vector);
-
-  MappingFE<dim> mapping_interpolation(Simplex::FE_P<dim>(1));
-  VectorTools::interpolate(mapping_interpolation,
-                           euler_dof_handler,
-                           Solution<dim>(),
-                           euler_vector);
+  VectorTools::get_position_vector(euler_dof_handler, euler_vector);
 
   MappingFEField<dim> mapping(euler_dof_handler, euler_vector);