Auto-detect update_jacobian_grads

diff --git a/include/deal.II/matrix_free/mapping_info.h b/include/deal.II/matrix_free/mapping_info.h
index ecabdcd953048a6796a132eb0f431b8240c41708..f9a8c73a806996d52619ee6ec3f64a25b0619761 100644 (file)
--- a/include/deal.II/matrix_free/mapping_info.h
+++ b/include/deal.II/matrix_free/mapping_info.h
@@ -74,7 +74,8 @@ namespace internal
         const UpdateFlags update_flags_cells,
         const UpdateFlags update_flags_boundary_faces,
         const UpdateFlags update_flags_inner_faces,
-        const UpdateFlags update_flags_faces_by_cells);
+        const UpdateFlags update_flags_faces_by_cells,
+        const bool        piola_transform = false);
 
       /**
        * Update the information in the given cells and faces that is the

diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h
index e9178bb5e2aedce1175a2907542522bef104dee6..e514424152ff7d95908321ac58ed7bd95b994846 100644 (file)
--- a/include/deal.II/matrix_free/mapping_info.templates.h
+++ b/include/deal.II/matrix_free/mapping_info.templates.h
@@ -68,7 +68,8 @@ namespace internal
       const UpdateFlags update_flags_cells,
       const UpdateFlags update_flags_boundary_faces,
       const UpdateFlags update_flags_inner_faces,
-      const UpdateFlags update_flags_faces_by_cells)
+      const UpdateFlags update_flags_faces_by_cells,
+      const bool        piola_transform)
     {
       clear();
       this->mapping_collection = mapping;
@@ -85,15 +86,18 @@ namespace internal
       // the mapping that are independent of the FE
       this->update_flags_cells =
         MappingInfoStorage<dim, dim, VectorizedArrayType>::compute_update_flags(
-          update_flags_cells, quad);
+          update_flags_cells, quad, piola_transform);
 
       this->update_flags_boundary_faces =
         ((update_flags_inner_faces | update_flags_boundary_faces) &
              update_quadrature_points ?
            update_quadrature_points :
            update_default) |
-        (((update_flags_inner_faces | update_flags_boundary_faces) &
-          (update_jacobian_grads | update_hessians)) != 0u ?
+        ((((update_flags_inner_faces | update_flags_boundary_faces) &
+           (update_jacobian_grads | update_hessians)) != 0u ||
+          (piola_transform &&
+           ((update_flags_inner_faces | update_flags_boundary_faces) &
+            update_gradients) != 0u)) ?
            update_jacobian_grads :
            update_default) |
         update_normal_vectors | update_JxW_values | update_jacobians;

diff --git a/include/deal.II/matrix_free/mapping_info_storage.h b/include/deal.II/matrix_free/mapping_info_storage.h
index 18bbba346251d5870ca0b16438b55322c5861b1a..67fd4ba04c094e631b166446f5dbf47863dc1a72 100644 (file)
--- a/include/deal.II/matrix_free/mapping_info_storage.h
+++ b/include/deal.II/matrix_free/mapping_info_storage.h
@@ -309,7 +309,8 @@ namespace internal
       compute_update_flags(
         const UpdateFlags                                     update_flags,
         const std::vector<dealii::hp::QCollection<spacedim>> &quads =
-          std::vector<dealii::hp::QCollection<spacedim>>());
+          std::vector<dealii::hp::QCollection<spacedim>>(),
+        const bool piola_transform = false);
 
       /**
        * Prints a detailed summary of memory consumption in the different

diff --git a/include/deal.II/matrix_free/mapping_info_storage.templates.h b/include/deal.II/matrix_free/mapping_info_storage.templates.h
index c9009d39d2f5dae14387db44f8299039a13d919d..62aa471fb8cce3dcbd98ac5cbb4eeb51cc7cb63f 100644 (file)
--- a/include/deal.II/matrix_free/mapping_info_storage.templates.h
+++ b/include/deal.II/matrix_free/mapping_info_storage.templates.h
@@ -124,7 +124,8 @@ namespace internal
     UpdateFlags
     MappingInfoStorage<structdim, spacedim, Number>::compute_update_flags(
       const UpdateFlags                                     update_flags,
-      const std::vector<dealii::hp::QCollection<spacedim>> &quads)
+      const std::vector<dealii::hp::QCollection<spacedim>> &quads,
+      const bool                                            piola_transform)
     {
       // this class is build around the evaluation of jacobians, so compute
       // them in any case. The Jacobians will be inverted manually. Since we
@@ -135,7 +136,8 @@ namespace internal
       // Jacobians (these two together will give use the gradients of the
       // inverse Jacobians, which is what we need)
       if ((update_flags & update_hessians) != 0u ||
-          (update_flags & update_jacobian_grads) != 0u)
+          (update_flags & update_jacobian_grads) != 0u ||
+          (piola_transform && (update_flags & update_gradients) != 0u))
         new_flags |= update_jacobian_grads;
 
       if ((update_flags & update_quadrature_points) != 0u)

diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h
index f3eb1b03dd9921b5b2236f3946b010fc4e2ea7cf..080a14b5631631961d77e0b7b64a97e4e736362d 100644 (file)
--- a/include/deal.II/matrix_free/matrix_free.templates.h
+++ b/include/deal.II/matrix_free/matrix_free.templates.h
@@ -35,6 +35,7 @@
 #include <deal.II/fe/fe_dgq.h>
 #include <deal.II/fe/fe_poly.h>
 #include <deal.II/fe/fe_q_dg0.h>
+#include <deal.II/fe/fe_raviart_thomas.h>
 
 #include <deal.II/hp/q_collection.h>
 
@@ -809,6 +810,23 @@ MatrixFree<dim, Number, VectorizedArrayType>::internal_reinit(
             }
         }
 
+      // Will the piola transform be used? If so we need to update
+      // the jacobian gradients in case of update_gradients on general cells.
+      bool piola_transform = false;
+      for (unsigned int no = 0, c = 0; no < dof_handler.size(); ++no)
+        for (unsigned int b = 0;
+             b < dof_handler[no]->get_fe(0).n_base_elements();
+             ++b, ++c)
+          for (unsigned int fe_no = 0;
+               fe_no < dof_handler[no]->get_fe_collection().size();
+               ++fe_no)
+            for (unsigned int nq = 0; nq < quad.size(); ++nq)
+              for (unsigned int q_no = 0; q_no < quad[nq].size(); ++q_no)
+                if (shape_info(c, nq, fe_no, q_no).element_type ==
+                    internal::MatrixFreeFunctions::ElementType::
+                      tensor_raviart_thomas)
+                  piola_transform = true;
+
       mapping_info.initialize(
         dof_handler[0]->get_triangulation(),
         cell_level_index,
@@ -821,7 +839,8 @@ MatrixFree<dim, Number, VectorizedArrayType>::internal_reinit(
         additional_data.mapping_update_flags,
         additional_data.mapping_update_flags_boundary_faces,
         additional_data.mapping_update_flags_inner_faces,
-        additional_data.mapping_update_flags_faces_by_cells);
+        additional_data.mapping_update_flags_faces_by_cells,
+        piola_transform);
 
       mapping_is_initialized = true;
     }