]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Make ShapeInfo compatible with Simplex::FE_P and Simplex::QGauss 11186/head
authorPeter Munch <peterrmuench@gmail.com>
Sun, 15 Nov 2020 09:12:40 +0000 (10:12 +0100)
committerPeter Munch <peterrmuench@gmail.com>
Tue, 17 Nov 2020 13:37:36 +0000 (14:37 +0100)
include/deal.II/matrix_free/shape_info.h
include/deal.II/matrix_free/shape_info.templates.h
tests/simplex/matrix_free_shape_info_01.cc [new file with mode: 0644]
tests/simplex/matrix_free_shape_info_01.with_simplex_support=on.output [new file with mode: 0644]

index bb896a9c54b24bbd094c644e5de27274021759b7..12f3b6eea1a13f3c9de6ac70fe77882f6ed950b8 100644 (file)
@@ -91,7 +91,12 @@ namespace internal
        * of the unit interval 0.5 that additionally add a constant shape
        * function according to FE_Q_DG0.
        */
-      tensor_symmetric_plus_dg0 = 5
+      tensor_symmetric_plus_dg0 = 5,
+
+      /**
+       * Shape functions without an tensor product properties.
+       */
+      tensor_none = 6
     };
 
 
index 873300c9e5e0cbfc48047e1da0c57908d93357fa..802607bfd6022c4566a57bf39a530ae5586d92ac 100644 (file)
@@ -34,6 +34,8 @@
 
 #include <deal.II/matrix_free/shape_info.h>
 
