]> https://gitweb.dealii.org/ - dealii.git/commitdiff
MatrixFree: set up shape_values_face and shape_gradients_face 11261/head
authorPeter Munch <peterrmuench@gmail.com>
Fri, 27 Nov 2020 08:08:58 +0000 (09:08 +0100)
committerPeter Munch <peterrmuench@gmail.com>
Sat, 28 Nov 2020 07:43:35 +0000 (08:43 +0100)
include/deal.II/matrix_free/mapping_info.templates.h
include/deal.II/matrix_free/matrix_free.templates.h
include/deal.II/matrix_free/shape_info.h
include/deal.II/matrix_free/shape_info.templates.h
include/deal.II/matrix_free/util.h [new file with mode: 0644]

index 984795d3c240144a7bf3c6d2e24d37f2eceff1c3..444d5c1b8c3a517fbf66241bfc9b8ca7294c3d26 100644 (file)
@@ -30,6 +30,7 @@
 
 #include <deal.II/matrix_free/evaluation_template_factory.h>
 #include <deal.II/matrix_free/mapping_info.h>
+#include <deal.II/matrix_free/util.h>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -71,7 +72,7 @@ namespace internal
       for (unsigned int i = 0; i < n_q_points; ++i)
         quadrature_weights[i] = quadrature.weight(i);
 
-      // note: quadrature_1 and tensor_quadrature_weights are not set up
+      // note: quadrature_1d and tensor_quadrature_weights are not set up
 
       // TODO: set up face_orientations
       (void)update_flags_inner_faces;
@@ -404,15 +405,17 @@ namespace internal
               if (flag == false)
                 {
 #ifdef DEAL_II_WITH_SIMPLEX_SUPPORT
-                  continue; // TODO: face data is not set up
+                  const auto quad_face = get_face_quadrature(quad[my_q][hpq]);
+                  face_data[my_q].descriptor[hpq].initialize(quad_face,
+                                                             update_default);
 #else
                   Assert(false, ExcNotImplemented());
 #endif
                 }
-
-              face_data[my_q].descriptor[hpq].initialize(
-                quad[my_q][hpq].get_tensor_basis()[0],
-                update_flags_boundary_faces);
+              else
+                face_data[my_q].descriptor[hpq].initialize(
+                  quad[my_q][hpq].get_tensor_basis()[0],
+                  update_flags_boundary_faces);
             }
 
           face_data_by_cells[my_q].descriptor.resize(n_hp_quads);
@@ -428,14 +431,16 @@ namespace internal
               if (flag == false)
                 {
 #ifdef DEAL_II_WITH_SIMPLEX_SUPPORT
-                  continue; // TODO: face data is not set up
+                  const auto quad_face = get_face_quadrature(quad[my_q][hpq]);
+                  face_data_by_cells[my_q].descriptor[hpq].initialize(
+                    quad_face, update_default);
 #else
                   Assert(false, ExcNotImplemented());
 #endif
                 }
-
-              face_data_by_cells[my_q].descriptor[hpq].initialize(
-                quad[my_q][hpq].get_tensor_basis()[0], update_default);
+              else
+                face_data_by_cells[my_q].descriptor[hpq].initialize(
+                  quad[my_q][hpq].get_tensor_basis()[0], update_default);
             }
         }
 
index f2e2325cf1941ad0603850804fac4f16032cb8b4..e2c0e006f36f28de026ac90090503e6e740ea239 100644 (file)
@@ -41,6 +41,8 @@
 #include <deal.II/matrix_free/face_setup_internal.h>
 #include <deal.II/matrix_free/matrix_free.h>
 
+#include <deal.II/simplex/quadrature_lib.h>
+
 #ifdef DEAL_II_WITH_TBB
 #  include <deal.II/base/parallel.h>
 
