]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Work on Work on MappingFEField for simplex meshes 11565/head
authorPeter Munch <peterrmuench@gmail.com>
Fri, 15 Jan 2021 19:47:58 +0000 (20:47 +0100)
committerPeter Munch <peterrmuench@gmail.com>
Sat, 16 Jan 2021 20:11:12 +0000 (21:11 +0100)
include/deal.II/fe/mapping_fe_field.h
include/deal.II/grid/reference_cell.h
source/fe/mapping_fe_field.cc
source/grid/reference_cell.cc
source/grid/reference_cell.inst.in
tests/simplex/mapping_fe_fields_01.cc [new file with mode: 0644]
tests/simplex/mapping_fe_fields_01.with_simplex_support=on.output [new file with mode: 0644]

index 494f5bc8ac2db16e8ff0d4a3c4b64a2abc1eb2cf..7f36f19cf96b992adea24920c33ef27e06e189ed 100644 (file)
@@ -13,8 +13,8 @@
 //
 // ---------------------------------------------------------------------
 
-#ifndef dealii_mapping_fe_h
-#define dealii_mapping_fe_h
+#ifndef dealii_mapping_fe_field_h
+#define dealii_mapping_fe_field_h
 
 
 #include <deal.II/base/config.h>
index 019d41062a1fa659d8157390f7842fbcf23efe1f..2d8161feef46aee3f7aaaa8bce4aa7cc87d0ed89 100644 (file)
@@ -506,6 +506,16 @@ namespace ReferenceCell
   get_gauss_type_quadrature(const Type &   reference_cell,
                             const unsigned n_points_1D);
 
