]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Make MatrixFree compile for arbitrary dim quadrature rules 11168/head
authorPeter Munch <peterrmuench@gmail.com>
Thu, 12 Nov 2020 13:01:23 +0000 (14:01 +0100)
committerPeter Munch <peterrmuench@gmail.com>
Fri, 13 Nov 2020 13:23:01 +0000 (14:23 +0100)
include/deal.II/hp/q_collection.h
include/deal.II/matrix_free/fe_evaluation.h
include/deal.II/matrix_free/mapping_info.h
include/deal.II/matrix_free/mapping_info.templates.h
include/deal.II/matrix_free/matrix_free.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
source/matrix_free/matrix_free.inst.in
source/matrix_free/shape_info.inst.in

index 0ead13dc43e4a4457ab2f58329ede8596966ae26..073bac7968435bca933b3dfbd3c576ed824058b8 100644 (file)
@@ -52,13 +52,20 @@ namespace hp
      */
     QCollection() = default;
 
+    /**
+     * Copy constructor.
+     */
+    template <int dim_in>
+    QCollection(const QCollection<dim_in> &other);
+
     /**
      * Conversion constructor. This constructor creates a QCollection from a
      * single quadrature rule. More quadrature formulas can be added with
      * push_back(), if desired, though it would probably be clearer to add all
      * mappings the same way.
      */
-    explicit QCollection(const Quadrature<dim> &quadrature);
+    template <int dim_in>
+    explicit QCollection(const Quadrature<dim_in> &quadrature);
 
     /**
      * Constructor. This constructor creates a QCollection from one or
@@ -93,8 +100,9 @@ namespace hp
      * is later destroyed by this object upon destruction of the entire
      * collection.
      */
+    template <int dim_in>
     void
-    push_back(const Quadrature<dim> &new_quadrature);
+    push_back(const Quadrature<dim_in> &new_quadrature);
 
     /**
      * Return a reference to the quadrature rule specified by the argument.
@@ -150,21 +158,41 @@ namespace hp
 
   /* --------------- inline functions ------------------- */
 
+  template <int dim>
+  template <int dim_in>
+  QCollection<dim>::QCollection(const QCollection<dim_in> &other)
+  {
+    for (unsigned int i = 0; i < other.size(); ++i)
+      push_back(other[i]);
+  }
+
+
+
   template <int dim>
   template <class... QTypes>
   QCollection<dim>::QCollection(const QTypes &... quadrature_objects)
   {
-    static_assert(is_base_of_all<Quadrature<dim>, QTypes...>::value,
-                  "Not all of the input arguments of this function "
-                  "are derived from Quadrature<dim>!");
-
     // loop over all of the given arguments and add the quadrature objects to
     // this collection. Inlining the definition of q_pointers causes internal
     // compiler errors on GCC 7.1.1 so we define it separately:
-    const auto q_pointers = {
-      (static_cast<const Quadrature<dim> *>(&quadrature_objects))...};
-    for (const auto p : q_pointers)
-      push_back(*p);
+    if (is_base_of_all<Quadrature<dim>, QTypes...>::value)
+      {
+        const auto q_pointers = {
+          (reinterpret_cast<const Quadrature<dim> *>(&quadrature_objects))...};
+        for (const auto p : q_pointers)
+          push_back(*p);
+      }
+    else if (is_base_of_all<Quadrature<1>, QTypes...>::value)
+      {
+        const auto q_pointers = {
+          (reinterpret_cast<const Quadrature<1> *>(&quadrature_objects))...};
+        for (const auto p : q_pointers)
+          push_back(*p);
+      }
+    else
+      {
+        Assert(false, ExcNotImplemented());
+      }
   }
 
 
@@ -223,7 +251,8 @@ namespace hp
 
 
   template <int dim>
-  inline QCollection<dim>::QCollection(const Quadrature<dim> &quadrature)
+  template <int dim_in>
+  inline QCollection<dim>::QCollection(const Quadrature<dim_in> &quadrature)
   {
     quadratures.push_back(std::make_shared<const Quadrature<dim>>(quadrature));
   }
@@ -239,8 +268,9 @@ namespace hp
 
 
   template <int dim>
+  template <int dim_in>
   inline void