+#include <deal.II/simplex/fe_lib.h>
+
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -134,7 +136,99 @@ namespace internal
                               const FiniteElement<dim> &fe_in,
                               const unsigned int        base_element_number)
     {
+#ifdef DEAL_II_WITH_SIMPLEX_SUPPORT
+      if (quad_in.is_tensor_product() == false ||
+          dynamic_cast<const Simplex::FE_P<dim> *>(&fe_in) ||
+          dynamic_cast<const Simplex::FE_DGP<dim> *>(&fe_in))
+        {
+          // specialization for arbitrary finite elements and quadrature rules
+          // as needed in the context, e.g., of simplices
+
+          AssertDimension(dim, dim_q);
+
+          const auto quad            = Quadrature<dim>(quad_in);
+          const auto fe              = &fe_in.base_element(base_element_number);
+          n_dimensions               = dim;
+          n_components               = fe_in.n_components();
+          n_q_points                 = quad.size();
+          dofs_per_component_on_cell = fe->n_dofs_per_cell();
+          n_q_points_face            = 0; // not implemented yet
+          dofs_per_component_on_face = 0; //
+
+          Assert(fe->n_components() == 1,
+                 ExcMessage(
+                   "FEEvaluation only works for scalar finite elements."));
+
+          data.resize(1);
+          UnivariateShapeData<Number> &univariate_shape_data = data.front();
+          data_access.reinit(n_dimensions, n_components);
+          data_access.fill(&univariate_shape_data);
+
+          // note: we cannot write `univariate_shape_data.quadrature = quad`,
+          // since the quadrature rule within UnivariateShapeData expects
+          // a 1D quadrature rule. However, in this case we are not able to
+          // define that rule anyway so other code cannot use this information.
+
+          univariate_shape_data.fe_degree     = fe->degree;
+          univariate_shape_data.n_q_points_1d = quad.size();
+
+          if ((fe->n_dofs_per_cell() == 0) || (quad.size() == 0))
+            return;
+
+          // grant write access to common univariate shape data
+          auto &shape_values    = univariate_shape_data.shape_values;
+          auto &shape_gradients = univariate_shape_data.shape_gradients;
+
+          const unsigned int n_dofs = fe->n_dofs_per_cell();
+
+          const unsigned int array_size = n_dofs * n_q_points;
+
+          shape_values.resize_fast(array_size);
+          shape_gradients.resize_fast(array_size * dim);
+
+          for (unsigned int i = 0; i < n_dofs; ++i)
+            for (unsigned int q = 0; q < n_q_points; ++q)
+              {
+                shape_values[i * n_q_points + q] =
+                  fe->shape_value(i, quad.point(q));
+
+                const auto grad = fe->shape_grad(i, quad.point(q));
+
+                for (int d = 0; d < dim; ++d)
+                  shape_gradients[d * n_dofs * n_q_points + i * n_q_points +
+                                  q] = grad[d];
+              }
+
+          // TODO: also fill shape_hessians, inverse_shape_values,
+          //   shape_data_on_face, quadrature_data_on_face,
+          //   values_within_subface, gradients_within_subface,
+          //   hessians_within_subface
+
+          // note: shape_gradients_collocation, shape_hessians_collocation,
+          //  shape_values_eo, shape_gradients_eo, shape_hessians_eo,
+          //  shape_gradients_collocation_eo, shape_hessians_collocation_eo,
+          //  inverse_shape_values_eo cannot be filled
+
+          // indicate that no tensor product properties could be exploited
+          element_type                       = tensor_none;
+          univariate_shape_data.element_type = this->element_type;
+
+          univariate_shape_data.nodal_at_cell_boundaries = true;
+
+          lexicographic_numbering.resize(n_dofs);
+
+          for (unsigned int i = 0; i < n_dofs; ++i)
+            lexicographic_numbering[i] = i;
+
+          // TODO: setup face_to_cell_index_nodal, face_to_cell_index_hermite,
+          //  face_orientations
+
+          return;
+        }
+#else
       Assert(quad_in.is_tensor_product(), ExcNotImplemented());
+#endif
+
       const auto quad = quad_in.get_tensor_basis()[0];
 
       const FiniteElement<dim> *fe = &fe_in.base_element(base_element_number);
diff --git a/tests/simplex/matrix_free_shape_info_01.cc b/tests/simplex/matrix_free_shape_info_01.cc
new file mode 100644 (file)
index 0000000..cfe1180
--- /dev/null
@@ -0,0 +1,194 @@
+// ---------------------------------------------------------------------
+//
+// 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 ShapeData for Simplex::FE_P and Simplex::QGauss
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/matrix_free/matrix_free.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 ExactSolution : public Function<dim>
+{
+public:
+  ExactSolution()
+  {}
+
+  virtual double
+  value(const Point<dim> &p, const unsigned int /*component*/ = 0) const
+  {
+    return p[0] + 2 * p[1];
+  }
+};
+
+template <int dim, typename Number>
+class FEEvaluationDummy
+{
+public:
+  FEEvaluationDummy(
+    const internal::MatrixFreeFunctions::ShapeInfo<Number> &shape_info)
+    : shape_info(shape_info)
+  {
+    values.resize(shape_info.n_q_points);
+    gradients.resize(shape_info.n_q_points * dim);
+  }
+
+  void
+  evaluate(const Vector<Number> &dof_values)
+  {
+    const unsigned int n_dofs     = shape_info.dofs_per_component_on_cell;
+    const unsigned int n_q_points = shape_info.n_q_points;
+
+    for (unsigned int q = 0, c = 0; q < shape_info.n_q_points; ++q)
+      {
+        values[q] = 0.0;
+        for (unsigned int i = 0; i < n_dofs; ++i, ++c)
+          values[q] += shape_info.data[0].shape_values[c] * dof_values[i];
+      }
+
+    for (unsigned int d = 0; d < dim; ++d)
+      for (unsigned int q = 0; q < n_q_points; ++q)
+        {
+          gradients[n_q_points * d + q] = 0.0;
+          for (unsigned int i = 0; i < n_dofs; ++i)
+            gradients[n_q_points * d + q] +=
+              shape_info.data[0]
+                .shape_gradients[d * n_dofs * n_q_points + i * n_q_points + q] *
+              dof_values[i];
+        }
+  }
+
+  Number
+  get_value(const unsigned int q_point) const
+  {
+    return values[q_point];
+  }
+
+  Tensor<1, dim, Number>
+  get_gradient(const unsigned int q_point) const
+  {
+    Tensor<1, dim, Number> out;
+
+    for (unsigned int c = 0; c < dim; ++c)
+      out[c] = gradients[shape_info.n_q_points * c + q_point];
+
+    return out;
+  }
+
+private:
+  const internal::MatrixFreeFunctions::ShapeInfo<Number> &shape_info;
+
+  AlignedVector<Number> values;
+  AlignedVector<Number> gradients;
+};
+
+template <int dim, int spacedim = dim, typename Number = double>
+void
+test(const FiniteElement<dim, spacedim> &fe)
+{
+  Triangulation<dim, spacedim> tria;
+  GridGenerator::subdivided_hyper_cube_with_simplices(tria, 1);
+
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+
+  MappingFE<dim> mapping(Simplex::FE_P<dim>(1));
+
+  Simplex::QGauss<dim> quadrature(1);
+
+  internal::MatrixFreeFunctions::ShapeInfo<Number> shape_info(quadrature, fe);
+
+  FEValues<dim> fe_values(mapping,
+                          fe,
+                          quadrature,
+                          update_values | update_gradients |
+                            update_inverse_jacobians);
+
+  Vector<Number> src(dof_handler.n_dofs());
+  VectorTools::interpolate(mapping, dof_handler, ExactSolution<dim>(), src);
+
+  FEEvaluationDummy<dim, Number> fe_eval(shape_info);
+
+  for (const auto &cell : dof_handler.active_cell_iterators())
+    {
+      fe_values.reinit(cell);
+
+      Vector<Number> dof_values(fe_values.dofs_per_cell);
+      cell->get_dof_values(src, dof_values);
+      fe_eval.evaluate(dof_values);
+
+      // check values
+      {
+        std::vector<Number> values(fe_values.n_quadrature_points);
+        fe_values.get_function_values(src, values);
+
+        for (const auto q : fe_values.quadrature_point_indices())
+          {
+            Assert(std::abs(values[q] - fe_eval.get_value(q)) < 1e-10,
+                   ExcMessage("Entries do not match!"));
+          }
+      }
+
+      // check gradients
+      {
+        std::vector<Tensor<1, dim, Number>> gradients(
+          fe_values.n_quadrature_points);
+        fe_values.get_function_gradients(src, gradients);
+
+        for (const auto q : fe_values.quadrature_point_indices())
+          {
+            // apply the transposed inverse Jacobian of the mapping
+            Tensor<1, dim, Number> temp = fe_eval.get_gradient(q);
+            Tensor<1, dim, Number> temp_transformed;
+            for (int d = 0; d < dim; ++d)
+              {
+                Number sum = 0;
+                for (int e = 0; e < dim; ++e)
+                  sum += fe_values.inverse_jacobian(q)[e][d] * temp[e];
+                temp_transformed[d] = sum;
+              }
+
+            for (int i = 0; i < dim; ++i)
+              Assert(std::abs(gradients[q][i] - temp_transformed[i]) < 1e-10,
+                     ExcMessage("Entries do not match!"));
+          }
+      }
+    }
+}
+
+int
+main()
+{
+  initlog();
+
+  test<2>(Simplex::FE_P<2>(2));
+}
diff --git a/tests/simplex/matrix_free_shape_info_01.with_simplex_support=on.output b/tests/simplex/matrix_free_shape_info_01.with_simplex_support=on.output
new file mode 100644 (file)
index 0000000..8b13789
--- /dev/null
@@ -0,0 +1 @@
+

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.