From c2421a907d4e0e845abcc8ab7ea6ba69fe995be8 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sun, 3 Nov 2024 19:40:17 +0100 Subject: [PATCH] FE_Poly: move content from .templates.h to .cc --- include/deal.II/fe/fe_poly.templates.h | 634 ------------------------- source/fe/fe_poly.cc | 600 ++++++++++++++++++++++- 2 files changed, 599 insertions(+), 635 deletions(-) delete mode 100644 include/deal.II/fe/fe_poly.templates.h diff --git a/include/deal.II/fe/fe_poly.templates.h b/include/deal.II/fe/fe_poly.templates.h deleted file mode 100644 index 06ffa6e611..0000000000 --- a/include/deal.II/fe/fe_poly.templates.h +++ /dev/null @@ -1,634 +0,0 @@ -// ------------------------------------------------------------------------ -// -// SPDX-License-Identifier: LGPL-2.1-or-later -// Copyright (C) 2006 - 2024 by the deal.II authors -// -// This file is part of the deal.II library. -// -// Part of the source code is dual licensed under Apache-2.0 WITH -// LLVM-exception OR LGPL-2.1-or-later. Detailed license information -// governing the source code and code contributions can be found in -// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. -// -// ------------------------------------------------------------------------ - -#ifndef dealii_fe_poly_templates_h -#define dealii_fe_poly_templates_h - - -#include - -#include -#include -#include -#include -#include -#include - -#include -#include -#include - - -DEAL_II_NAMESPACE_OPEN - -template -FE_Poly::FE_Poly( - const ScalarPolynomialsBase &poly_space, - const FiniteElementData &fe_data, - const std::vector &restriction_is_additive_flags, - const std::vector &nonzero_components) - : FiniteElement(fe_data, - restriction_is_additive_flags, - nonzero_components) - , poly_space(poly_space.clone()) -{} - - -template -unsigned int -FE_Poly::get_degree() const -{ - return this->degree; -} - - -template -double -FE_Poly::shape_value(const unsigned int i, - const Point &p) const -{ - AssertIndexRange(i, this->n_dofs_per_cell()); - return poly_space->compute_value(i, p); -} - - -template -double -FE_Poly::shape_value_component( - const unsigned int i, - const Point &p, - const unsigned int component) const -{ - (void)component; - AssertIndexRange(i, this->n_dofs_per_cell()); - AssertIndexRange(component, 1); - return poly_space->compute_value(i, p); -} - - - -template -Tensor<1, dim> -FE_Poly::shape_grad(const unsigned int i, - const Point &p) const -{ - AssertIndexRange(i, this->n_dofs_per_cell()); - return poly_space->template compute_derivative<1>(i, p); -} - - - -template -Tensor<1, dim> -FE_Poly::shape_grad_component(const unsigned int i, - const Point &p, - const unsigned int component) const -{ - (void)component; - AssertIndexRange(i, this->n_dofs_per_cell()); - AssertIndexRange(component, 1); - return poly_space->template compute_derivative<1>(i, p); -} - - - -template -Tensor<2, dim> -FE_Poly::shape_grad_grad(const unsigned int i, - const Point &p) const -{ - AssertIndexRange(i, this->n_dofs_per_cell()); - return poly_space->template compute_derivative<2>(i, p); -} - - - -template -Tensor<2, dim> -FE_Poly::shape_grad_grad_component( - const unsigned int i, - const Point &p, - const unsigned int component) const -{ - (void)component; - AssertIndexRange(i, this->n_dofs_per_cell()); - AssertIndexRange(component, 1); - return poly_space->template compute_derivative<2>(i, p); -} - - - -template -Tensor<3, dim> -FE_Poly::shape_3rd_derivative(const unsigned int i, - const Point &p) const -{ - AssertIndexRange(i, this->n_dofs_per_cell()); - return poly_space->template compute_derivative<3>(i, p); -} - - - -template -Tensor<3, dim> -FE_Poly::shape_3rd_derivative_component( - const unsigned int i, - const Point &p, - const unsigned int component) const -{ - (void)component; - AssertIndexRange(i, this->n_dofs_per_cell()); - AssertIndexRange(component, 1); - return poly_space->template compute_derivative<3>(i, p); -} - - - -template -Tensor<4, dim> -FE_Poly::shape_4th_derivative(const unsigned int i, - const Point &p) const -{ - AssertIndexRange(i, this->n_dofs_per_cell()); - return poly_space->template compute_derivative<4>(i, p); -} - - - -template -Tensor<4, dim> -FE_Poly::shape_4th_derivative_component( - const unsigned int i, - const Point &p, - const unsigned int component) const -{ - (void)component; - AssertIndexRange(i, this->n_dofs_per_cell()); - AssertIndexRange(component, 1); - return poly_space->template compute_derivative<4>(i, p); -} - - - -//--------------------------------------------------------------------------- -// Auxiliary functions -//--------------------------------------------------------------------------- - - -template -UpdateFlags -FE_Poly::requires_update_flags(const UpdateFlags flags) const -{ - UpdateFlags out = update_default; - - if (flags & update_values) - out |= update_values; - if (flags & update_gradients) - out |= update_gradients | update_covariant_transformation; - if (flags & update_hessians) - out |= update_hessians | update_covariant_transformation | - update_gradients | update_jacobian_pushed_forward_grads; - if (flags & update_3rd_derivatives) - out |= update_3rd_derivatives | update_covariant_transformation | - update_hessians | update_gradients | - update_jacobian_pushed_forward_grads | - update_jacobian_pushed_forward_2nd_derivatives; - if (flags & update_normal_vectors) - out |= update_normal_vectors | update_JxW_values; - - return out; -} - - - -//--------------------------------------------------------------------------- -// Fill data of FEValues -//--------------------------------------------------------------------------- - - - -/** - * Returns whether we need to correct the Hessians and third derivatives with - * the derivatives of the Jacobian. This is determined by checking if - * the jacobian_pushed_forward are zero. - * - * Especially for the third derivatives, the correction term is very expensive, - * which is why we check if the derivatives are zero before computing the - * correction. - */ -template -bool -higher_derivatives_need_correcting( - const Mapping &mapping, - const internal::FEValuesImplementation::MappingRelatedData - &mapping_data, - const unsigned int n_q_points, - const UpdateFlags update_flags) -{ - // If higher derivatives weren't requested we don't need to correct them. - const bool update_higher_derivatives = - (update_flags & update_hessians) || (update_flags & update_3rd_derivatives); - if (!update_higher_derivatives) - return false; - - // If we have a Cartesian mapping, we know that jacoban_pushed_forward_grads - // are identically zero. - if (dynamic_cast *>(&mapping)) - return false; - - // Here, we should check if jacobian_pushed_forward_grads are zero at the - // quadrature points. This is yet to be implemented. - (void)mapping_data; - (void)n_q_points; - - return true; -} - - - -template -void -FE_Poly::fill_fe_values( - const typename Triangulation::cell_iterator &, - const CellSimilarity::Similarity cell_similarity, - const Quadrature &quadrature, - const Mapping &mapping, - const typename Mapping::InternalDataBase &mapping_internal, - const internal::FEValuesImplementation::MappingRelatedData - &mapping_data, - const typename FiniteElement::InternalDataBase &fe_internal, - dealii::internal::FEValuesImplementation::FiniteElementRelatedData - &output_data) const -{ - // convert data object to internal - // data for this class. fails with - // an exception if that is not - // possible - Assert(dynamic_cast(&fe_internal) != nullptr, - ExcInternalError()); - const InternalData &fe_data = static_cast(fe_internal); - - const UpdateFlags flags(fe_data.update_each); - - const bool need_to_correct_higher_derivatives = - higher_derivatives_need_correcting(mapping, - mapping_data, - quadrature.size(), - flags); - - // transform gradients and higher derivatives. there is nothing to do - // for values since we already emplaced them into output_data when - // we were in get_data() - if ((flags & update_gradients) && - (cell_similarity != CellSimilarity::translation)) - for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k) - mapping.transform(make_array_view(fe_data.shape_gradients, k), - mapping_covariant, - mapping_internal, - make_array_view(output_data.shape_gradients, k)); - - if ((flags & update_hessians) && - (cell_similarity != CellSimilarity::translation)) - { - for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k) - mapping.transform(make_array_view(fe_data.shape_hessians, k), - mapping_covariant_gradient, - mapping_internal, - make_array_view(output_data.shape_hessians, k)); - - if (need_to_correct_higher_derivatives) - correct_hessians(output_data, mapping_data, quadrature.size()); - } - - if ((flags & update_3rd_derivatives) && - (cell_similarity != CellSimilarity::translation)) - { - for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k) - mapping.transform(make_array_view(fe_data.shape_3rd_derivatives, k), - mapping_covariant_hessian, - mapping_internal, - make_array_view(output_data.shape_3rd_derivatives, - k)); - - if (need_to_correct_higher_derivatives) - correct_third_derivatives(output_data, mapping_data, quadrature.size()); - } -} - - - -template -void -FE_Poly::fill_fe_face_values( - const typename Triangulation::cell_iterator &cell, - const unsigned int face_no, - const hp::QCollection &quadrature, - const Mapping &mapping, - const typename Mapping::InternalDataBase &mapping_internal, - const internal::FEValuesImplementation::MappingRelatedData - &mapping_data, - const typename FiniteElement::InternalDataBase &fe_internal, - dealii::internal::FEValuesImplementation::FiniteElementRelatedData - &output_data) const -{ - const unsigned int n_q_points = - quadrature[quadrature.size() == 1 ? 0 : face_no].size(); - - // convert data object to internal - // data for this class. fails with - // an exception if that is not - // possible - Assert(dynamic_cast(&fe_internal) != nullptr, - ExcInternalError()); - const InternalData &fe_data = static_cast(fe_internal); - - // offset determines which data set - // to take (all data sets for all - // faces are stored contiguously) - - const auto offset = - QProjector::DataSetDescriptor::face(this->reference_cell(), - face_no, - cell->face_orientation(face_no), - cell->face_flip(face_no), - cell->face_rotation(face_no), - quadrature); - - const UpdateFlags flags(fe_data.update_each); - - const bool need_to_correct_higher_derivatives = - higher_derivatives_need_correcting(mapping, - mapping_data, - n_q_points, - flags); - - // transform gradients and higher derivatives. we also have to copy - // the values (unlike in the case of fill_fe_values()) since - // we need to take into account the offsets - if (flags & update_values) - for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k) - for (unsigned int i = 0; i < n_q_points; ++i) - output_data.shape_values(k, i) = fe_data.shape_values[k][i + offset]; - - if (flags & update_gradients) - for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k) - mapping.transform( - make_array_view(fe_data.shape_gradients, k, offset, n_q_points), - mapping_covariant, - mapping_internal, - make_array_view(output_data.shape_gradients, k)); - - if (flags & update_hessians) - { - for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k) - mapping.transform( - make_array_view(fe_data.shape_hessians, k, offset, n_q_points), - mapping_covariant_gradient, - mapping_internal, - make_array_view(output_data.shape_hessians, k)); - - if (need_to_correct_higher_derivatives) - correct_hessians(output_data, mapping_data, n_q_points); - } - - if (flags & update_3rd_derivatives) - { - for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k) - mapping.transform( - make_array_view(fe_data.shape_3rd_derivatives, k, offset, n_q_points), - mapping_covariant_hessian, - mapping_internal, - make_array_view(output_data.shape_3rd_derivatives, k)); - - if (need_to_correct_higher_derivatives) - correct_third_derivatives(output_data, mapping_data, n_q_points); - } -} - - - -template -void -FE_Poly::fill_fe_subface_values( - const typename Triangulation::cell_iterator &cell, - const unsigned int face_no, - const unsigned int sub_no, - const Quadrature &quadrature, - const Mapping &mapping, - const typename Mapping::InternalDataBase &mapping_internal, - const internal::FEValuesImplementation::MappingRelatedData - &mapping_data, - const typename FiniteElement::InternalDataBase &fe_internal, - dealii::internal::FEValuesImplementation::FiniteElementRelatedData - &output_data) const -{ - // convert data object to internal - // data for this class. fails with - // an exception if that is not - // possible - Assert(dynamic_cast(&fe_internal) != nullptr, - ExcInternalError()); - const InternalData &fe_data = static_cast(fe_internal); - - // offset determines which data set - // to take (all data sets for all - // sub-faces are stored contiguously) - - const auto offset = - QProjector::DataSetDescriptor::subface(this->reference_cell(), - face_no, - sub_no, - cell->face_orientation(face_no), - cell->face_flip(face_no), - cell->face_rotation(face_no), - quadrature.size(), - cell->subface_case(face_no)); - - const UpdateFlags flags(fe_data.update_each); - - const bool need_to_correct_higher_derivatives = - higher_derivatives_need_correcting(mapping, - mapping_data, - quadrature.size(), - flags); - - // transform gradients and higher derivatives. we also have to copy - // the values (unlike in the case of fill_fe_values()) since - // we need to take into account the offsets - if (flags & update_values) - for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k) - for (unsigned int i = 0; i < quadrature.size(); ++i) - output_data.shape_values(k, i) = fe_data.shape_values[k][i + offset]; - - if (flags & update_gradients) - for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k) - mapping.transform( - make_array_view(fe_data.shape_gradients, k, offset, quadrature.size()), - mapping_covariant, - mapping_internal, - make_array_view(output_data.shape_gradients, k)); - - if (flags & update_hessians) - { - for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k) - mapping.transform( - make_array_view(fe_data.shape_hessians, k, offset, quadrature.size()), - mapping_covariant_gradient, - mapping_internal, - make_array_view(output_data.shape_hessians, k)); - - if (need_to_correct_higher_derivatives) - correct_hessians(output_data, mapping_data, quadrature.size()); - } - - if (flags & update_3rd_derivatives) - { - for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k) - mapping.transform(make_array_view(fe_data.shape_3rd_derivatives, - k, - offset, - quadrature.size()), - mapping_covariant_hessian, - mapping_internal, - make_array_view(output_data.shape_3rd_derivatives, - k)); - - if (need_to_correct_higher_derivatives) - correct_third_derivatives(output_data, mapping_data, quadrature.size()); - } -} - - - -template -inline void -FE_Poly::correct_hessians( - internal::FEValuesImplementation::FiniteElementRelatedData - &output_data, - const internal::FEValuesImplementation::MappingRelatedData - &mapping_data, - const unsigned int n_q_points) const -{ - for (unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof) - for (unsigned int i = 0; i < n_q_points; ++i) - for (unsigned int j = 0; j < spacedim; ++j) - output_data.shape_hessians[dof][i] -= - mapping_data.jacobian_pushed_forward_grads[i][j] * - output_data.shape_gradients[dof][i][j]; -} - - - -template -inline void -FE_Poly::correct_third_derivatives( - internal::FEValuesImplementation::FiniteElementRelatedData - &output_data, - const internal::FEValuesImplementation::MappingRelatedData - &mapping_data, - const unsigned int n_q_points) const -{ - for (unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof) - for (unsigned int i = 0; i < n_q_points; ++i) - for (unsigned int j = 0; j < spacedim; ++j) - for (unsigned int k = 0; k < spacedim; ++k) - for (unsigned int l = 0; l < spacedim; ++l) - for (unsigned int m = 0; m < spacedim; ++m) - { - output_data.shape_3rd_derivatives[dof][i][j][k][l] -= - (mapping_data.jacobian_pushed_forward_grads[i][m][j][l] * - output_data.shape_hessians[dof][i][k][m]) + - (mapping_data.jacobian_pushed_forward_grads[i][m][k][l] * - output_data.shape_hessians[dof][i][j][m]) + - (mapping_data.jacobian_pushed_forward_grads[i][m][j][k] * - output_data.shape_hessians[dof][i][l][m]) + - (mapping_data - .jacobian_pushed_forward_2nd_derivatives[i][m][j][k][l] * - output_data.shape_gradients[dof][i][m]); - } -} - - - -template -inline const ScalarPolynomialsBase & -FE_Poly::get_poly_space() const -{ - return *poly_space; -} - - - -template -std::vector -FE_Poly::get_poly_space_numbering() const -{ - auto *const space_tensor_prod = - dynamic_cast *>(this->poly_space.get()); - if (space_tensor_prod != nullptr) - return space_tensor_prod->get_numbering(); - - auto *const space_tensor_prod_aniso = - dynamic_cast *>(this->poly_space.get()); - if (space_tensor_prod_aniso != nullptr) - return space_tensor_prod_aniso->get_numbering(); - - auto *const space_tensor_prod_piecewise = dynamic_cast< - TensorProductPolynomials> *>( - this->poly_space.get()); - if (space_tensor_prod_piecewise != nullptr) - return space_tensor_prod_piecewise->get_numbering(); - - auto *const space_tensor_prod_bubbles = - dynamic_cast *>( - this->poly_space.get()); - if (space_tensor_prod_bubbles != nullptr) - return space_tensor_prod_bubbles->get_numbering(); - - auto *const space_tensor_prod_const = - dynamic_cast *>(this->poly_space.get()); - if (space_tensor_prod_const != nullptr) - return space_tensor_prod_const->get_numbering(); - - DEAL_II_NOT_IMPLEMENTED(); - return std::vector(); -} - - - -template -std::vector -FE_Poly::get_poly_space_numbering_inverse() const -{ - return Utilities::invert_permutation(get_poly_space_numbering()); -} - - - -template -std::size_t -FE_Poly::memory_consumption() const -{ - return FiniteElement::memory_consumption() + - poly_space->memory_consumption(); -} - - - -DEAL_II_NAMESPACE_CLOSE - -#endif diff --git a/source/fe/fe_poly.cc b/source/fe/fe_poly.cc index 9de0950e5d..6a164dbf82 100644 --- a/source/fe/fe_poly.cc +++ b/source/fe/fe_poly.cc @@ -13,6 +13,9 @@ // ------------------------------------------------------------------------ +#include + +#include #include #include #include @@ -22,8 +25,8 @@ #include #include -#include #include +#include DEAL_II_NAMESPACE_OPEN @@ -35,6 +38,601 @@ FE_Poly::FE_Poly(const FE_Poly &fe) , poly_space(fe.poly_space->clone()) {} +template +FE_Poly::FE_Poly( + const ScalarPolynomialsBase &poly_space, + const FiniteElementData &fe_data, + const std::vector &restriction_is_additive_flags, + const std::vector &nonzero_components) + : FiniteElement(fe_data, + restriction_is_additive_flags, + nonzero_components) + , poly_space(poly_space.clone()) +{} + + +template +unsigned int +FE_Poly::get_degree() const +{ + return this->degree; +} + + +template +double +FE_Poly::shape_value(const unsigned int i, + const Point &p) const +{ + AssertIndexRange(i, this->n_dofs_per_cell()); + return poly_space->compute_value(i, p); +} + + +template +double +FE_Poly::shape_value_component( + const unsigned int i, + const Point &p, + const unsigned int component) const +{ + (void)component; + AssertIndexRange(i, this->n_dofs_per_cell()); + AssertIndexRange(component, 1); + return poly_space->compute_value(i, p); +} + + + +template +Tensor<1, dim> +FE_Poly::shape_grad(const unsigned int i, + const Point &p) const +{ + AssertIndexRange(i, this->n_dofs_per_cell()); + return poly_space->template compute_derivative<1>(i, p); +} + + + +template +Tensor<1, dim> +FE_Poly::shape_grad_component(const unsigned int i, + const Point &p, + const unsigned int component) const +{ + (void)component; + AssertIndexRange(i, this->n_dofs_per_cell()); + AssertIndexRange(component, 1); + return poly_space->template compute_derivative<1>(i, p); +} + + + +template +Tensor<2, dim> +FE_Poly::shape_grad_grad(const unsigned int i, + const Point &p) const +{ + AssertIndexRange(i, this->n_dofs_per_cell()); + return poly_space->template compute_derivative<2>(i, p); +} + + + +template +Tensor<2, dim> +FE_Poly::shape_grad_grad_component( + const unsigned int i, + const Point &p, + const unsigned int component) const +{ + (void)component; + AssertIndexRange(i, this->n_dofs_per_cell()); + AssertIndexRange(component, 1); + return poly_space->template compute_derivative<2>(i, p); +} + + + +template +Tensor<3, dim> +FE_Poly::shape_3rd_derivative(const unsigned int i, + const Point &p) const +{ + AssertIndexRange(i, this->n_dofs_per_cell()); + return poly_space->template compute_derivative<3>(i, p); +} + + + +template +Tensor<3, dim> +FE_Poly::shape_3rd_derivative_component( + const unsigned int i, + const Point &p, + const unsigned int component) const +{ + (void)component; + AssertIndexRange(i, this->n_dofs_per_cell()); + AssertIndexRange(component, 1); + return poly_space->template compute_derivative<3>(i, p); +} + + + +template +Tensor<4, dim> +FE_Poly::shape_4th_derivative(const unsigned int i, + const Point &p) const +{ + AssertIndexRange(i, this->n_dofs_per_cell()); + return poly_space->template compute_derivative<4>(i, p); +} + + + +template +Tensor<4, dim> +FE_Poly::shape_4th_derivative_component( + const unsigned int i, + const Point &p, + const unsigned int component) const +{ + (void)component; + AssertIndexRange(i, this->n_dofs_per_cell()); + AssertIndexRange(component, 1); + return poly_space->template compute_derivative<4>(i, p); +} + + + +//--------------------------------------------------------------------------- +// Auxiliary functions +//--------------------------------------------------------------------------- + + +template +UpdateFlags +FE_Poly::requires_update_flags(const UpdateFlags flags) const +{ + UpdateFlags out = update_default; + + if (flags & update_values) + out |= update_values; + if (flags & update_gradients) + out |= update_gradients | update_covariant_transformation; + if (flags & update_hessians) + out |= update_hessians | update_covariant_transformation | + update_gradients | update_jacobian_pushed_forward_grads; + if (flags & update_3rd_derivatives) + out |= update_3rd_derivatives | update_covariant_transformation | + update_hessians | update_gradients | + update_jacobian_pushed_forward_grads | + update_jacobian_pushed_forward_2nd_derivatives; + if (flags & update_normal_vectors) + out |= update_normal_vectors | update_JxW_values; + + return out; +} + + + +//--------------------------------------------------------------------------- +// Fill data of FEValues +//--------------------------------------------------------------------------- + + + +/** + * Returns whether we need to correct the Hessians and third derivatives with + * the derivatives of the Jacobian. This is determined by checking if + * the jacobian_pushed_forward are zero. + * + * Especially for the third derivatives, the correction term is very expensive, + * which is why we check if the derivatives are zero before computing the + * correction. + */ +template +bool +higher_derivatives_need_correcting( + const Mapping &mapping, + const internal::FEValuesImplementation::MappingRelatedData + &mapping_data, + const unsigned int n_q_points, + const UpdateFlags update_flags) +{ + // If higher derivatives weren't requested we don't need to correct them. + const bool update_higher_derivatives = + (update_flags & update_hessians) || (update_flags & update_3rd_derivatives); + if (!update_higher_derivatives) + return false; + + // If we have a Cartesian mapping, we know that jacoban_pushed_forward_grads + // are identically zero. + if (dynamic_cast *>(&mapping)) + return false; + + // Here, we should check if jacobian_pushed_forward_grads are zero at the + // quadrature points. This is yet to be implemented. + (void)mapping_data; + (void)n_q_points; + + return true; +} + + + +template +void +FE_Poly::fill_fe_values( + const typename Triangulation::cell_iterator &, + const CellSimilarity::Similarity cell_similarity, + const Quadrature &quadrature, + const Mapping &mapping, + const typename Mapping::InternalDataBase &mapping_internal, + const internal::FEValuesImplementation::MappingRelatedData + &mapping_data, + const typename FiniteElement::InternalDataBase &fe_internal, + dealii::internal::FEValuesImplementation::FiniteElementRelatedData + &output_data) const +{ + // convert data object to internal + // data for this class. fails with + // an exception if that is not + // possible + Assert(dynamic_cast(&fe_internal) != nullptr, + ExcInternalError()); + const InternalData &fe_data = static_cast(fe_internal); + + const UpdateFlags flags(fe_data.update_each); + + const bool need_to_correct_higher_derivatives = + higher_derivatives_need_correcting(mapping, + mapping_data, + quadrature.size(), + flags); + + // transform gradients and higher derivatives. there is nothing to do + // for values since we already emplaced them into output_data when + // we were in get_data() + if ((flags & update_gradients) && + (cell_similarity != CellSimilarity::translation)) + for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k) + mapping.transform(make_array_view(fe_data.shape_gradients, k), + mapping_covariant, + mapping_internal, + make_array_view(output_data.shape_gradients, k)); + + if ((flags & update_hessians) && + (cell_similarity != CellSimilarity::translation)) + { + for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k) + mapping.transform(make_array_view(fe_data.shape_hessians, k), + mapping_covariant_gradient, + mapping_internal, + make_array_view(output_data.shape_hessians, k)); + + if (need_to_correct_higher_derivatives) + correct_hessians(output_data, mapping_data, quadrature.size()); + } + + if ((flags & update_3rd_derivatives) && + (cell_similarity != CellSimilarity::translation)) + { + for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k) + mapping.transform(make_array_view(fe_data.shape_3rd_derivatives, k), + mapping_covariant_hessian, + mapping_internal, + make_array_view(output_data.shape_3rd_derivatives, + k)); + + if (need_to_correct_higher_derivatives) + correct_third_derivatives(output_data, mapping_data, quadrature.size()); + } +} + + + +template +void +FE_Poly::fill_fe_face_values( + const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const hp::QCollection &quadrature, + const Mapping &mapping, + const typename Mapping::InternalDataBase &mapping_internal, + const internal::FEValuesImplementation::MappingRelatedData + &mapping_data, + const typename FiniteElement::InternalDataBase &fe_internal, + dealii::internal::FEValuesImplementation::FiniteElementRelatedData + &output_data) const +{ + const unsigned int n_q_points = + quadrature[quadrature.size() == 1 ? 0 : face_no].size(); + + // convert data object to internal + // data for this class. fails with + // an exception if that is not + // possible + Assert(dynamic_cast(&fe_internal) != nullptr, + ExcInternalError()); + const InternalData &fe_data = static_cast(fe_internal); + + // offset determines which data set + // to take (all data sets for all + // faces are stored contiguously) + + const auto offset = + QProjector::DataSetDescriptor::face(this->reference_cell(), + face_no, + cell->face_orientation(face_no), + cell->face_flip(face_no), + cell->face_rotation(face_no), + quadrature); + + const UpdateFlags flags(fe_data.update_each); + + const bool need_to_correct_higher_derivatives = + higher_derivatives_need_correcting(mapping, + mapping_data, + n_q_points, + flags); + + // transform gradients and higher derivatives. we also have to copy + // the values (unlike in the case of fill_fe_values()) since + // we need to take into account the offsets + if (flags & update_values) + for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k) + for (unsigned int i = 0; i < n_q_points; ++i) + output_data.shape_values(k, i) = fe_data.shape_values[k][i + offset]; + + if (flags & update_gradients) + for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k) + mapping.transform( + make_array_view(fe_data.shape_gradients, k, offset, n_q_points), + mapping_covariant, + mapping_internal, + make_array_view(output_data.shape_gradients, k)); + + if (flags & update_hessians) + { + for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k) + mapping.transform( + make_array_view(fe_data.shape_hessians, k, offset, n_q_points), + mapping_covariant_gradient, + mapping_internal, + make_array_view(output_data.shape_hessians, k)); + + if (need_to_correct_higher_derivatives) + correct_hessians(output_data, mapping_data, n_q_points); + } + + if (flags & update_3rd_derivatives) + { + for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k) + mapping.transform( + make_array_view(fe_data.shape_3rd_derivatives, k, offset, n_q_points), + mapping_covariant_hessian, + mapping_internal, + make_array_view(output_data.shape_3rd_derivatives, k)); + + if (need_to_correct_higher_derivatives) + correct_third_derivatives(output_data, mapping_data, n_q_points); + } +} + + + +template +void +FE_Poly::fill_fe_subface_values( + const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const unsigned int sub_no, + const Quadrature &quadrature, + const Mapping &mapping, + const typename Mapping::InternalDataBase &mapping_internal, + const internal::FEValuesImplementation::MappingRelatedData + &mapping_data, + const typename FiniteElement::InternalDataBase &fe_internal, + dealii::internal::FEValuesImplementation::FiniteElementRelatedData + &output_data) const +{ + // convert data object to internal + // data for this class. fails with + // an exception if that is not + // possible + Assert(dynamic_cast(&fe_internal) != nullptr, + ExcInternalError()); + const InternalData &fe_data = static_cast(fe_internal); + + // offset determines which data set + // to take (all data sets for all + // sub-faces are stored contiguously) + + const auto offset = + QProjector::DataSetDescriptor::subface(this->reference_cell(), + face_no, + sub_no, + cell->face_orientation(face_no), + cell->face_flip(face_no), + cell->face_rotation(face_no), + quadrature.size(), + cell->subface_case(face_no)); + + const UpdateFlags flags(fe_data.update_each); + + const bool need_to_correct_higher_derivatives = + higher_derivatives_need_correcting(mapping, + mapping_data, + quadrature.size(), + flags); + + // transform gradients and higher derivatives. we also have to copy + // the values (unlike in the case of fill_fe_values()) since + // we need to take into account the offsets + if (flags & update_values) + for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k) + for (unsigned int i = 0; i < quadrature.size(); ++i) + output_data.shape_values(k, i) = fe_data.shape_values[k][i + offset]; + + if (flags & update_gradients) + for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k) + mapping.transform( + make_array_view(fe_data.shape_gradients, k, offset, quadrature.size()), + mapping_covariant, + mapping_internal, + make_array_view(output_data.shape_gradients, k)); + + if (flags & update_hessians) + { + for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k) + mapping.transform( + make_array_view(fe_data.shape_hessians, k, offset, quadrature.size()), + mapping_covariant_gradient, + mapping_internal, + make_array_view(output_data.shape_hessians, k)); + + if (need_to_correct_higher_derivatives) + correct_hessians(output_data, mapping_data, quadrature.size()); + } + + if (flags & update_3rd_derivatives) + { + for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k) + mapping.transform(make_array_view(fe_data.shape_3rd_derivatives, + k, + offset, + quadrature.size()), + mapping_covariant_hessian, + mapping_internal, + make_array_view(output_data.shape_3rd_derivatives, + k)); + + if (need_to_correct_higher_derivatives) + correct_third_derivatives(output_data, mapping_data, quadrature.size()); + } +} + + + +template +inline void +FE_Poly::correct_hessians( + internal::FEValuesImplementation::FiniteElementRelatedData + &output_data, + const internal::FEValuesImplementation::MappingRelatedData + &mapping_data, + const unsigned int n_q_points) const +{ + for (unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof) + for (unsigned int i = 0; i < n_q_points; ++i) + for (unsigned int j = 0; j < spacedim; ++j) + output_data.shape_hessians[dof][i] -= + mapping_data.jacobian_pushed_forward_grads[i][j] * + output_data.shape_gradients[dof][i][j]; +} + + + +template +inline void +FE_Poly::correct_third_derivatives( + internal::FEValuesImplementation::FiniteElementRelatedData + &output_data, + const internal::FEValuesImplementation::MappingRelatedData + &mapping_data, + const unsigned int n_q_points) const +{ + for (unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof) + for (unsigned int i = 0; i < n_q_points; ++i) + for (unsigned int j = 0; j < spacedim; ++j) + for (unsigned int k = 0; k < spacedim; ++k) + for (unsigned int l = 0; l < spacedim; ++l) + for (unsigned int m = 0; m < spacedim; ++m) + { + output_data.shape_3rd_derivatives[dof][i][j][k][l] -= + (mapping_data.jacobian_pushed_forward_grads[i][m][j][l] * + output_data.shape_hessians[dof][i][k][m]) + + (mapping_data.jacobian_pushed_forward_grads[i][m][k][l] * + output_data.shape_hessians[dof][i][j][m]) + + (mapping_data.jacobian_pushed_forward_grads[i][m][j][k] * + output_data.shape_hessians[dof][i][l][m]) + + (mapping_data + .jacobian_pushed_forward_2nd_derivatives[i][m][j][k][l] * + output_data.shape_gradients[dof][i][m]); + } +} + + + +template +inline const ScalarPolynomialsBase & +FE_Poly::get_poly_space() const +{ + return *poly_space; +} + + + +template +std::vector +FE_Poly::get_poly_space_numbering() const +{ + auto *const space_tensor_prod = + dynamic_cast *>(this->poly_space.get()); + if (space_tensor_prod != nullptr) + return space_tensor_prod->get_numbering(); + + auto *const space_tensor_prod_aniso = + dynamic_cast *>(this->poly_space.get()); + if (space_tensor_prod_aniso != nullptr) + return space_tensor_prod_aniso->get_numbering(); + + auto *const space_tensor_prod_piecewise = dynamic_cast< + TensorProductPolynomials> *>( + this->poly_space.get()); + if (space_tensor_prod_piecewise != nullptr) + return space_tensor_prod_piecewise->get_numbering(); + + auto *const space_tensor_prod_bubbles = + dynamic_cast *>( + this->poly_space.get()); + if (space_tensor_prod_bubbles != nullptr) + return space_tensor_prod_bubbles->get_numbering(); + + auto *const space_tensor_prod_const = + dynamic_cast *>(this->poly_space.get()); + if (space_tensor_prod_const != nullptr) + return space_tensor_prod_const->get_numbering(); + + DEAL_II_NOT_IMPLEMENTED(); + return std::vector(); +} + + + +template +std::vector +FE_Poly::get_poly_space_numbering_inverse() const +{ + return Utilities::invert_permutation(get_poly_space_numbering()); +} + + + +template +std::size_t +FE_Poly::memory_consumption() const +{ + return FiniteElement::memory_consumption() + + poly_space->memory_consumption(); +} + #endif #include "fe_poly.inst" -- 2.39.5