+  /**
+   * Return a quadrature rule with the support points of the given reference
+   * cell.
+   *
+   * @note The weights are not filled.
+   */
+  template <int dim>
+  Quadrature<dim> &
+  get_nodal_type_quadrature(const Type &reference_cell);
+
   namespace internal
   {
     /**
index 1d75c345e107c02239f4ec8a9328d4abd921529c..d440769a9a95f559a1d925d38fb059ab816cc58e 100644 (file)
@@ -192,28 +192,6 @@ MappingFEField<dim, spacedim, VectorType, void>::InternalData::
 
 
 
-namespace
-{
-  // Construct a quadrature formula containing the vertices of the reference
-  // cell in dimension dim (with invalid weights)
-  template <int dim>
-  Quadrature<dim> &
-  get_vertex_quadrature()
-  {
-    static Quadrature<dim> quad;
-    if (quad.size() == 0)
-      {
-        std::vector<Point<dim>> points(GeometryInfo<dim>::vertices_per_cell);
-        for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
-          points[i] = GeometryInfo<dim>::unit_cell_vertex(i);
-        quad = Quadrature<dim>(points);
-      }
-    return quad;
-  }
-} // namespace
-
-
-
 template <int dim, int spacedim, typename VectorType>
 MappingFEField<dim, spacedim, VectorType, void>::MappingFEField(
   const DoFHandler<dim, spacedim> &euler_dof_handler,
@@ -229,7 +207,8 @@ MappingFEField<dim, spacedim, VectorType, void>::MappingFEField(
                 true))
   , fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int)
   , fe_values(this->euler_dof_handler->get_fe(),
-              get_vertex_quadrature<dim>(),
+              ReferenceCell::get_nodal_type_quadrature<dim>(
+                this->euler_dof_handler->get_fe().reference_cell_type()),
               update_values)
 {
   unsigned int size = 0;
@@ -257,7 +236,8 @@ MappingFEField<dim, spacedim, VectorType, void>::MappingFEField(
                 true))
   , fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int)
   , fe_values(this->euler_dof_handler->get_fe(),
-              get_vertex_quadrature<dim>(),
+              ReferenceCell::get_nodal_type_quadrature<dim>(
+                this->euler_dof_handler->get_fe().reference_cell_type()),
               update_values)
 {
   unsigned int size = 0;
@@ -296,7 +276,8 @@ MappingFEField<dim, spacedim, VectorType, void>::MappingFEField(
                 true))
   , fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int)
   , fe_values(this->euler_dof_handler->get_fe(),
-              get_vertex_quadrature<dim>(),
+              ReferenceCell::get_nodal_type_quadrature<dim>(
+                this->euler_dof_handler->get_fe().reference_cell_type()),
               update_values)
 {
   unsigned int size = 0;
@@ -332,7 +313,8 @@ MappingFEField<dim, spacedim, VectorType, void>::MappingFEField(
   , fe_mask(mapping.fe_mask)
   , fe_to_real(mapping.fe_to_real)
   , fe_values(mapping.euler_dof_handler->get_fe(),
-              get_vertex_quadrature<dim>(),
+              ReferenceCell::get_nodal_type_quadrature<dim>(
+                this->euler_dof_handler->get_fe().reference_cell_type()),
               update_values)
 {}
 
@@ -372,8 +354,7 @@ MappingFEField<dim, spacedim, VectorType, void>::get_vertices(
     *cell, euler_dof_handler);
 
   Assert(uses_level_dofs || dof_cell->is_active() == true, ExcInactiveCell());
-  AssertDimension(GeometryInfo<dim>::vertices_per_cell,
-                  fe_values.n_quadrature_points);
+  AssertDimension(cell->n_vertices(), fe_values.n_quadrature_points);
   AssertDimension(fe_to_real.size(),
                   euler_dof_handler->get_fe().n_components());
   if (uses_level_dofs)
@@ -402,7 +383,7 @@ MappingFEField<dim, spacedim, VectorType, void>::get_vertices(
 
   boost::container::small_vector<Point<spacedim>,
                                  GeometryInfo<dim>::vertices_per_cell>
-    vertices(GeometryInfo<dim>::vertices_per_cell);
+    vertices(cell->n_vertices());
   for (unsigned int i = 0; i < dofs_per_cell; ++i)
     {
       const unsigned int comp = fe_to_real
@@ -412,7 +393,7 @@ MappingFEField<dim, spacedim, VectorType, void>::get_vertices(
           typename VectorType::value_type value =
             internal::ElementAccess<VectorType>::get(vector, dof_indices[i]);
           if (euler_dof_handler->get_fe().is_primitive(i))
-            for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
+            for (const unsigned int v : cell->vertex_indices())
               vertices[v][comp] += fe_values.shape_value(i, v) * value;
           else
             Assert(false, ExcNotImplemented());
@@ -585,23 +566,30 @@ MappingFEField<dim, spacedim, VectorType, void>::compute_face_data(
           data.aux.resize(
             dim - 1, std::vector<Tensor<1, spacedim>>(n_original_q_points));
 
+
+          // TODO: only a single reference cell type possible...
+          const auto reference_cell_type =
+            this->euler_dof_handler->get_fe().reference_cell_type();
+          const auto n_faces =
+            ReferenceCell::internal::Info::get_cell(reference_cell_type)
+              .n_faces();
+
           // Compute tangentials to the unit cell.
-          for (const unsigned int i : GeometryInfo<dim>::face_indices())
+          for (unsigned int i = 0; i < n_faces; ++i)
             {
               data.unit_tangentials[i].resize(n_original_q_points);
               std::fill(data.unit_tangentials[i].begin(),
                         data.unit_tangentials[i].end(),
-                        GeometryInfo<dim>::unit_tangential_vectors[i][0]);
+                        ReferenceCell::unit_tangential_vectors<dim>(
+                          reference_cell_type, i, 0));
               if (dim > 2)
                 {
-                  data.unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
-                    .resize(n_original_q_points);
-                  std::fill(
-                    data.unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
-                      .begin(),
-                    data.unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
-                      .end(),
-                    GeometryInfo<dim>::unit_tangential_vectors[i][1]);
+                  data.unit_tangentials[n_faces + i].resize(
+                    n_original_q_points);
+                  std::fill(data.unit_tangentials[n_faces + i].begin(),
+                            data.unit_tangentials[n_faces + i].end(),
+                            ReferenceCell::unit_tangential_vectors<dim>(
+                              reference_cell_type, i, 1));
                 }
             }
         }
@@ -1273,22 +1261,17 @@ namespace internal
             // 0.
             for (unsigned int d = 0; d != dim - 1; ++d)
               {
-                Assert(face_no + GeometryInfo<dim>::faces_per_cell * d <
+                Assert(face_no + cell->n_faces() * d <
                          data.unit_tangentials.size(),
                        ExcInternalError());
                 Assert(
                   data.aux[d].size() <=
-                    data
-                      .unit_tangentials[face_no +
-                                        GeometryInfo<dim>::faces_per_cell * d]
-                      .size(),
+                    data.unit_tangentials[face_no + cell->n_faces() * d].size(),
                   ExcInternalError());
 
                 mapping.transform(
                   make_array_view(
-                    data
-                      .unit_tangentials[face_no +
-                                        GeometryInfo<dim>::faces_per_cell * d]),
+                    data.unit_tangentials[face_no + cell->n_faces() * d]),
                   mapping_contravariant,
                   data,
                   make_array_view(data.aux[d]));
@@ -1371,6 +1354,7 @@ namespace internal
 
                       if (subface_no != numbers::invalid_unsigned_int)
                         {
+                          // TODO
                           const double area_ratio =
                             GeometryInfo<dim>::subface_ratio(
                               cell->subface_case(face_no), subface_no);
@@ -2096,6 +2080,7 @@ MappingFEField<dim, spacedim, VectorType, void>::transform_real_to_unit_cell(
         initial_p_unit[d] = 0.5;
     }
 
+  // TODO
   initial_p_unit = GeometryInfo<dim>::project_to_unit_cell(initial_p_unit);
 
   // for (unsigned int d=0; d<dim; ++d)
index 816da89dad9309d7ab5b24b85446feb8bc88eb7e..f670dde1f8fd8fb04c027f2e6b65214f6e10d405 100644 (file)
@@ -181,6 +181,47 @@ namespace ReferenceCell
     return Quadrature<dim>(); // never reached
   }
 
+  template <int dim>
+  Quadrature<dim> &
+  get_nodal_type_quadrature(const Type &reference_cell)
+  {
+    AssertDimension(dim, get_dimension(reference_cell));
+
+    const auto create_quadrature = [](const Type &reference_cell) {
+      Triangulation<dim> tria;
+      make_triangulation(reference_cell, tria);
+
+      return Quadrature<dim>(tria.get_vertices());
+    };
+
+    if (reference_cell == get_hypercube(dim))
+      {
+        static Quadrature<dim> quadrature = create_quadrature(reference_cell);
+        return quadrature;
+      }
+    else if (reference_cell == Type::Tri || reference_cell == Type::Tet)
+      {
+        static Quadrature<dim> quadrature = create_quadrature(reference_cell);
+        return quadrature;
+      }
+    else if (reference_cell == Type::Pyramid)
+      {
+        static Quadrature<dim> quadrature = create_quadrature(reference_cell);
+        return quadrature;
+      }
+    else if (reference_cell == Type::Wedge)
+      {
+        static Quadrature<dim> quadrature = create_quadrature(reference_cell);
+        return quadrature;
+      }
+    else
+      Assert(false, ExcNotImplemented());
+
+    static Quadrature<dim> dummy;
+
+    return dummy; // never reached
+  }
+
 #include "reference_cell.inst"
 
 } // namespace ReferenceCell
index 198079f766afc9e81eb39629a7feedf3cacd9980..37f930981911fb7d04e8b197241493db938daeac 100644 (file)
@@ -34,4 +34,6 @@ for (deal_II_dimension : DIMENSIONS)
   {
     template Quadrature<deal_II_dimension> get_gauss_type_quadrature(
       const Type &reference_cell, const unsigned n_points_1D);
+    template Quadrature<deal_II_dimension> &get_nodal_type_quadrature(
+      const Type &reference_cell);
   }
diff --git a/tests/simplex/mapping_fe_fields_01.cc b/tests/simplex/mapping_fe_fields_01.cc
new file mode 100644 (file)
index 0000000..4f96c4a
--- /dev/null
@@ -0,0 +1,121 @@
+// ---------------------------------------------------------------------
+//
+// 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 MappingFEField for simplex meshes.
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_fe.h>
+#include <deal.II/fe/mapping_fe_field.h>
+
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <deal.II/simplex/fe_lib.h>
+#include <deal.II/simplex/grid_generator.h>
+#include <deal.II/simplex/quadrature_lib.h>
+
+#include "../tests.h"
+
+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()
+{
+  const int dim = 2;
+
+  Triangulation<dim> tria;
+  GridGenerator::subdivided_hyper_cube_with_simplices(tria, 1);
+
+  Simplex::FE_P<dim> fe(2);
+  FESystem<dim>      euler_fe(fe, dim);
+
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+
+  DoFHandler<dim> euler_dof_handler(tria);
+  euler_dof_handler.distribute_dofs(euler_fe);
+
+  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);
+
+  MappingFEField<dim> mapping(euler_dof_handler, euler_vector);
+
+  Simplex::QGauss<dim> quadrature_formula(1);
+
+  FEValues<dim> fe_values(mapping,
+                          fe,
+                          quadrature_formula,
+                          update_values | update_gradients);
+
+  for (const auto &cell : dof_handler.active_cell_iterators())
+    {
+      fe_values.reinit(cell);
+    }
+
+  {
+    DataOut<dim> data_out;
+
+    data_out.attach_dof_handler(dof_handler);
+
+    Vector<double> solution(dof_handler.n_dofs());
+    data_out.add_data_vector(solution, "solution");
+
+    data_out.build_patches(mapping, 2);
+
+#if false
+    std::ofstream output("test.vtk");
+    data_out.write_vtk(output);
+#else
+    data_out.write_vtk(deallog.get_file_stream());
+#endif
+  }
+}
+
+int
+main()
+{
+  initlog();
+
+  test();
+}
\ No newline at end of file
diff --git a/tests/simplex/mapping_fe_fields_01.with_simplex_support=on.output b/tests/simplex/mapping_fe_fields_01.with_simplex_support=on.output
new file mode 100644 (file)
index 0000000..146f92e
--- /dev/null
@@ -0,0 +1,30 @@
+
+# vtk DataFile Version 3.0
+#This file was generated 
+ASCII
+DATASET UNSTRUCTURED_GRID
+
+POINTS 12 double
+0.00000 0.00000 0
+1.00000 0.00000 0
+0.00000 1.00000 0
+0.500000 0.00000 0
+0.500000 0.500000 0
+0.00000 0.500000 0
+1.00000 1.00000 0
+0.00000 1.00000 0
+1.00000 0.00000 0
+0.500000 1.00000 0
+0.500000 0.500000 0
+1.00000 0.500000 0
+
+CELLS 2 14
+       6       0       1       2       3       4       5
+       6       6       7       8       9       10      11
+
+CELL_TYPES 2
+ 22 22
+POINT_DATA 12
+SCALARS solution double 1
+LOOKUP_TABLE default
+0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 

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.