-  QCollection<dim>::push_back(const Quadrature<dim> &new_quadrature)
+  QCollection<dim>::push_back(const Quadrature<dim_in> &new_quadrature)
   {
     quadratures.push_back(
       std::make_shared<const Quadrature<dim>>(new_quadrature));
index 6446d407ddfcfc7ef2691856216ec100d6ccdc2e..f37f136841d4b23033412c5baf5398cec27a43b3 100644 (file)
@@ -3331,7 +3331,7 @@ inline FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
   ,
   // select the correct base element from the given FE component
   data(new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>(
-    quadrature,
+    Quadrature<dim - is_face>(quadrature),
     fe,
     fe.component_to_base_index(first_selected_component).first))
   , jacobian(nullptr)
index b2ba970de7c856701cea70eb23f9592401192267..e4852ea9eb8661b99cbec03b516cfe2ed24ede0f 100644 (file)
@@ -320,17 +320,17 @@ namespace internal
        * CellIterator::level() and CellIterator::index(), in order to allow
        * for different kinds of iterators, e.g. standard DoFHandler,
        * multigrid, etc.)  on a fixed Triangulation. In addition, a mapping
-       * and several quadrature formulas are given.
+       * and several 1D quadrature formulas are given.
        */
       void
       initialize(
         const dealii::Triangulation<dim> &                        tria,
         const std::vector<std::pair<unsigned int, unsigned int>> &cells,
         const FaceInfo<VectorizedArrayType::size()> &             faces,
-        const std::vector<unsigned int> &              active_fe_index,
-        const Mapping<dim> &                           mapping,
-        const std::vector<dealii::hp::QCollection<1>> &quad,
-        const UpdateFlags                              update_flags_cells,
+        const std::vector<unsigned int> &                active_fe_index,
+        const Mapping<dim> &                             mapping,
+        const std::vector<dealii::hp::QCollection<dim>> &quad,
+        const UpdateFlags                                update_flags_cells,
         const UpdateFlags update_flags_boundary_faces,
         const UpdateFlags update_flags_inner_faces,
         const UpdateFlags update_flags_faces_by_cells);
@@ -506,9 +506,10 @@ namespace internal
        * internal functions to initialize all data as requested by the user.
        */
       static UpdateFlags
-      compute_update_flags(const UpdateFlags update_flags,
-                           const std::vector<dealii::hp::QCollection<1>> &quad =
-                             std::vector<dealii::hp::QCollection<1>>());
+      compute_update_flags(
+        const UpdateFlags                                update_flags,
+        const std::vector<dealii::hp::QCollection<dim>> &quad =
+          std::vector<dealii::hp::QCollection<dim>>());
     };
 
 