@@ -1335,7 +1337,8 @@ MatrixFree<dim, Number, VectorizedArrayType>::initialize_indices(
   Table<2, internal::MatrixFreeFunctions::ShapeInfo<double>> shape_info_dummy(
     shape_info.size(0), shape_info.size(2));
   {
-    QGauss<dim> quad(1);
+    Quadrature<dim> quad(QGauss<dim>(1));
+    Quadrature<dim> quad_simplex(Simplex::QGauss<dim>(1));
     for (unsigned int no = 0, c = 0; no < dof_handlers.size(); no++)
       for (unsigned int b = 0;
            b < dof_handlers[no]->get_fe(0).n_base_elements();
@@ -1343,9 +1346,13 @@ MatrixFree<dim, Number, VectorizedArrayType>::initialize_indices(
         for (unsigned int fe_no = 0;
              fe_no < dof_handlers[no]->get_fe_collection().size();
              ++fe_no)
-          shape_info_dummy(c, fe_no).reinit(quad,
-                                            dof_handlers[no]->get_fe(fe_no),
-                                            b);
+          shape_info_dummy(c, fe_no).reinit(
+            dof_handlers[no]->get_fe(fe_no).reference_cell_type() ==
+                ReferenceCell::get_hypercube(dim) ?
+              quad :
+              quad_simplex,
+            dof_handlers[no]->get_fe(fe_no),
+            b);
   }
 
   const unsigned int n_lanes     = VectorizedArrayType::size();
index 8b7eba41b8128bca49d961ab30a63bbe4df8dedc..23d07de2bbd639b3ba95be07f55e6ae999ce18fa 100644 (file)
@@ -282,6 +282,20 @@ namespace internal
        * end points of the unit cell.
        */
       bool nodal_at_cell_boundaries;
+
+      /**
+       * Stores the shape values of the finite element evaluated on all
+       * quadrature points for all faces and orientations (no tensor-product
+       * structure exploited).
+       */
+      Table<3, Number> shape_values_face;
+
+      /**
+       * Stores the shape gradients of the finite element evaluated on all
+       * quadrature points for all faces, orientations, and directions
+       * (no tensor-product structure  exploited).
+       */
+      Table<4, Number> shape_gradients_face;
     };
 
 
index 81280e4f480239b4dd9c6019e9080d31d182f49d..90f7cea9c39aa04e170d791686d472a020b75f3b 100644 (file)
@@ -22,6 +22,7 @@
 #include <deal.II/base/memory_consumption.h>
 #include <deal.II/base/polynomial.h>
 #include <deal.II/base/polynomials_piecewise.h>
+#include <deal.II/base/qprojector.h>
 #include <deal.II/base/tensor_product_polynomials.h>
 #include <deal.II/base/utilities.h>
 
 #include <deal.II/fe/fe_poly.h>
 #include <deal.II/fe/fe_q_dg0.h>
 
+#include <deal.II/grid/reference_cell.h>
+
 #include <deal.II/lac/householder.h>
 
 #include <deal.II/matrix_free/shape_info.h>
+#include <deal.II/matrix_free/util.h>
 
 #include <deal.II/simplex/fe_lib.h>
 
@@ -184,8 +188,11 @@ namespace internal
             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;
+          auto &shape_values      = univariate_shape_data.shape_values;
+          auto &shape_values_face = univariate_shape_data.shape_values_face;
+          auto &shape_gradients   = univariate_shape_data.shape_gradients;
+          auto &shape_gradients_face =
+            univariate_shape_data.shape_gradients_face;
 
           const unsigned int n_dofs = fe->n_dofs_per_cell();
 
@@ -207,6 +214,64 @@ namespace internal
                                   q] = grad[d];
               }
 
+          const auto reference_cell_type = ReferenceCell::get_simplex(dim);
+
+          try
+            {
+              const auto quad_face  = get_face_quadrature(quad);
+              this->n_q_points_face = quad_face.size();
+
+              const unsigned int n_face_orientations = dim == 2 ? 2 : 6;
+              const unsigned int n_faces =
+                ReferenceCell::internal::Info::get_cell(reference_cell_type)
+                  .n_faces();
+
+              const auto projected_quad_face =
+                QProjector<dim>::project_to_all_faces(reference_cell_type,
+                                                      quad_face);
+
+              shape_values_face.reinit(
+                {n_faces, n_face_orientations, n_dofs * n_q_points_face});
+
+              shape_gradients_face.reinit(
+                {n_faces, n_face_orientations, dim, n_dofs * n_q_points_face});
+
+              for (unsigned int f = 0; f < n_faces; ++f)
+                for (unsigned int o = 0; o < n_face_orientations; ++o)
+                  {
+                    const auto offset =
+                      QProjector<dim>::DataSetDescriptor::face(
+                        reference_cell_type,
+                        f,
+                        (o ^ 1) & 1,  // face_orientation
+                        (o >> 1) & 1, // face_flip
+                        (o >> 2) & 1, // face_rotation
+                        n_q_points_face);
+
+                    for (unsigned int i = 0; i < n_dofs; ++i)
+                      for (unsigned int q = 0; q < n_q_points_face; ++q)
+                        {
+                          const auto point =
+                            projected_quad_face.point(q + offset);
+
+                          shape_values_face(f, o, i * n_q_points_face + q) =
+                            fe->shape_value(i, point);
+
+                          const auto grad = fe->shape_grad(i, point);
+
+                          for (int d = 0; d < dim; ++d)
+                            shape_gradients_face(
+                              f, o, d, i * n_q_points_face + q) = grad[d];
+                        }
+                  }
+            }
+          catch (...)
+            {
+              // TODO: this might happen if the quadrature rule and the
+              // the FE do not match
+              this->n_q_points_face = 0;
+            }
+
           // TODO: also fill shape_hessians, inverse_shape_values,
           //   shape_data_on_face, quadrature_data_on_face,
           //   values_within_subface, gradients_within_subface,
diff --git a/include/deal.II/matrix_free/util.h b/include/deal.II/matrix_free/util.h
new file mode 100644 (file)
index 0000000..b46f3c2
--- /dev/null
@@ -0,0 +1,53 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+#ifndef dealii_matrix_free_util_h
+#define dealii_matrix_free_util_h
+
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/quadrature.h>
+
+#include <deal.II/simplex/quadrature_lib.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace internal
+{
+  namespace MatrixFreeFunctions
+  {
+    template <int dim>
+    inline Quadrature<dim - 1>
+    get_face_quadrature(const Quadrature<dim> &quad)
+    {
+      if (dim == 2 || dim == 3)
+        for (unsigned int i = 1; i <= 3; ++i)
+          if (quad == Simplex::QGauss<dim>(i))
+            return Simplex::QGauss<dim - 1>(i);
+
+      AssertThrow(false, ExcNotImplemented());
+
+      return Quadrature<dim - 1>();
+    }
+
+  } // end of namespace MatrixFreeFunctions
+} // end of namespace internal
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif

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.