From: Simon Sticko Date: Wed, 17 Nov 2021 09:42:24 +0000 (+0100) Subject: Add class NonMatching::FEImmersedSurfaceValues X-Git-Tag: v9.4.0-rc1~822^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F12024%2Fhead;p=dealii.git Add class NonMatching::FEImmersedSurfaceValues That has the same shape function values/derivatives as FEValues, but has JxW-values corresponding to integration over a surface immersed in the cell. This is achieved by passing an ImmersedSurfaceQuadrature instead of a Quadrature to the constructor. Add a virtual function fill_fe_immersed_surface_values on Mapping to be able to compute the correct JxW-values and surface normals. Make it throw ExcNotImplemented by default, but override it with an implementation on MappingCartesian and MappingQ. --- diff --git a/include/deal.II/fe/fe.h b/include/deal.II/fe/fe.h index 4536a68747..453e5e519c 100644 --- a/include/deal.II/fe/fe.h +++ b/include/deal.II/fe/fe.h @@ -40,6 +40,11 @@ template class FEFaceValues; template class FESubfaceValues; +namespace NonMatching +{ + template + class FEImmersedSurfaceValues; +} template class FESystem; @@ -3068,6 +3073,7 @@ protected: friend class FEValues; friend class FEFaceValues; friend class FESubfaceValues; + friend class NonMatching::FEImmersedSurfaceValues; friend class FESystem; // explicitly check for sensible template arguments, but not on windows diff --git a/include/deal.II/fe/mapping.h b/include/deal.II/fe/mapping.h index 6c38099873..f5510b58a2 100644 --- a/include/deal.II/fe/mapping.h +++ b/include/deal.II/fe/mapping.h @@ -28,6 +28,8 @@ #include +#include + #include #include #include @@ -48,6 +50,11 @@ template class FEFaceValues; template class FESubfaceValues; +namespace NonMatching +{ + template + class FEImmersedSurfaceValues; +} /** @@ -981,6 +988,20 @@ protected: dealii::internal::FEValuesImplementation::MappingRelatedData &output_data) const = 0; + /** + * The equivalent of Mapping::fill_fe_values(), but for the case that the + * quadrature is an ImmersedSurfaceQuadrature. See there for a comprehensive + * description of the input parameters. This function is called by + * FEImmersedSurfaceValues::reinit(). + */ + 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; + /** * @} */ @@ -1282,6 +1303,7 @@ public: friend class FEValues; friend class FEFaceValues; friend class FESubfaceValues; + friend class NonMatching::FEImmersedSurfaceValues; }; diff --git a/include/deal.II/fe/mapping_cartesian.h b/include/deal.II/fe/mapping_cartesian.h index 127d570f25..e2838213d7 100644 --- a/include/deal.II/fe/mapping_cartesian.h +++ b/include/deal.II/fe/mapping_cartesian.h @@ -287,6 +287,15 @@ private: internal::FEValuesImplementation::MappingRelatedData &output_data) const override; + // documentation can be found in Mapping::fill_fe_immersed_surface_values() + 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/fe/mapping_q.h b/include/deal.II/fe/mapping_q.h index 29eccf2bad..7cf2303251 100644 --- a/include/deal.II/fe/mapping_q.h +++ b/include/deal.II/fe/mapping_q.h @@ -615,6 +615,15 @@ protected: dealii::internal::FEValuesImplementation::MappingRelatedData &output_data) const override; + // documentation can be found in Mapping::fill_fe_immersed_surface_values() + 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 new file mode 100644 index 0000000000..daf0034fdf --- /dev/null +++ b/include/deal.II/non_matching/fe_immersed_values.h @@ -0,0 +1,153 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 - 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +#ifndef dealii_non_matching_fe_immersed_values_h +#define dealii_non_matching_fe_immersed_values_h + +#include + +#include + +#include +#include +#include +#include + +#include + +#include + +DEAL_II_NAMESPACE_OPEN +namespace NonMatching +{ + /** + * Finite element evaluated in the quadrature points of an + * ImmersedSurfaceQuadrature of a cell. + * + * The shape functions values and their derivatives are the same as for an + * FEValues-object, but the JxW-values are computed with the transformation + * described in the documentation of ImmersedSurfaceQuadrature. Further, the + * normal_vector-function returns the normal to the immersed surface. + * + * The reinit-function of this class exist mostly to be consistent with the + * other FEValues-like classes. The immersed quadrature rule will typically + * vary between each cell of the triangulation. Thus, an + * FEImmersedSurfaceValues object can, typically, not be reused for different + * cells. + * + * See also documentation in FEValuesBase. + * + * @ingroup feaccess + */ + template + class FEImmersedSurfaceValues : public FEValuesBase + { + public: + /** + * 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. + */ + FEImmersedSurfaceValues(const Mapping & mapping, + const FiniteElement & element, + const ImmersedSurfaceQuadrature &quadrature, + const UpdateFlags update_flags); + + /** + * Reinitialize quantities (normals, JxW-values, etc) for the given cell of + * type "iterator into a Triangulation object". + */ + void + reinit(const typename Triangulation::cell_iterator &cell); + + /** + * Reinitialize quantitites (shape function values, gradients, etc) for the + * given cell of type "iterator into a DoFHandler object", and the finite + * element associated with this object. + */ + template + void + reinit( + const TriaIterator> &cell); + + /** + * Returns the surface gradient of the shape function with index + * @p function_no at the quadrature point with index @p quadrature_point. + * + * The surface gradient is defined as the projection of the gradient to the + * tangent plane of the surface: + * $ \nabla u - (n \cdot \nabla u) n $, + * where $n$ is the unit normal to the surface. + * + * @dealiiRequiresUpdateFlags{update_gradients | update_normal_vectors} + */ + Tensor<1, dim> + shape_surface_grad(const unsigned int function_no, + const unsigned int quadrature_point) const; + + /** + * Return one vector component of the surface gradient of the shape function + * at a quadrature point. See the definition of the surface gradient in the + * shape_surface_grad function. + * + * @p function_no Index of the shape function to be evaluated. + * + * @p point_no Index of the quadrature point at which the function is to be + * evaluated. + * + * @p component Vector component to be evaluated. + * + * @dealiiRequiresUpdateFlags{update_gradients | update_normal_vectors} + */ + Tensor<1, dim> + shape_surface_grad_component(const unsigned int function_no, + const unsigned int quadrature_point, + const unsigned int component) const; + + /** + * Return a reference to the copy of the quadrature formula stored by this + * object. + */ + const NonMatching::ImmersedSurfaceQuadrature & + get_quadrature() const; + + protected: + /** + * Do work common to the constructors. + */ + void + initialize(const UpdateFlags update_flags); + + /** + * The reinit() functions do only that part of the work that requires + * knowledge of the type of iterator. After setting present_cell(), they + * pass on to this function, which does the real work, and which is + * independent of the actual type of the cell iterator. + */ + void + do_reinit(); + + /** + * Copy of the quadrature formula that was passed to the constructor. + */ + const ImmersedSurfaceQuadrature quadrature; + }; + +} // namespace NonMatching +DEAL_II_NAMESPACE_CLOSE + +#endif /* dealii_non_matching_fe_immersed_values_h */ diff --git a/source/fe/mapping.cc b/source/fe/mapping.cc index b075dd6b65..0aeb35b77e 100644 --- a/source/fe/mapping.cc +++ b/source/fe/mapping.cc @@ -82,6 +82,20 @@ Mapping::get_bounding_box( +template +void +Mapping::fill_fe_immersed_surface_values( + const typename Triangulation::cell_iterator &, + const NonMatching::ImmersedSurfaceQuadrature &, + const typename Mapping::InternalDataBase &, + dealii::internal::FEValuesImplementation::MappingRelatedData &) + const +{ + AssertThrow(false, ExcNotImplemented()); +} + + + template void Mapping::transform_points_real_to_unit_cell( diff --git a/source/fe/mapping_cartesian.cc b/source/fe/mapping_cartesian.cc index 2d29cbc65b..84cc8a3ddc 100644 --- a/source/fe/mapping_cartesian.cc +++ b/source/fe/mapping_cartesian.cc @@ -637,6 +637,69 @@ MappingCartesian::fill_fe_subface_values( +template +void +MappingCartesian::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); + + // Convert data object to internal data for this class. Fails with an + // exception if that is not possible. + Assert(dynamic_cast(&internal_data) != nullptr, + ExcInternalError()); + const InternalData &data = static_cast(internal_data); + + + update_cell_extents(cell, CellSimilarity::none, data); + + maybe_update_cell_quadrature_points(cell, + data, + output_data.quadrature_points); + + if (data.update_each & update_normal_vectors) + for (unsigned int i = 0; i < output_data.normal_vectors.size(); ++i) + { + // The normals are n = J^{-T} * \hat{n} before normalizing. + Tensor<1, dim> normal; + const Tensor<1, dim> &ref_space_normal = quadrature.normal_vector(i); + for (unsigned int d = 0; d < dim; ++d) + { + normal[d] = ref_space_normal[d] / data.cell_extents[d]; + } + normal /= normal.norm(); + output_data.normal_vectors[i] = normal; + } + + if (data.update_each & update_JxW_values) + for (unsigned int i = 0; i < output_data.JxW_values.size(); ++i) + { + const Tensor<1, dim> &ref_space_normal = quadrature.normal_vector(i); + + // J^{-T} \times \hat{n} + Tensor<1, dim> invJTxNormal; + double det_jacobian = 1.; + for (unsigned int d = 0; d < dim; ++d) + { + det_jacobian *= data.cell_extents[d]; + invJTxNormal[d] = ref_space_normal[d] / data.cell_extents[d]; + } + output_data.JxW_values[i] = + det_jacobian * invJTxNormal.norm() * quadrature.weight(i); + } + + maybe_update_volume_elements(data); + maybe_update_jacobians(data, CellSimilarity::none, output_data); + maybe_update_jacobian_derivatives(data, CellSimilarity::none, output_data); + maybe_update_inverse_jacobians(data, CellSimilarity::none, output_data); +} + + + template void MappingCartesian::transform( diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc index 2f4c30aeb4..ed4bfe9d0e 100644 --- a/source/fe/mapping_q.cc +++ b/source/fe/mapping_q.cc @@ -1259,6 +1259,139 @@ MappingQ::fill_fe_subface_values( +template +void +MappingQ::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 +{ + Assert(dim == spacedim, ExcNotImplemented()); + + // ensure that the following static_cast is really correct: + Assert(dynamic_cast(&internal_data) != nullptr, + ExcInternalError()); + const InternalData &data = static_cast(internal_data); + + const unsigned int n_q_points = quadrature.size(); + + data.mapping_support_points = this->compute_mapping_support_points(cell); + data.cell_of_current_support_points = cell; + + internal::MappingQImplementation::maybe_compute_q_points( + QProjector::DataSetDescriptor::cell(), + data, + output_data.quadrature_points); + + internal::MappingQImplementation::maybe_update_Jacobians( + CellSimilarity::none, QProjector::DataSetDescriptor::cell(), data); + + internal::MappingQImplementation::maybe_update_jacobian_grads( + CellSimilarity::none, + QProjector::DataSetDescriptor::cell(), + data, + output_data.jacobian_grads); + + internal::MappingQImplementation::maybe_update_jacobian_pushed_forward_grads< + dim, + spacedim>(CellSimilarity::none, + QProjector::DataSetDescriptor::cell(), + data, + output_data.jacobian_pushed_forward_grads); + + internal::MappingQImplementation::maybe_update_jacobian_2nd_derivatives< + dim, + spacedim>(CellSimilarity::none, + QProjector::DataSetDescriptor::cell(), + data, + output_data.jacobian_2nd_derivatives); + + internal::MappingQImplementation:: + maybe_update_jacobian_pushed_forward_2nd_derivatives( + CellSimilarity::none, + QProjector::DataSetDescriptor::cell(), + data, + output_data.jacobian_pushed_forward_2nd_derivatives); + + internal::MappingQImplementation::maybe_update_jacobian_3rd_derivatives< + dim, + spacedim>(CellSimilarity::none, + QProjector::DataSetDescriptor::cell(), + data, + output_data.jacobian_3rd_derivatives); + + internal::MappingQImplementation:: + maybe_update_jacobian_pushed_forward_3rd_derivatives( + CellSimilarity::none, + QProjector::DataSetDescriptor::cell(), + data, + output_data.jacobian_pushed_forward_3rd_derivatives); + + 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) + { + 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(); + } +} + + + template void MappingQ::fill_mapping_data_for_generic_points( diff --git a/source/non_matching/CMakeLists.txt b/source/non_matching/CMakeLists.txt index aec821f02a..aa3d565b8b 100644 --- a/source/non_matching/CMakeLists.txt +++ b/source/non_matching/CMakeLists.txt @@ -16,13 +16,15 @@ INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR}) SET(_src + fe_immersed_values.cc mesh_classifier.cc quadrature_generator.cc coupling.cc immersed_surface_quadrature.cc ) -SET(_inst + SET(_inst + fe_immersed_values.inst.in mesh_classifier.inst.in quadrature_generator.inst.in coupling.inst.in diff --git a/source/non_matching/fe_immersed_values.cc b/source/non_matching/fe_immersed_values.cc new file mode 100644 index 0000000000..a844126daf --- /dev/null +++ b/source/non_matching/fe_immersed_values.cc @@ -0,0 +1,216 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 - 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#include + +#include + + +DEAL_II_NAMESPACE_OPEN +namespace NonMatching +{ + template + FEImmersedSurfaceValues::FEImmersedSurfaceValues( + const Mapping & mapping, + const FiniteElement & element, + const ImmersedSurfaceQuadrature &quadrature, + const UpdateFlags update_flags) + : FEValuesBase(quadrature.size(), + element.dofs_per_cell, + update_default, + mapping, + element) + , quadrature(quadrature) + { + initialize(update_flags); + } + + + + template + void + FEImmersedSurfaceValues::reinit( + const typename Triangulation::cell_iterator &cell) + { + // Check that mapping and reference cell type are compatible: + Assert(this->get_mapping().is_compatible_with(cell->reference_cell()), + ExcMessage( + "You are trying to call FEValues::reinit() with a cell of type " + + cell->reference_cell().to_string() + + " with a Mapping that is not compatible with it.")); + + // No FE in this cell, so no assertion necessary here. + this->maybe_invalidate_previous_present_cell(cell); + this->present_cell = {cell}; + + // This was the part of the work that is dependent on the actual data type + // of the iterator. Now pass on to the function doing the real work. + do_reinit(); + } + + + + template + template + void + FEImmersedSurfaceValues::reinit( + const TriaIterator> &cell) + { + // Check that mapping and reference cell type are compatible: + Assert(this->get_mapping().is_compatible_with(cell->reference_cell()), + ExcMessage( + "You are trying to call FEValues::reinit() with a cell of type " + + cell->reference_cell().to_string() + + " with a Mapping that is not compatible with it.")); + + // Assert that the finite elements passed to the constructor and used by the + // DoFHandler used by this cell, are the same + Assert(static_cast &>(*this->fe) == + static_cast &>(cell->get_fe()), + (typename FEValuesBase::ExcFEDontMatch())); + + this->maybe_invalidate_previous_present_cell(cell); + this->present_cell = {cell}; + + // This was the part of the work that is dependent on the actual data type + // of the iterator. Now pass on to the function doing the real work. + do_reinit(); + } + + + + template + void + FEImmersedSurfaceValues::do_reinit() + { + // First call the mapping and let it generate the data specific to the + // mapping. + if (this->update_flags & update_mapping) + { + this->get_mapping().fill_fe_immersed_surface_values( + this->present_cell, + quadrature, + *this->mapping_data, + this->mapping_output); + } + + // Call the finite element and, with the data already filled by the mapping, + // let it compute the data for the mapped shape function values, gradients + // etc. + this->get_fe().fill_fe_values(this->present_cell, + CellSimilarity::none, + this->quadrature, + this->get_mapping(), + *this->mapping_data, + this->mapping_output, + *this->fe_data, + this->finite_element_output); + } + + + + template + Tensor<1, dim> + FEImmersedSurfaceValues::shape_surface_grad( + const unsigned int function_no, + const unsigned int quadrature_point) const + { + const unsigned int component = 0; + return shape_surface_grad_component(function_no, + quadrature_point, + component); + } + + + + template + Tensor<1, dim> + FEImmersedSurfaceValues::shape_surface_grad_component( + const unsigned int function_no, + const unsigned int quadrature_point, + const unsigned int component) const + { + const Tensor<1, dim> gradient = + this->shape_grad_component(function_no, quadrature_point, component); + const Tensor<1, dim> &normal = this->normal_vector(quadrature_point); + + return gradient - (normal * gradient) * normal; + } + + + + template + const ImmersedSurfaceQuadrature & + FEImmersedSurfaceValues::get_quadrature() const + { + return quadrature; + } + + + + template + inline void + FEImmersedSurfaceValues::initialize(const UpdateFlags update_flags) + { + UpdateFlags flags = this->compute_update_flags(update_flags); + + if (flags & (update_JxW_values | update_normal_vectors)) + flags |= update_covariant_transformation; + + // Initialize the base classes. + if (flags & update_mapping) + this->mapping_output.initialize(this->n_quadrature_points, flags); + this->finite_element_output.initialize(this->n_quadrature_points, + *this->fe, + flags); + + // Then get objects into which the FE and the Mapping can store + // intermediate data used across calls to reinit. We can do this in + // parallel. + Threads::Task< + std::unique_ptr::InternalDataBase>> + fe_get_data = Threads::new_task(&FiniteElement::get_data, + *this->fe, + flags, + *this->mapping, + this->quadrature, + this->finite_element_output); + + Threads::Task::InternalDataBase>> + mapping_get_data; + if (flags & update_mapping) + mapping_get_data = Threads::new_task(&Mapping::get_data, + *this->mapping, + flags, + this->quadrature); + + this->update_flags = flags; + + // Then collect answers from the two task above. + this->fe_data = std::move(fe_get_data.return_value()); + if (flags & update_mapping) + this->mapping_data = std::move(mapping_get_data.return_value()); + else + this->mapping_data = + std::make_unique::InternalDataBase>(); + } + + +#ifndef DOXYGEN +# include "fe_immersed_values.inst" +#endif + +} // namespace NonMatching +DEAL_II_NAMESPACE_CLOSE diff --git a/source/non_matching/fe_immersed_values.inst.in b/source/non_matching/fe_immersed_values.inst.in new file mode 100644 index 0000000000..639964cd86 --- /dev/null +++ b/source/non_matching/fe_immersed_values.inst.in @@ -0,0 +1,26 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 - 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +for (deal_II_dimension : DIMENSIONS) + { + template class FEImmersedSurfaceValues; + } + +for (deal_II_dimension : DIMENSIONS; lda : BOOL) + { + template void FEImmersedSurfaceValues::reinit( + const TriaIterator< + DoFCellAccessor> &); + } diff --git a/tests/non_matching/fe_immersed_surface_values.cc b/tests/non_matching/fe_immersed_surface_values.cc new file mode 100644 index 0000000000..9e67f21e82 --- /dev/null +++ b/tests/non_matching/fe_immersed_surface_values.cc @@ -0,0 +1,305 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 - 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +#include +#include + +#include + +#include +#include +#include + +#include +#include + +#include + +#include +#include + +#include "../tests.h" + + +using namespace dealii; +using NonMatching::FEImmersedSurfaceValues; + +// Set up a triangulation with a single Cartesian cell and a finite element +// space with FE_Q-elements. Test that we can construct an +// FEImmersedSurfaceValues object and compute known values when using a +// Cartesian Mapping. +template +class Test +{ +public: + Test(const Mapping &mapping); + + void + run(); + +private: + // Construct a triangulation with a single cartesian cell, with extents. + // [0,2] in 1D, + // [0,2]x[0,3] in 2D, and + // [0,2]x[0,3]x[0,4] in 3D. + // Use extends that are different in all directions so that it's easier to + // check that the various quantities are mapped correctly. + void + setup_single_cartesian_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_Q element; + Triangulation triangulation; + DoFHandler dof_handler; + const SmartPointer> mapping; + NonMatching::ImmersedSurfaceQuadrature quadrature; +}; + + + +template +Test::Test(const Mapping &mapping) + : element(1) + , dof_handler(triangulation) + , mapping(&mapping) +{} + + + +template +void +Test::run() +{ + setup_single_cartesian_cell_triangulation(); + dof_handler.distribute_dofs(element); + 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_cartesian_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); +} + + + +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(*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(*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(*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(*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(*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( + *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; + + const MappingCartesian mapping_cartesian; + const unsigned int polynomial_degree = 1; + const MappingQ mapping_q(polynomial_degree); + + const std::vector mapping_names = {"MappingCartesian", + "MappingQ"}; + const hp::MappingCollection mappings(mapping_cartesian, mapping_q); + + for (unsigned int i = 0; i < mappings.size(); i++) + { + deallog << mapping_names.at(i) << std::endl; + Test test(mappings[i]); + 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.output b/tests/non_matching/fe_immersed_surface_values.output new file mode 100644 index 0000000000..797c582463 --- /dev/null +++ b/tests/non_matching/fe_immersed_surface_values.output @@ -0,0 +1,58 @@ + +DEAL::dim = 1 +DEAL::MappingCartesian +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::MappingQ +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::MappingCartesian +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::MappingQ +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::MappingCartesian +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:: +DEAL::MappingQ +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::