From: Niklas Wik Date: Tue, 7 Jun 2022 16:05:00 +0000 (+0200) Subject: Auto-detect update_jacobian_grads X-Git-Tag: v9.5.0-rc1~1174^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f7b8952c60683052e652cf5294d4240e1cb2dcac;p=dealii.git 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 ecabdcd953..f9a8c73a80 100644 --- 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 e9178bb5e2..e514424152 100644 --- 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::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 18bbba3462..67fd4ba04c 100644 --- 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> &quads = - std::vector>()); + std::vector>(), + 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 c9009d39d2..62aa471fb8 100644 --- 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::compute_update_flags( const UpdateFlags update_flags, - const std::vector> &quads) + const std::vector> &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 f3eb1b03dd..080a14b563 100644 --- 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 #include #include +#include #include @@ -809,6 +810,23 @@ MatrixFree::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::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; }