index 85d747d18ce67e1992af6073b09fd0ff03e8fae3..629b1d6166e204b55b4905a2709cfc9e65f7baac 100644 (file)
@@ -252,8 +252,8 @@ namespace internal
     template <int dim, typename Number, typename VectorizedArrayType>
     UpdateFlags
     MappingInfo<dim, Number, VectorizedArrayType>::compute_update_flags(
-      const UpdateFlags                              update_flags,
-      const std::vector<dealii::hp::QCollection<1>> &quad)
+      const UpdateFlags                                update_flags,
+      const std::vector<dealii::hp::QCollection<dim>> &quad)
     {
       // this class is build around the evaluation of jacobians, so compute
       // them in any case. The Jacobians will be inverted manually. Since we
@@ -304,7 +304,7 @@ namespace internal
       const FaceInfo<VectorizedArrayType::size()> &             face_info,
       const std::vector<unsigned int> &                         active_fe_index,
       const Mapping<dim> &                                      mapping,
-      const std::vector<dealii::hp::QCollection<1>> &           quad,
+      const std::vector<dealii::hp::QCollection<dim>> &         quad,
       const UpdateFlags update_flags_cells,
       const UpdateFlags update_flags_boundary_faces,
       const UpdateFlags update_flags_inner_faces,
@@ -337,18 +337,49 @@ namespace internal
           AssertIndexRange(0, n_hp_quads);
           cell_data[my_q].descriptor.resize(n_hp_quads);
           for (unsigned int q = 0; q < n_hp_quads; ++q)
-            cell_data[my_q].descriptor[q].initialize(quad[my_q][q],
-                                                     update_default);
+            {
+              Assert(quad[my_q][q].is_tensor_product(), ExcNotImplemented());
+              for (unsigned int i = 1; i < dim; ++i)
+                {
+                  Assert(quad[my_q][q].get_tensor_basis()[0] ==
+                           quad[my_q][q].get_tensor_basis()[i],
+                         ExcNotImplemented());
+                }
+
+              cell_data[my_q].descriptor[q].initialize(
+                quad[my_q][q].get_tensor_basis()[0], update_default);
+            }
 
           face_data[my_q].descriptor.resize(n_hp_quads);
           for (unsigned int hpq = 0; hpq < n_hp_quads; ++hpq)
-            face_data[my_q].descriptor[hpq].initialize(
-              quad[my_q][hpq], update_flags_boundary_faces);
+            {
+              Assert(quad[my_q][hpq].is_tensor_product(), ExcNotImplemented());
+              for (unsigned int i = 1; i < dim; ++i)
+                {
+                  Assert(quad[my_q][hpq].get_tensor_basis()[0] ==
+                           quad[my_q][hpq].get_tensor_basis()[i],
+                         ExcNotImplemented());
+                }
+
+              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);
           for (unsigned int hpq = 0; hpq < n_hp_quads; ++hpq)
-            face_data_by_cells[my_q].descriptor[hpq].initialize(quad[my_q][hpq],
-                                                                update_default);
+            {
+              Assert(quad[my_q][hpq].is_tensor_product(), ExcNotImplemented());
+              for (unsigned int i = 1; i < dim; ++i)
+                {
+                  Assert(quad[my_q][hpq].get_tensor_basis()[0] ==
+                           quad[my_q][hpq].get_tensor_basis()[i],
+                         ExcNotImplemented());
+                }
+
+              face_data_by_cells[my_q].descriptor[hpq].initialize(
+                quad[my_q][hpq].get_tensor_basis()[0], update_default);
+            }
         }
 
       // In case we have no hp adaptivity (active_fe_index is empty), we have
