From: Marco Feder Date: Sat, 29 Oct 2022 16:51:30 +0000 (+0200) Subject: Allow construction of FEImmersedSurfaceValues with a MappingFEField X-Git-Tag: v9.5.0-rc1~887^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a3d2ec1eefd0380c5ec754dcf301d6749c77cc88;p=dealii.git Allow construction of FEImmersedSurfaceValues with a MappingFEField --- diff --git a/doc/news/changes/minor/20221031Feder b/doc/news/changes/minor/20221031Feder new file mode 100644 index 0000000000..a01350b990 --- /dev/null +++ b/doc/news/changes/minor/20221031Feder @@ -0,0 +1,4 @@ +Improved: FEImmersedSurfaceValues can now be constructed using a MappingFEField. +
+(Marco Feder, 2022/10/31) + diff --git a/include/deal.II/fe/mapping_fe_field.h b/include/deal.II/fe/mapping_fe_field.h index 9562157930..676369469d 100644 --- a/include/deal.II/fe/mapping_fe_field.h +++ b/include/deal.II/fe/mapping_fe_field.h @@ -534,6 +534,14 @@ protected: internal::FEValuesImplementation::MappingRelatedData &output_data) const override; + virtual void + fill_fe_immersed_surface_values( + const typename Triangulation::cell_iterator &cell, + const NonMatching::ImmersedSurfaceQuadrature & quadrature, + const typename Mapping::InternalDataBase & internal_data, + dealii::internal::FEValuesImplementation::MappingRelatedData + &output_data) const override; + /** * @} */ diff --git a/include/deal.II/non_matching/fe_immersed_values.h b/include/deal.II/non_matching/fe_immersed_values.h index 8ee37ac395..4889720678 100644 --- a/include/deal.II/non_matching/fe_immersed_values.h +++ b/include/deal.II/non_matching/fe_immersed_values.h @@ -60,8 +60,8 @@ namespace NonMatching * Constructor. Gets cell-independent data from mapping and finite element * objects, matching the quadrature rule and update flags. * - * @note Currently this class is only implemented for MappingCartesian and - * MappingQ. + * @note Currently this class is only implemented for MappingCartesian, + * MappingQ and MappingFEField. */ FEImmersedSurfaceValues(const Mapping & mapping, const FiniteElement & element, diff --git a/source/fe/mapping_fe_field.cc b/source/fe/mapping_fe_field.cc index 749783d9dc..58be345c39 100644 --- a/source/fe/mapping_fe_field.cc +++ b/source/fe/mapping_fe_field.cc @@ -1771,6 +1771,170 @@ MappingFEField::fill_fe_subface_values( } + +template +void +MappingFEField::fill_fe_immersed_surface_values( + const typename Triangulation::cell_iterator &cell, + const NonMatching::ImmersedSurfaceQuadrature & quadrature, + const typename Mapping::InternalDataBase & internal_data, + dealii::internal::FEValuesImplementation::MappingRelatedData + &output_data) const +{ + AssertDimension(dim, spacedim); + Assert(dynamic_cast(&internal_data) != nullptr, + ExcInternalError()); + const InternalData &data = static_cast(internal_data); + + const unsigned int n_q_points = quadrature.size(); + + update_internal_dofs(cell, data); + + internal::MappingFEFieldImplementation:: + maybe_compute_q_points( + QProjector::DataSetDescriptor::cell(), + data, + euler_dof_handler->get_fe(), + fe_mask, + fe_to_real, + output_data.quadrature_points); + + internal::MappingFEFieldImplementation:: + maybe_update_Jacobians( + QProjector::DataSetDescriptor::cell(), + data, + euler_dof_handler->get_fe(), + fe_mask, + fe_to_real); + + const UpdateFlags update_flags = data.update_each; + const std::vector &weights = quadrature.get_weights(); + + if (update_flags & (update_normal_vectors | update_JxW_values)) + { + AssertDimension(output_data.JxW_values.size(), n_q_points); + + Assert(!(update_flags & update_normal_vectors) || + (output_data.normal_vectors.size() == n_q_points), + ExcDimensionMismatch(output_data.normal_vectors.size(), + n_q_points)); + + + for (unsigned int point = 0; point < n_q_points; ++point) + { + const double det = data.contravariant[point].determinant(); + + // check for distorted cells. + + // TODO: this allows for anisotropies of up to 1e6 in 3D and + // 1e12 in 2D. might want to find a finer + // (dimension-independent) criterion + Assert(det > 1e-12 * Utilities::fixed_power( + cell->diameter() / std::sqrt(double(dim))), + (typename Mapping::ExcDistortedMappedCell( + cell->center(), det, point))); + + // The normals are n = J^{-T} * \hat{n} before normalizing. + Tensor<1, spacedim> normal; + for (unsigned int d = 0; d < spacedim; d++) + normal[d] = + data.covariant[point][d] * quadrature.normal_vector(point); + + output_data.JxW_values[point] = weights[point] * det * normal.norm(); + + if ((update_flags & update_normal_vectors) != 0u) + { + normal /= normal.norm(); + output_data.normal_vectors[point] = normal; + } + } + + // copy values from InternalData to vector given by reference + if (update_flags & update_jacobians) + { + AssertDimension(output_data.jacobians.size(), n_q_points); + for (unsigned int point = 0; point < n_q_points; ++point) + output_data.jacobians[point] = data.contravariant[point]; + } + + // copy values from InternalData to vector given by reference + if (update_flags & update_inverse_jacobians) + { + AssertDimension(output_data.inverse_jacobians.size(), n_q_points); + for (unsigned int point = 0; point < n_q_points; ++point) + output_data.inverse_jacobians[point] = + data.covariant[point].transpose(); + } + + // calculate derivatives of the Jacobians + internal::MappingFEFieldImplementation:: + maybe_update_jacobian_grads( + QProjector::DataSetDescriptor::cell(), + data, + euler_dof_handler->get_fe(), + fe_mask, + fe_to_real, + output_data.jacobian_grads); + + // calculate derivatives of the Jacobians pushed forward to real cell + // coordinates + internal::MappingFEFieldImplementation:: + maybe_update_jacobian_pushed_forward_grads( + QProjector::DataSetDescriptor::cell(), + data, + euler_dof_handler->get_fe(), + fe_mask, + fe_to_real, + output_data.jacobian_pushed_forward_grads); + + // calculate hessians of the Jacobians + internal::MappingFEFieldImplementation:: + maybe_update_jacobian_2nd_derivatives( + QProjector::DataSetDescriptor::cell(), + data, + euler_dof_handler->get_fe(), + fe_mask, + fe_to_real, + output_data.jacobian_2nd_derivatives); + + // calculate hessians of the Jacobians pushed forward to real cell + // coordinates + internal::MappingFEFieldImplementation:: + maybe_update_jacobian_pushed_forward_2nd_derivatives( + QProjector::DataSetDescriptor::cell(), + data, + euler_dof_handler->get_fe(), + fe_mask, + fe_to_real, + output_data.jacobian_pushed_forward_2nd_derivatives); + + // calculate gradients of the hessians of the Jacobians + internal::MappingFEFieldImplementation:: + maybe_update_jacobian_3rd_derivatives( + QProjector::DataSetDescriptor::cell(), + data, + euler_dof_handler->get_fe(), + fe_mask, + fe_to_real, + output_data.jacobian_3rd_derivatives); + + // calculate gradients of the hessians of the Jacobians pushed forward to + // real cell coordinates + internal::MappingFEFieldImplementation:: + maybe_update_jacobian_pushed_forward_3rd_derivatives( + QProjector::DataSetDescriptor::cell(), + data, + euler_dof_handler->get_fe(), + fe_mask, + fe_to_real, + output_data.jacobian_pushed_forward_3rd_derivatives); + } +} + namespace internal { namespace MappingFEFieldImplementation diff --git a/tests/non_matching/fe_immersed_surface_values_02.cc b/tests/non_matching/fe_immersed_surface_values_02.cc new file mode 100644 index 0000000000..3affc77db7 --- /dev/null +++ b/tests/non_matching/fe_immersed_surface_values_02.cc @@ -0,0 +1,292 @@ +#include +#include +#include + +#include + +#include +#include +#include + +#include +#include +#include + +#include + +#include +#include + +#include + +#include "../tests.h" + + +using namespace dealii; +using NonMatching::FEImmersedSurfaceValues; + +// Set up a triangulation that we can construct an +// FEImmersedSurfaceValues object and compute known values when using a +// MappingFEField. + + +template +class Test +{ +public: + Test(); + + void + run(); + +private: + // Construct a triangulation with a single cell + void + setup_single_cell_triangulation(); + + // Make the ImmersedSurfaceQuadraturewith contain a single point in the center + // of the cell, with weight .5 and a normal different in all components. + void + setup_single_point_quadrature(); + + // Test that the constant n_quadrature_points corresponds to what we sent in. + void + test_n_quadrature_points(); + + // Test that we can use the get_quadrature function to get the quadrature + // we passed to the FEImmersedSurfaceValues constructor. + void + test_get_quadrature(); + + // Print the quadrature point in real space to make sure that it's mapped to + // real space correctly. + void + test_point_mapped_correctly(); + + // Print the normal in real space to make sure that it's mapped to real space + // correctly. + void + test_normal(); + + // Print JxW to check that the value is correct. + void + test_JxW(); + + // Test that we can call the shape_surface_grad function. + void + test_shape_surface_grad(); + + const FE_Bernstein element; + const FESystem fe_system; + Triangulation triangulation; + DoFHandler dof_handler; + DoFHandler euler_dof_handler; + Vector euler_vector; + std::unique_ptr> euler_mapping; + NonMatching::ImmersedSurfaceQuadrature quadrature; +}; + + + +template +Test::Test() + : element(1) + , fe_system(element, dim) + , dof_handler(triangulation) + , euler_dof_handler(triangulation) +{} + + + +template +void +Test::run() +{ + setup_single_cell_triangulation(); + + setup_single_point_quadrature(); + + test_n_quadrature_points(); + test_get_quadrature(); + test_point_mapped_correctly(); + test_normal(); + test_JxW(); + test_shape_surface_grad(); +} + + + +template +void +Test::setup_single_cell_triangulation() +{ + const Point lower_left; + Point upper_right; + for (unsigned int d = 0; d < dim; ++d) + { + upper_right[d] = 2 + d; + } + + GridGenerator::hyper_rectangle(triangulation, lower_left, upper_right); + dof_handler.distribute_dofs(element); + euler_dof_handler.distribute_dofs(fe_system); + euler_vector.reinit(euler_dof_handler.n_dofs()); + + const ComponentMask mask(dim, true); + VectorTools::get_position_vector(euler_dof_handler, euler_vector, mask); + euler_mapping = std::make_unique>(euler_dof_handler, + euler_vector, + mask); +} + + + +template +void +Test::setup_single_point_quadrature() +{ + Point point; + const double weight = .5; + Tensor<1, dim> normal; + for (unsigned int d = 0; d < dim; ++d) + { + point[d] = .5; + normal[d] = d + 1; + } + normal /= normal.norm(); + quadrature.push_back(point, weight, normal); +} + + + +template +void +Test::test_n_quadrature_points() +{ + FEImmersedSurfaceValues fe_values(*euler_mapping, + element, + quadrature, + update_default); + fe_values.reinit(triangulation.begin_active()); + + deallog << "n_quadrature_points = " << fe_values.n_quadrature_points + << std::endl; +} + + + +template +void +Test::test_get_quadrature() +{ + FEImmersedSurfaceValues fe_values(*euler_mapping, + element, + quadrature, + update_default); + fe_values.reinit(triangulation.begin_active()); + + const NonMatching::ImmersedSurfaceQuadrature &stored_quadrature = + fe_values.get_quadrature(); + + for (unsigned int q = 0; q < stored_quadrature.size(); q++) + { + deallog << "(point, weight, normal) = ([" << stored_quadrature.point(q) + << "], " << stored_quadrature.weight(q) << ", [" + << stored_quadrature.normal_vector(q) << "])" << std::endl; + } +} + + + +template +void +Test::test_point_mapped_correctly() +{ + FEImmersedSurfaceValues fe_values(*euler_mapping, + element, + quadrature, + update_quadrature_points); + fe_values.reinit(triangulation.begin_active()); + + deallog << "point = " << fe_values.quadrature_point(0) << std::endl; +} + + + +template +void +Test::test_normal() +{ + FEImmersedSurfaceValues fe_values(*euler_mapping, + element, + quadrature, + update_normal_vectors); + fe_values.reinit(triangulation.begin_active()); + + deallog << "normal = " << fe_values.normal_vector(0) << std::endl; +} + + + +template +void +Test::test_JxW() +{ + FEImmersedSurfaceValues fe_values(*euler_mapping, + element, + quadrature, + update_JxW_values); + fe_values.reinit(triangulation.begin_active()); + + deallog << "JxW = " << fe_values.JxW(0) << std::endl; +} + + + +template +void +Test::test_shape_surface_grad() +{ + FEImmersedSurfaceValues fe_values(*euler_mapping, + element, + quadrature, + update_gradients | + update_normal_vectors); + fe_values.reinit(dof_handler.begin_active()); + + const unsigned int function_index = 0; + const unsigned int q_index = 0; + const unsigned int component = 0; + + deallog << "shape_surface_grad = " + << fe_values.shape_surface_grad(function_index, q_index) << std::endl; + + deallog << "shape_surface_grad_component = " + << fe_values.shape_surface_grad_component(function_index, + q_index, + component) + << std::endl; +} + + + +template +void +run_test() +{ + deallog << "dim = " << dim << std::endl; + Test test; + test.run(); + deallog << std::endl; +} + + + +int +main() +{ + initlog(); + run_test<1>(); + run_test<2>(); + run_test<3>(); +} diff --git a/tests/non_matching/fe_immersed_surface_values_02.output b/tests/non_matching/fe_immersed_surface_values_02.output new file mode 100644 index 0000000000..4d4f0d8178 --- /dev/null +++ b/tests/non_matching/fe_immersed_surface_values_02.output @@ -0,0 +1,28 @@ + +DEAL::dim = 1 +DEAL::n_quadrature_points = 1 +DEAL::(point, weight, normal) = ([0.500000], 0.500000, [1.00000]) +DEAL::point = 1.00000 +DEAL::normal = 1.00000 +DEAL::JxW = 0.500000 +DEAL::shape_surface_grad = 0.00000 +DEAL::shape_surface_grad_component = 0.00000 +DEAL:: +DEAL::dim = 2 +DEAL::n_quadrature_points = 1 +DEAL::(point, weight, normal) = ([0.500000 0.500000], 0.500000, [0.447214 0.894427]) +DEAL::point = 1.00000 1.50000 +DEAL::normal = 0.600000 0.800000 +DEAL::JxW = 1.11803 +DEAL::shape_surface_grad = -0.0800000 0.0600000 +DEAL::shape_surface_grad_component = -0.0800000 0.0600000 +DEAL:: +DEAL::dim = 3 +DEAL::n_quadrature_points = 1 +DEAL::(point, weight, normal) = ([0.500000 0.500000 0.500000], 0.500000, [0.267261 0.534522 0.801784]) +DEAL::point = 1.00000 1.50000 2.00000 +DEAL::normal = 0.445976 0.594635 0.668965 +DEAL::JxW = 3.59563 +DEAL::shape_surface_grad = -0.0593923 0.00414365 0.0359116 +DEAL::shape_surface_grad_component = -0.0593923 0.00414365 0.0359116 +DEAL::