@@ -2588,7 +2619,7 @@ namespace internal
       {
         FE_DGQ<dim> fe_geometry(mapping_degree);
         for (unsigned int my_q = 0; my_q < cell_data.size(); ++my_q)
-          shape_infos[my_q].reinit(cell_data[my_q].descriptor[0].quadrature_1d,
+          shape_infos[my_q].reinit(cell_data[my_q].descriptor[0].quadrature,
                                    fe_geometry);
       }
 
index 7dadca5886381be04e40b4f6baf6056b940ad603..ea4576eb496ce1415e5193ad9fceef9e71e2edf6 100644 (file)
@@ -2025,14 +2025,14 @@ private:
    * This is the actual reinit function that sets up the indices for the
    * DoFHandler case.
    */
-  template <typename number2>
+  template <typename number2, int q_dim>
   void
   internal_reinit(
     const Mapping<dim> &                                   mapping,
     const std::vector<const DoFHandler<dim, dim> *> &      dof_handlers,
     const std::vector<const AffineConstraints<number2> *> &constraint,
     const std::vector<IndexSet> &                          locally_owned_set,
-    const std::vector<hp::QCollection<1>> &                quad,
+    const std::vector<hp::QCollection<q_dim>> &            quad,
     const AdditionalData &                                 additional_data);
 
   /**
@@ -2843,7 +2843,7 @@ MatrixFree<dim, Number, VectorizedArrayType>::reinit(
     internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
       dof_handlers, additional_data.mg_level);
 
-  std::vector<hp::QCollection<1>> quad_hp;
+  std::vector<hp::QCollection<dim>> quad_hp;
   quad_hp.emplace_back(quad);
 
   internal_reinit(StaticMappingQ1<dim>::mapping,
@@ -2877,7 +2877,7 @@ MatrixFree<dim, Number, VectorizedArrayType>::reinit(
     internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
       dof_handlers, additional_data.mg_level);
 
-  std::vector<hp::QCollection<1>> quad_hp;
+  std::vector<hp::QCollection<dim>> quad_hp;
   quad_hp.emplace_back(quad);
 
   internal_reinit(mapping,
@@ -2903,7 +2903,7 @@ MatrixFree<dim, Number, VectorizedArrayType>::reinit(
   std::vector<IndexSet> locally_owned_set =
     internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
       dof_handler, additional_data.mg_level);
-  std::vector<hp::QCollection<1>> quad_hp;
+  std::vector<hp::QCollection<dim>> quad_hp;
   for (unsigned int q = 0; q < quad.size(); ++q)
     quad_hp.emplace_back(quad[q]);
   internal_reinit(StaticMappingQ1<dim>::mapping,
@@ -2949,7 +2949,7 @@ MatrixFree<dim, Number, VectorizedArrayType>::reinit(
   std::vector<IndexSet> locally_owned_set =
     internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
       dof_handler, additional_data.mg_level);
-  std::vector<hp::QCollection<1>> quad_hp;
+  std::vector<hp::QCollection<dim>> quad_hp;
   quad_hp.emplace_back(quad);
   internal_reinit(StaticMappingQ1<dim>::mapping,
                   dof_handler,
@@ -2994,7 +2994,7 @@ MatrixFree<dim, Number, VectorizedArrayType>::reinit(
   std::vector<IndexSet> locally_owned_set =
     internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
       dof_handler, additional_data.mg_level);
-  std::vector<hp::QCollection<1>> quad_hp;
+  std::vector<hp::QCollection<dim>> quad_hp;
   quad_hp.emplace_back(quad);
   internal_reinit(mapping,
                   dof_handler,
@@ -3020,7 +3020,7 @@ MatrixFree<dim, Number, VectorizedArrayType>::reinit(
   std::vector<IndexSet> locally_owned_set =
     internal::MatrixFreeImplementation::extract_locally_owned_index_sets(
       dof_handler, additional_data.mg_level);
-  std::vector<hp::QCollection<1>> quad_hp;
+  std::vector<hp::QCollection<dim>> quad_hp;
   for (unsigned int q = 0; q < quad.size(); ++q)
     quad_hp.emplace_back(quad[q]);
   internal_reinit(mapping,
index 03f224613cbe3a63e18757b7a9bbd92b94e654a0..fdd4febe56aeb804801b02e14a1d315634de9c3e 100644 (file)
@@ -275,14 +275,14 @@ MatrixFree<dim, Number, VectorizedArrayType>::copy_from(
 
 
 template <int dim, typename Number, typename VectorizedArrayType>
-template <typename number2>
+template <typename number2, int q_dim>
 void
 MatrixFree<dim, Number, VectorizedArrayType>::internal_reinit(
   const Mapping<dim> &                                   mapping,
   const std::vector<const DoFHandler<dim, dim> *> &      dof_handler,
   const std::vector<const AffineConstraints<number2> *> &constraint,
   const std::vector<IndexSet> &                          locally_owned_dofs,
-  const std::vector<hp::QCollection<1>> &                quad,
+  const std::vector<hp::QCollection<q_dim>> &            quad,
   const typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
     &additional_data)
 {
@@ -1293,7 +1293,7 @@ MatrixFree<dim, Number, VectorizedArrayType>::initialize_indices(
   Table<2, internal::MatrixFreeFunctions::ShapeInfo<double>> shape_info_dummy(
     shape_info.size(0), shape_info.size(2));
   {
-    QGauss<1> quad(1);
+    QGauss<dim> quad(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();
index 64daf31106fba2ff3de9e477a0b20a3c35de0857..bb896a9c54b24bbd094c644e5de27274021759b7 100644 (file)
@@ -312,8 +312,8 @@ namespace internal
       /**
        * Constructor that initializes the data fields using the reinit method.
        */
-      template <int dim>
-      ShapeInfo(const Quadrature<1> &     quad,
+      template <int dim, int dim_q>
+      ShapeInfo(const Quadrature<dim_q> & quad,
                 const FiniteElement<dim> &fe,
                 const unsigned int        base_element = 0);
 
@@ -325,9 +325,9 @@ namespace internal
        * dimensional element by a tensor product and that the zeroth shape
        * function in zero evaluates to one.
        */
-      template <int dim>
+      template <int dim, int dim_q>
       void
-      reinit(const Quadrature<1> &     quad,
+      reinit(const Quadrature<dim_q> & quad,
              const FiniteElement<dim> &fe_dim,
              const unsigned int        base_element = 0);
 
@@ -519,8 +519,8 @@ namespace internal
     // ------------------------------------------ inline functions
 
     template <typename Number>
-    template <int dim>
-    inline ShapeInfo<Number>::ShapeInfo(const Quadrature<1> &     quad,
+    template <int dim, int dim_q>
+    inline ShapeInfo<Number>::ShapeInfo(const Quadrature<dim_q> & quad,
                                         const FiniteElement<dim> &fe_in,
                                         const unsigned int base_element_number)
       : element_type(tensor_general)
index 83e73620986b67f8b5c7036cbc5fbb722bce3c60..873300c9e5e0cbfc48047e1da0c57908d93357fa 100644 (file)
@@ -128,12 +128,15 @@ namespace internal
 
 
     template <typename Number>
-    template <int dim>
+    template <int dim, int dim_q>
     void
-    ShapeInfo<Number>::reinit(const Quadrature<1> &     quad,
+    ShapeInfo<Number>::reinit(const Quadrature<dim_q> & quad_in,
                               const FiniteElement<dim> &fe_in,
                               const unsigned int        base_element_number)
     {
+      Assert(quad_in.is_tensor_product(), ExcNotImplemented());
+      const auto quad = quad_in.get_tensor_basis()[0];
+
       const FiniteElement<dim> *fe = &fe_in.base_element(base_element_number);
       n_dimensions                 = dim;
       n_components                 = fe_in.n_components();
index 5fae952f73f6032555a8b30b1907be05a082e315..b256093dcb7f438abaea9dd1cc965ceaa3a65fc3 100644 (file)
@@ -40,7 +40,7 @@ for (deal_II_dimension : DIMENSIONS;
         const std::vector<
           const AffineConstraints<deal_II_scalar_vectorized::value_type> *> &,
         const std::vector<IndexSet> &,
-        const std::vector<hp::QCollection<1>> &,
+        const std::vector<hp::QCollection<deal_II_dimension>> &,
         const AdditionalData &);
 
     template const DoFHandler<deal_II_dimension> &
@@ -62,7 +62,7 @@ for (deal_II_dimension : DIMENSIONS;
         const std::vector<const DoFHandler<deal_II_dimension> *> &,
         const std::vector<const AffineConstraints<double> *> &,
         const std::vector<IndexSet> &,
-        const std::vector<hp::QCollection<1>> &,
+        const std::vector<hp::QCollection<deal_II_dimension>> &,
         const AdditionalData &);
   }
 
index 293099d241c7573c6c7cf865325f78ac10c61bfa..30be5ee7aaf869e6228d01639e8030c783d11eb3 100644 (file)
 
 for (deal_II_dimension : DIMENSIONS; deal_II_scalar : REAL_SCALARS)
   {
+    template void internal::MatrixFreeFunctions::ShapeInfo<deal_II_scalar>::
+      reinit<deal_II_dimension>(
+        const Quadrature<deal_II_dimension> &,
+        const FiniteElement<deal_II_dimension, deal_II_dimension> &,
+        const unsigned int);
+
+#if deal_II_dimension > 1
     template void internal::MatrixFreeFunctions::ShapeInfo<deal_II_scalar>::
       reinit<deal_II_dimension>(
         const Quadrature<1> &,
         const FiniteElement<deal_II_dimension, deal_II_dimension> &,
         const unsigned int);
+#endif
   }
 
 for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
@@ -37,11 +45,19 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
 for (deal_II_dimension : DIMENSIONS;
      deal_II_scalar_vectorized : REAL_SCALARS_VECTORIZED)
   {
+    template void internal::MatrixFreeFunctions::
+      ShapeInfo<deal_II_scalar_vectorized>::reinit<deal_II_dimension>(
+        const Quadrature<deal_II_dimension> &,
+        const FiniteElement<deal_II_dimension, deal_II_dimension> &,
+        const unsigned int);
+
+#if deal_II_dimension > 1
     template void internal::MatrixFreeFunctions::
       ShapeInfo<deal_II_scalar_vectorized>::reinit<deal_II_dimension>(
         const Quadrature<1> &,
         const FiniteElement<deal_II_dimension, deal_II_dimension> &,
         const unsigned int);
+#endif
   }
 
 for (deal_II_scalar_vectorized : REAL_SCALARS_VECTORIZED)

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.