From 9a3edfead47b355eee7f45d8e5aaca6965240dc9 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Sat, 6 Apr 2024 12:12:34 +0200 Subject: [PATCH] More efficient access via new Mapping::InternalDataBase::reinit --- include/deal.II/fe/mapping.h | 17 +- include/deal.II/fe/mapping_cartesian.h | 25 +-- include/deal.II/fe/mapping_fe.h | 18 +- include/deal.II/fe/mapping_fe_field.h | 43 ++-- include/deal.II/fe/mapping_manifold.h | 47 +---- include/deal.II/fe/mapping_q.h | 12 +- include/deal.II/non_matching/mapping_info.h | 68 +++--- source/fe/fe_hermite.cc | 14 +- source/fe/mapping.cc | 10 + source/fe/mapping_cartesian.cc | 98 +++++---- source/fe/mapping_fe.cc | 19 +- source/fe/mapping_fe_field.cc | 222 ++++++++++---------- source/fe/mapping_manifold.cc | 65 +++++- source/fe/mapping_q.cc | 14 +- 14 files changed, 331 insertions(+), 341 deletions(-) diff --git a/include/deal.II/fe/mapping.h b/include/deal.II/fe/mapping.h index 74e0b11267..baba788e1c 100644 --- a/include/deal.II/fe/mapping.h +++ b/include/deal.II/fe/mapping.h @@ -59,6 +59,8 @@ namespace NonMatching template class ComputeMappingDataHelper; } + template + class MappingInfo; } // namespace NonMatching @@ -655,6 +657,18 @@ public: */ InternalDataBase(const InternalDataBase &) = delete; + /** + * This function initializes the data fields related to evaluation of the + * mapping on cells, implemented by (derived) classes. This function is + * used both when setting up a field of this class for the first time or + * when a new Quadrature formula should be considered without creating an + * entirely new object. This is used when the number of evaluation points + * is different on each cell, e.g. when using FEPointEvaluation for + * handling particles or with certain non-matching problem settings. + */ + virtual void + reinit(const UpdateFlags update_flags, const Quadrature &quadrature); + /** * Virtual destructor for derived classes */ @@ -683,7 +697,6 @@ public: memory_consumption() const; }; - protected: /** * Given a set of update flags, compute which other quantities also @@ -1323,6 +1336,8 @@ public: friend class FESubfaceValues; friend class NonMatching::FEImmersedSurfaceValues; friend class NonMatching::internal::ComputeMappingDataHelper; + template + friend class NonMatching::MappingInfo; }; diff --git a/include/deal.II/fe/mapping_cartesian.h b/include/deal.II/fe/mapping_cartesian.h index 73f1f2006f..14dbaf6788 100644 --- a/include/deal.II/fe/mapping_cartesian.h +++ b/include/deal.II/fe/mapping_cartesian.h @@ -218,6 +218,11 @@ public: */ InternalData(const Quadrature &quadrature); + // Documentation see Mapping::InternalDataBase. + virtual void + reinit(const UpdateFlags update_flags, + const Quadrature &quadrature) override; + /** * Return an estimate (in bytes) for the memory consumption of this object. */ @@ -243,7 +248,9 @@ public: mutable double volume_element; /** - * Vector of all quadrature points. Especially, all points on all faces. + * Location of quadrature points of faces or subfaces in 3d with all + * possible orientations. Can be accessed with the correct offset provided + * via QProjector::DataSetDescriptor. Not needed/used for cells. */ std::vector> quadrature_points; }; @@ -335,7 +342,8 @@ private: maybe_update_cell_quadrature_points( const typename Triangulation::cell_iterator &cell, const InternalData &data, - std::vector> &quadrature_points) const; + const ArrayView> &unit_quadrature_points, + std::vector> &quadrature_points) const; /** * Compute the quadrature points if the UpdateFlags of the incoming @@ -364,19 +372,6 @@ private: const InternalData &data, std::vector> &quadrature_points) const; - /** - * Transform quadrature points in InternalData to real space by scaling unit - * coordinates with cell_extents in each direction. - * - * Called from the various maybe_update_*_quadrature_points functions. - */ - void - transform_quadrature_points( - const typename Triangulation::cell_iterator &cell, - const InternalData &data, - const typename QProjector::DataSetDescriptor &offset, - std::vector> &quadrature_points) const; - /** * Compute the normal vectors if the UpdateFlags of the incoming InternalData * object say that they should be updated. diff --git a/include/deal.II/fe/mapping_fe.h b/include/deal.II/fe/mapping_fe.h index ff265d1c55..48c9e18bd3 100644 --- a/include/deal.II/fe/mapping_fe.h +++ b/include/deal.II/fe/mapping_fe.h @@ -185,23 +185,15 @@ public: */ InternalData(const FiniteElement &fe); - /** - * Initialize the object's member variables related to cell data based on - * the given arguments. - * - * The function also calls compute_shape_function_values() to actually set - * the member variables related to the values and derivatives of the - * mapping shape functions. - */ - void - initialize(const UpdateFlags update_flags, - const Quadrature &quadrature, - const unsigned int n_original_q_points); + // Documentation see Mapping::InternalDataBase. + virtual void + reinit(const UpdateFlags update_flags, + const Quadrature &quadrature) override; /** * Initialize the object's member variables related to cell and face data * based on the given arguments. In order to initialize cell data, this - * function calls initialize(). + * function calls reinit(). */ void initialize_face(const UpdateFlags update_flags, diff --git a/include/deal.II/fe/mapping_fe_field.h b/include/deal.II/fe/mapping_fe_field.h index fd842faa55..1f90d82033 100644 --- a/include/deal.II/fe/mapping_fe_field.h +++ b/include/deal.II/fe/mapping_fe_field.h @@ -289,6 +289,11 @@ public: InternalData(const FiniteElement &fe, const ComponentMask &mask); + // Documentation see Mapping::InternalDataBase. + virtual void + reinit(const UpdateFlags update_flags, + const Quadrature &quadrature) override; + /** * Shape function at quadrature point. Shape functions are in tensor * product order, so vertices must be reordered to obtain transformation. @@ -359,6 +364,11 @@ public: virtual std::size_t memory_consumption() const override; + /** + * A pointer to the underlying finite element. + */ + SmartPointer> fe; + /** * Values of shape functions. Access by function @p shape. * @@ -590,7 +600,6 @@ private: Point do_transform_unit_to_real_cell(const InternalData &mdata) const; - /** * Transform the point @p p on the real cell to the corresponding point on * the unit cell @p cell by a Newton iteration. @@ -598,10 +607,8 @@ private: * Takes a reference to an @p InternalData that is assumed to be previously * created by the @p get_data function with @p UpdateFlags including @p * update_transformation_values and @p update_transformation_gradients and a - * one point Quadrature that includes the given initial guess for the - * transformation @p initial_p_unit. Hence this function assumes that @p - * mdata already includes the transformation shape values and gradients - * computed at @p initial_p_unit. + * one point Quadrature that includes the given initial guess specified + * through the given @p point_quadrature object. * * @p mdata will be changed by this function. */ @@ -609,8 +616,8 @@ private: do_transform_real_to_unit_cell( const typename Triangulation::cell_iterator &cell, const Point &p, - const Point &initial_p_unit, - InternalData &mdata) const; + Quadrature &point_quadrature, + InternalData &mdata) const; /** * Update internal degrees of freedom. @@ -621,15 +628,6 @@ private: const typename MappingFEField::InternalData &data) const; - /** - * See the documentation of the base class for detailed information. - */ - virtual void - compute_shapes_virtual( - const std::vector> &unit_points, - typename MappingFEField::InternalData &data) - const; - /* * Which components to use for the mapping. */ @@ -658,17 +656,8 @@ private: mutable Threads::Mutex fe_values_mutex; void - compute_data(const UpdateFlags update_flags, - const Quadrature &q, - const unsigned int n_original_q_points, - InternalData &data) const; - - void - compute_face_data(const UpdateFlags update_flags, - const Quadrature &q, - const unsigned int n_original_q_points, - InternalData &data) const; - + compute_face_data(const unsigned int n_original_q_points, + InternalData &data) const; // Declare other MappingFEField classes friends. template diff --git a/include/deal.II/fe/mapping_manifold.h b/include/deal.II/fe/mapping_manifold.h index db10995ce0..b7de835a49 100644 --- a/include/deal.II/fe/mapping_manifold.h +++ b/include/deal.II/fe/mapping_manifold.h @@ -167,18 +167,10 @@ public: */ InternalData() = default; - /** - * Initialize the object's member variables related to cell data based on - * the given arguments. - * - * The function also calls compute_shape_function_values() to actually set - * the member variables related to the values and derivatives of the - * mapping shape functions. - */ - void - initialize(const UpdateFlags update_flags, - const Quadrature &quadrature, - const unsigned int n_original_q_points); + // Documentation see Mapping::InternalDataBase. + virtual void + reinit(const UpdateFlags update_flags, + const Quadrature &quadrature) override; /** * Initialize the object's member variables related to cell and face data @@ -391,37 +383,6 @@ private: #ifndef DOXYGEN -template -inline void -MappingManifold::InternalData::store_vertices( - const typename Triangulation::cell_iterator &cell) const -{ - vertices.resize(GeometryInfo::vertices_per_cell); - for (const unsigned int i : GeometryInfo::vertex_indices()) - vertices[i] = cell->vertex(i); - this->cell = cell; -} - - -template -inline void -MappingManifold::InternalData:: - compute_manifold_quadrature_weights(const Quadrature &quad) -{ - cell_manifold_quadrature_weights.resize( - quad.size(), std::vector(GeometryInfo::vertices_per_cell)); - for (unsigned int q = 0; q < quad.size(); ++q) - { - for (const unsigned int i : GeometryInfo::vertex_indices()) - { - cell_manifold_quadrature_weights[q][i] = - GeometryInfo::d_linear_shape_function(quad.point(q), i); - } - } -} - - - template inline bool MappingManifold::preserves_vertex_locations() const diff --git a/include/deal.II/fe/mapping_q.h b/include/deal.II/fe/mapping_q.h index e2da64c35c..a440434ae6 100644 --- a/include/deal.II/fe/mapping_q.h +++ b/include/deal.II/fe/mapping_q.h @@ -309,14 +309,10 @@ public: */ InternalData(const unsigned int polynomial_degree); - /** - * Initialize the object's member variables related to cell data based on - * the given arguments. - */ - void - initialize(const UpdateFlags update_flags, - const Quadrature &quadrature, - const unsigned int n_original_q_points); + // Documentation see Mapping::InternalDataBase. + virtual void + reinit(const UpdateFlags update_flags, + const Quadrature &quadrature) override; /** * Initialize the object's member variables related to cell and face data diff --git a/include/deal.II/non_matching/mapping_info.h b/include/deal.II/non_matching/mapping_info.h index e931f59a1a..93672de532 100644 --- a/include/deal.II/non_matching/mapping_info.h +++ b/include/deal.II/non_matching/mapping_info.h @@ -70,24 +70,7 @@ namespace NonMatching MappingData &mapping_data) { mapping_data.initialize(quadrature.size(), update_flags_mapping); - - // reuse internal_mapping_data for MappingQ to avoid memory allocations - if (const MappingQ *mapping_q = - dynamic_cast *>(mapping.get())) - { - (void)mapping_q; - auto &data = - static_cast::InternalData &>( - *internal_mapping_data); - data.initialize(update_flags_mapping, - quadrature, - quadrature.size()); - } - else - { - internal_mapping_data = - mapping->get_data(update_flags_mapping, quadrature); - } + internal_mapping_data->reinit(update_flags_mapping, quadrature); cell_similarity = mapping->fill_fe_values(cell, cell_similarity, @@ -110,21 +93,7 @@ namespace NonMatching { mapping_data.initialize(quadrature.size(), update_flags_mapping); - // reuse internal_mapping_data for MappingQ to avoid memory allocations - if (dynamic_cast *>(&(*mapping))) - { - auto &data = - dynamic_cast::InternalData &>( - *internal_mapping_data); - data.initialize(update_flags_mapping, - quadrature, - quadrature.size()); - } - else - { - internal_mapping_data = - mapping->get_data(update_flags_mapping, quadrature); - } + internal_mapping_data->reinit(update_flags_mapping, quadrature); mapping->fill_fe_immersed_surface_values(cell, quadrature, @@ -241,6 +210,12 @@ namespace NonMatching * template argument, e.g. VectorizedArray) provides both mapping data * and unit points in vectorized format. The Number template parameter of * MappingInfo and FEPointEvaluation has to be identical. + * + * @note This class cannot be copied or copy-constructed. The intention of + * this class is to be available as a single use, e.g. inside + * FEPointEvaluation or as a storage container for the whole mesh. Use a + * shared or unique pointer of this class in case several objects should be + * held. */ template class MappingInfo : public Subscriptor @@ -300,6 +275,17 @@ namespace NonMatching const UpdateFlags update_flags, const AdditionalData additional_data = AdditionalData()); + /** + * Do not allow making copies. + */ + MappingInfo(const MappingInfo &) = delete; + + /** + * Do not allow copy assignment of this class. + */ + MappingInfo & + operator=(const MappingInfo &) = delete; + /** * Compute the mapping information for the incoming cell and unit * points. This overload is needed to resolve ambiguity. @@ -665,7 +651,9 @@ namespace NonMatching Quadrature quadrature; /** - * An array containing the data fields of this class. + * Helper class that temporarily holds the data requested for one cell by + * Mapping::fill_fe_values() before it is filled into appropriate data + * structures for consumption by FEPointEvaluation. */ MappingData mapping_data; @@ -837,15 +825,9 @@ namespace NonMatching internal::ComputeMappingDataHelper::required_update_flags( this->mapping, update_flags_mapping); - // construct internal_mapping_data for MappingQ to be able to reuse it in - // reinit() calls to avoid memory allocations - if (const MappingQ *mapping_q = - dynamic_cast *>(&mapping)) - { - internal_mapping_data = - std::make_unique::InternalData>( - mapping_q->get_degree()); - } + // construct internal_mapping_data for mappings for reuse in reinit() + // calls to avoid frequent memory allocations + internal_mapping_data = mapping.get_data(update_flags, Quadrature()); } diff --git a/source/fe/fe_hermite.cc b/source/fe/fe_hermite.cc index 47687f18ed..5e6d97c3ec 100644 --- a/source/fe/fe_hermite.cc +++ b/source/fe/fe_hermite.cc @@ -247,8 +247,7 @@ namespace internal const typename Mapping<1, spacedim>::InternalDataBase &mapping_data, Table<2, Number> &value_list) { - unsigned int n_q_points; - double cell_extent = 1.0; + double cell_extent = 1.0; // Check mapping_data is associated with a compatible mapping class if (dynamic_cast::InternalData *>( @@ -257,7 +256,6 @@ namespace internal const typename MappingCartesian<1>::InternalData *data = dynamic_cast::InternalData *>( &mapping_data); - n_q_points = data->quadrature_points.size(); cell_extent = data->cell_extents[0]; } else @@ -270,8 +268,6 @@ namespace internal AssertDimension(value_list.size(0), n_dofs_per_cell); AssertDimension(n_dofs_per_cell, 2 * regularity + 2); - AssertIndexRange(n_q_points_out, n_q_points + 1); - (void)n_q_points; std::vector l2h = hermite_lexicographic_to_hierarchic_numbering<1>(regularity); @@ -306,7 +302,6 @@ namespace internal const typename Mapping<2, spacedim>::InternalDataBase &mapping_data, Table<2, Number> &value_list) { - unsigned int n_q_points; Tensor<1, 2> cell_extents; // Check mapping_data is associated with a compatible mapping class @@ -316,7 +311,6 @@ namespace internal const typename MappingCartesian<2>::InternalData *data = dynamic_cast::InternalData *>( &mapping_data); - n_q_points = data->quadrature_points.size(); cell_extents = data->cell_extents; } else @@ -330,8 +324,6 @@ namespace internal AssertDimension(value_list.size(0), n_dofs_per_cell); AssertDimension(n_dofs_per_dim * n_dofs_per_dim, n_dofs_per_cell); - AssertIndexRange(n_q_points_out, n_q_points + 1); - (void)n_q_points; std::vector l2h = hermite_lexicographic_to_hierarchic_numbering<2>(regularity); @@ -380,7 +372,6 @@ namespace internal const typename Mapping<3, spacedim>::InternalDataBase &mapping_data, Table<2, Number> &value_list) { - unsigned int n_q_points; Tensor<1, 3> cell_extents; // Check mapping_data is associated with a compatible mapping class @@ -390,7 +381,6 @@ namespace internal const typename MappingCartesian<3>::InternalData *data = dynamic_cast::InternalData *>( &mapping_data); - n_q_points = data->quadrature_points.size(); cell_extents = data->cell_extents; } else @@ -405,8 +395,6 @@ namespace internal AssertDimension(value_list.size(0), n_dofs_per_cell); AssertDimension(Utilities::pow(n_dofs_per_dim, 3), n_dofs_per_cell); - AssertIndexRange(n_q_points_out, n_q_points + 1); - (void)n_q_points; std::vector l2h = hermite_lexicographic_to_hierarchic_numbering<3>(regularity); diff --git a/source/fe/mapping.cc b/source/fe/mapping.cc index 4487a54854..314a2d18a8 100644 --- a/source/fe/mapping.cc +++ b/source/fe/mapping.cc @@ -269,6 +269,16 @@ Mapping::InternalDataBase::InternalDataBase() +template +void +Mapping::InternalDataBase::reinit(const UpdateFlags, + const Quadrature &) +{ + DEAL_II_ASSERT_UNREACHABLE(); +} + + + template std::size_t Mapping::InternalDataBase::memory_consumption() const diff --git a/source/fe/mapping_cartesian.cc b/source/fe/mapping_cartesian.cc index 07b185eddf..2bb5977b7e 100644 --- a/source/fe/mapping_cartesian.cc +++ b/source/fe/mapping_cartesian.cc @@ -98,6 +98,20 @@ MappingCartesian::InternalData::InternalData( +template +void +MappingCartesian::InternalData::reinit( + const UpdateFlags update_flags, + const Quadrature &) +{ + // store the flags in the internal data object so we can access them + // in fill_fe_*_values(). use the transitive hull of the required + // flags + this->update_each = update_flags; +} + + + template std::size_t MappingCartesian::InternalData::memory_consumption() const @@ -105,8 +119,7 @@ MappingCartesian::InternalData::memory_consumption() const return (Mapping::InternalDataBase::memory_consumption() + MemoryConsumption::memory_consumption(cell_extents) + MemoryConsumption::memory_consumption(inverse_cell_extents) + - MemoryConsumption::memory_consumption(volume_element) + - MemoryConsumption::memory_consumption(quadrature_points)); + MemoryConsumption::memory_consumption(volume_element)); } @@ -162,13 +175,8 @@ MappingCartesian::get_data(const UpdateFlags update_flags, const Quadrature &q) const { std::unique_ptr::InternalDataBase> data_ptr = - std::make_unique(q); - auto &data = dynamic_cast(*data_ptr); - - // store the flags in the internal data object so we can access them - // in fill_fe_*_values(). use the transitive hull of the required - // flags - data.update_each = requires_update_flags(update_flags); + std::make_unique(); + data_ptr->reinit(requires_update_flags(update_flags), q); return data_ptr; } @@ -252,18 +260,46 @@ MappingCartesian::update_cell_extents( +namespace +{ + template + void + transform_quadrature_points( + const Tensor<1, dim> first_vertex, + const Tensor<1, dim> cell_extents, + const ArrayView> &unit_quadrature_points, + const typename QProjector::DataSetDescriptor &offset, + std::vector> &quadrature_points) + { + for (unsigned int i = 0; i < quadrature_points.size(); ++i) + { + quadrature_points[i] = first_vertex; + for (unsigned int d = 0; d < dim; ++d) + quadrature_points[i][d] += + cell_extents[d] * unit_quadrature_points[i + offset][d]; + } + } +} // namespace + + + template void MappingCartesian::maybe_update_cell_quadrature_points( const typename Triangulation::cell_iterator &cell, const InternalData &data, - std::vector> &quadrature_points) const + const ArrayView> &unit_quadrature_points, + std::vector> &quadrature_points) const { if (data.update_each & update_quadrature_points) { const auto offset = QProjector::DataSetDescriptor::cell(); - transform_quadrature_points(cell, data, offset, quadrature_points); + transform_quadrature_points(cell->vertex(0), + data.cell_extents, + unit_quadrature_points, + offset, + quadrature_points); } } @@ -288,7 +324,11 @@ MappingCartesian::maybe_update_face_quadrature_points( quadrature_points.size()); - transform_quadrature_points(cell, data, offset, quadrature_points); + transform_quadrature_points(cell->vertex(0), + data.cell_extents, + make_array_view(data.quadrature_points), + offset, + quadrature_points); } } @@ -320,30 +360,11 @@ MappingCartesian::maybe_update_subface_quadrature_points( quadrature_points.size(), cell->subface_case(face_no)); - transform_quadrature_points(cell, data, offset, quadrature_points); - } -} - - - -template -void -MappingCartesian::transform_quadrature_points( - const typename Triangulation::cell_iterator &cell, - const InternalData &data, - const typename QProjector::DataSetDescriptor &offset, - std::vector> &quadrature_points) const -{ - // let @p{start} be the origin of a local coordinate system. it is chosen - // as the (lower) left vertex - const Point start = cell->vertex(0); - - for (unsigned int i = 0; i < quadrature_points.size(); ++i) - { - quadrature_points[i] = start; - for (unsigned int d = 0; d < dim; ++d) - quadrature_points[i][d] += - data.cell_extents[d] * data.quadrature_points[i + offset][d]; + transform_quadrature_points(cell->vertex(0), + data.cell_extents, + make_array_view(data.quadrature_points), + offset, + quadrature_points); } } @@ -504,6 +525,7 @@ MappingCartesian::fill_fe_values( maybe_update_cell_quadrature_points(cell, data, + quadrature.get_points(), output_data.quadrature_points); // compute Jacobian determinant. all values are equal and are the @@ -553,13 +575,12 @@ MappingCartesian::fill_mapping_data_for_generic_points( InternalData data; data.update_each = update_flags; - data.quadrature_points = - std::vector>(unit_points.begin(), unit_points.end()); update_cell_extents(cell, CellSimilarity::none, data); maybe_update_cell_quadrature_points(cell, data, + unit_points, output_data.quadrature_points); maybe_update_jacobians(data, CellSimilarity::none, output_data); @@ -703,6 +724,7 @@ MappingCartesian::fill_fe_immersed_surface_values( maybe_update_cell_quadrature_points(cell, data, + quadrature.get_points(), output_data.quadrature_points); if (data.update_each & update_normal_vectors) diff --git a/source/fe/mapping_fe.cc b/source/fe/mapping_fe.cc index 9707bb7a81..7856db406f 100644 --- a/source/fe/mapping_fe.cc +++ b/source/fe/mapping_fe.cc @@ -74,12 +74,11 @@ MappingFE::InternalData::memory_consumption() const } + template void -MappingFE::InternalData::initialize( - const UpdateFlags update_flags, - const Quadrature &q, - const unsigned int n_original_q_points) +MappingFE::InternalData::reinit(const UpdateFlags update_flags, + const Quadrature &q) { // store the flags in the internal data object so we can access them // in fill_fe_*_values() @@ -88,13 +87,13 @@ MappingFE::InternalData::initialize( const unsigned int n_q_points = q.size(); if (this->update_each & update_covariant_transformation) - covariant.resize(n_original_q_points); + covariant.resize(n_q_points); if (this->update_each & update_contravariant_transformation) - contravariant.resize(n_original_q_points); + contravariant.resize(n_q_points); if (this->update_each & update_volume_elements) - volume_elements.resize(n_original_q_points); + volume_elements.resize(n_q_points); // see if we need the (transformation) shape function values // and/or gradients and resize the necessary arrays @@ -139,7 +138,7 @@ MappingFE::InternalData::initialize_face( const Quadrature &q, const unsigned int n_original_q_points) { - initialize(update_flags, q, n_original_q_points); + reinit(update_flags, q); if (this->update_each & (update_boundary_forms | update_normal_vectors | update_jacobians | @@ -1059,9 +1058,7 @@ MappingFE::get_data(const UpdateFlags update_flags, { std::unique_ptr::InternalDataBase> data_ptr = std::make_unique(*this->fe); - auto &data = dynamic_cast(*data_ptr); - data.initialize(this->requires_update_flags(update_flags), q, q.size()); - + data_ptr->reinit(this->requires_update_flags(update_flags), q); return data_ptr; } diff --git a/source/fe/mapping_fe_field.cc b/source/fe/mapping_fe_field.cc index 5fa44c8548..50fd675c30 100644 --- a/source/fe/mapping_fe_field.cc +++ b/source/fe/mapping_fe_field.cc @@ -59,7 +59,8 @@ template MappingFEField::InternalData::InternalData( const FiniteElement &fe, const ComponentMask &mask) - : unit_tangentials() + : fe(&fe) + , unit_tangentials() , n_shape_functions(fe.n_dofs_per_cell()) , mask(mask) , local_dof_indices(fe.n_dofs_per_cell()) @@ -68,6 +69,86 @@ MappingFEField::InternalData::InternalData( +template +void +MappingFEField::InternalData::reinit( + const UpdateFlags update_flags, + const Quadrature &quadrature) +{ + // store the flags in the internal data object so we can access them + // in fill_fe_*_values(). use the transitive hull of the required + // flags + this->update_each = update_flags; + + const unsigned int n_q_points = quadrature.size(); + const std::vector> &points = quadrature.get_points(); + + // see if we need the (transformation) shape function values + // and/or gradients and resize the necessary arrays + if (update_flags & update_quadrature_points) + { + shape_values.resize(n_shape_functions * n_q_points); + for (unsigned int point = 0; point < n_q_points; ++point) + for (unsigned int i = 0; i < n_shape_functions; ++i) + shape(point, i) = fe->shape_value(i, points[point]); + } + + if (update_flags & + (update_covariant_transformation | update_contravariant_transformation | + update_JxW_values | update_boundary_forms | update_normal_vectors | + update_jacobians | update_jacobian_grads | update_inverse_jacobians)) + { + shape_derivatives.resize(n_shape_functions * n_q_points); + for (unsigned int point = 0; point < n_q_points; ++point) + for (unsigned int i = 0; i < n_shape_functions; ++i) + derivative(point, i) = fe->shape_grad(i, points[point]); + } + + if (update_flags & update_covariant_transformation) + covariant.resize(n_q_points); + + if (update_flags & update_contravariant_transformation) + contravariant.resize(n_q_points); + + if (update_flags & update_volume_elements) + volume_elements.resize(n_q_points); + + if (update_flags & + (update_jacobian_grads | update_jacobian_pushed_forward_grads)) + { + shape_second_derivatives.resize(n_shape_functions * n_q_points); + for (unsigned int point = 0; point < n_q_points; ++point) + for (unsigned int i = 0; i < n_shape_functions; ++i) + second_derivative(point, i) = fe->shape_grad_grad(i, points[point]); + } + + if (update_flags & (update_jacobian_2nd_derivatives | + update_jacobian_pushed_forward_2nd_derivatives)) + { + shape_third_derivatives.resize(n_shape_functions * n_q_points); + for (unsigned int point = 0; point < n_q_points; ++point) + for (unsigned int i = 0; i < n_shape_functions; ++i) + third_derivative(point, i) = + fe->shape_3rd_derivative(i, points[point]); + } + + if (update_flags & (update_jacobian_3rd_derivatives | + update_jacobian_pushed_forward_3rd_derivatives)) + { + shape_fourth_derivatives.resize(n_shape_functions * n_q_points); + for (unsigned int point = 0; point < n_q_points; ++point) + for (unsigned int i = 0; i < n_shape_functions; ++i) + fourth_derivative(point, i) = + fe->shape_4th_derivative(i, points[point]); + } + + // This (for face values and simplices) can be different for different + // calls, so always copy + quadrature_weights = quadrature.get_weights(); +} + + + template std::size_t MappingFEField::InternalData::memory_consumption() @@ -420,43 +501,6 @@ MappingFEField::get_vertices( -template -void -MappingFEField::compute_shapes_virtual( - const std::vector> &unit_points, - typename MappingFEField::InternalData &data) const -{ - const auto fe = &euler_dof_handler->get_fe(); - const unsigned int n_points = unit_points.size(); - - for (unsigned int point = 0; point < n_points; ++point) - { - if (data.shape_values.size() != 0) - for (unsigned int i = 0; i < data.n_shape_functions; ++i) - data.shape(point, i) = fe->shape_value(i, unit_points[point]); - - if (data.shape_derivatives.size() != 0) - for (unsigned int i = 0; i < data.n_shape_functions; ++i) - data.derivative(point, i) = fe->shape_grad(i, unit_points[point]); - - if (data.shape_second_derivatives.size() != 0) - for (unsigned int i = 0; i < data.n_shape_functions; ++i) - data.second_derivative(point, i) = - fe->shape_grad_grad(i, unit_points[point]); - - if (data.shape_third_derivatives.size() != 0) - for (unsigned int i = 0; i < data.n_shape_functions; ++i) - data.third_derivative(point, i) = - fe->shape_3rd_derivative(i, unit_points[point]); - - if (data.shape_fourth_derivatives.size() != 0) - for (unsigned int i = 0; i < data.n_shape_functions; ++i) - data.fourth_derivative(point, i) = - fe->shape_4th_derivative(i, unit_points[point]); - } -} - - template UpdateFlags MappingFEField::requires_update_flags( @@ -511,33 +555,14 @@ MappingFEField::requires_update_flags( } - template void -MappingFEField::compute_data( - const UpdateFlags update_flags, - const Quadrature &q, - const unsigned int n_original_q_points, - InternalData &data) const +MappingFEField::compute_face_data( + const unsigned int n_original_q_points, + InternalData &data) const { - // store the flags in the internal data object so we can access them - // in fill_fe_*_values(). use the transitive hull of the required - // flags - data.update_each = requires_update_flags(update_flags); - - const unsigned int n_q_points = q.size(); - - // see if we need the (transformation) shape function values - // and/or gradients and resize the necessary arrays - if (data.update_each & update_quadrature_points) - data.shape_values.resize(data.n_shape_functions * n_q_points); - - if (data.update_each & - (update_covariant_transformation | update_contravariant_transformation | - update_JxW_values | update_boundary_forms | update_normal_vectors | - update_jacobians | update_jacobian_grads | update_inverse_jacobians)) - data.shape_derivatives.resize(data.n_shape_functions * n_q_points); - + // Set to the size of a single quadrature object for faces, as the size set + // in in reinit() is for all points if (data.update_each & update_covariant_transformation) data.covariant.resize(n_original_q_points); @@ -547,36 +572,6 @@ MappingFEField::compute_data( if (data.update_each & update_volume_elements) data.volume_elements.resize(n_original_q_points); - if (data.update_each & - (update_jacobian_grads | update_jacobian_pushed_forward_grads)) - data.shape_second_derivatives.resize(data.n_shape_functions * n_q_points); - - if (data.update_each & (update_jacobian_2nd_derivatives | - update_jacobian_pushed_forward_2nd_derivatives)) - data.shape_third_derivatives.resize(data.n_shape_functions * n_q_points); - - if (data.update_each & (update_jacobian_3rd_derivatives | - update_jacobian_pushed_forward_3rd_derivatives)) - data.shape_fourth_derivatives.resize(data.n_shape_functions * n_q_points); - - compute_shapes_virtual(q.get_points(), data); - - // This (for face values and simplices) can be different for different calls, - // so always copy - data.quadrature_weights = q.get_weights(); -} - - -template -void -MappingFEField::compute_face_data( - const UpdateFlags update_flags, - const Quadrature &q, - const unsigned int n_original_q_points, - InternalData &data) const -{ - compute_data(update_flags, q, n_original_q_points, data); - if (dim > 1) { if (data.update_each & update_boundary_forms) @@ -611,6 +606,7 @@ MappingFEField::compute_face_data( } + template typename std::unique_ptr::InternalDataBase> MappingFEField::get_data( @@ -619,8 +615,7 @@ MappingFEField::get_data( { std::unique_ptr::InternalDataBase> data_ptr = std::make_unique(euler_dof_handler->get_fe(), fe_mask); - auto &data = dynamic_cast(*data_ptr); - this->compute_data(update_flags, quadrature, quadrature.size(), data); + data_ptr->reinit(requires_update_flags(update_flags), quadrature); return data_ptr; } @@ -637,10 +632,12 @@ MappingFEField::get_face_data( std::unique_ptr::InternalDataBase> data_ptr = std::make_unique(euler_dof_handler->get_fe(), fe_mask); - auto &data = dynamic_cast(*data_ptr); + auto &data = dynamic_cast(*data_ptr); + const Quadrature q( QProjector::project_to_all_faces(reference_cell, quadrature[0])); - this->compute_face_data(update_flags, q, quadrature[0].size(), data); + data.reinit(requires_update_flags(update_flags), q); + this->compute_face_data(quadrature[0].size(), data); return data_ptr; } @@ -654,10 +651,12 @@ MappingFEField::get_subface_data( { std::unique_ptr::InternalDataBase> data_ptr = std::make_unique(euler_dof_handler->get_fe(), fe_mask); - auto &data = dynamic_cast(*data_ptr); + auto &data = dynamic_cast(*data_ptr); + const Quadrature q( QProjector::project_to_all_subfaces(reference_cell, quadrature)); - this->compute_face_data(update_flags, q, quadrature.size(), data); + data.reinit(requires_update_flags(update_flags), q); + this->compute_face_data(quadrature.size(), data); return data_ptr; } @@ -2198,9 +2197,9 @@ MappingFEField::transform_unit_to_real_cell( Assert(dynamic_cast(mdata.get()) != nullptr, ExcInternalError()); - update_internal_dofs(cell, dynamic_cast(*mdata)); + update_internal_dofs(cell, static_cast(*mdata)); - return do_transform_unit_to_real_cell(dynamic_cast(*mdata)); + return do_transform_unit_to_real_cell(static_cast(*mdata)); } @@ -2250,7 +2249,7 @@ MappingFEField::transform_real_to_unit_cell( initial_p_unit = cell->reference_cell().closest_point(initial_p_unit); - const Quadrature point_quadrature(initial_p_unit); + Quadrature point_quadrature(initial_p_unit); UpdateFlags update_flags = update_quadrature_points | update_jacobians; if (spacedim > dim) @@ -2260,12 +2259,12 @@ MappingFEField::transform_real_to_unit_cell( Assert(dynamic_cast(mdata.get()) != nullptr, ExcInternalError()); - update_internal_dofs(cell, dynamic_cast(*mdata)); + update_internal_dofs(cell, static_cast(*mdata)); return do_transform_real_to_unit_cell(cell, p, - initial_p_unit, - dynamic_cast(*mdata)); + point_quadrature, + static_cast(*mdata)); } @@ -2274,7 +2273,7 @@ Point MappingFEField::do_transform_real_to_unit_cell( const typename Triangulation::cell_iterator &cell, const Point &p, - const Point &initial_p_unit, + Quadrature &point_quadrature, InternalData &mdata) const { const unsigned int n_shapes = mdata.shape_values.size(); @@ -2291,10 +2290,12 @@ MappingFEField::do_transform_real_to_unit_cell( // The shape values and derivatives // of the mapping at this point are // previously computed. - // f(x) - Point p_unit = initial_p_unit; + + AssertDimension(point_quadrature.size(), 1); + Point p_unit = point_quadrature.point(0); Point f; - compute_shapes_virtual(std::vector>(1, p_unit), mdata); + mdata.reinit(mdata.update_each, point_quadrature); + Point p_real(do_transform_unit_to_real_cell(mdata)); Tensor<1, spacedim> p_minus_F = p - p_real; const double eps = 1.e-12 * cell->diameter(); @@ -2342,8 +2343,9 @@ MappingFEField::do_transform_real_to_unit_cell( p_unit_trial[i] -= step_length * delta[i]; // shape values and derivatives // at new p_unit point - compute_shapes_virtual(std::vector>(1, p_unit_trial), - mdata); + point_quadrature.initialize( + ArrayView>(&p_unit_trial, 1)); + mdata.reinit(mdata.update_each, point_quadrature); // f(x) Point p_real_trial = do_transform_unit_to_real_cell(mdata); const Tensor<1, spacedim> f_trial = p - p_real_trial; diff --git a/source/fe/mapping_manifold.cc b/source/fe/mapping_manifold.cc index 5ae5ffb3e6..551b83b87d 100644 --- a/source/fe/mapping_manifold.cc +++ b/source/fe/mapping_manifold.cc @@ -63,19 +63,21 @@ MappingManifold::InternalData::memory_consumption() const } + template void -MappingManifold::InternalData::initialize( +MappingManifold::InternalData::reinit( const UpdateFlags update_flags, - const Quadrature &q, - const unsigned int n_original_q_points) + const Quadrature &q) { // store the flags in the internal data object so we can access them // in fill_fe_*_values() this->update_each = update_flags; + const unsigned int n_q_points = q.size(); + // Store the quadrature - this->quad = q; + this->quad.initialize(q.get_points(), q.get_weights()); // Resize the weights this->vertex_weights.resize(GeometryInfo::vertices_per_cell); @@ -87,16 +89,17 @@ MappingManifold::InternalData::initialize( compute_manifold_quadrature_weights(q); if (this->update_each & update_covariant_transformation) - covariant.resize(n_original_q_points); + covariant.resize(n_q_points); if (this->update_each & update_contravariant_transformation) - contravariant.resize(n_original_q_points); + contravariant.resize(n_q_points); if (this->update_each & update_volume_elements) - volume_elements.resize(n_original_q_points); + volume_elements.resize(n_q_points); } + template void MappingManifold::InternalData::initialize_face( @@ -104,7 +107,18 @@ MappingManifold::InternalData::initialize_face( const Quadrature &q, const unsigned int n_original_q_points) { - initialize(update_flags, q, n_original_q_points); + reinit(update_flags, q); + + // Set to the size of a single quadrature object for faces, as the size set + // in in reinit() is for all points + if (this->update_each & update_covariant_transformation) + covariant.resize(n_original_q_points); + + if (this->update_each & update_contravariant_transformation) + contravariant.resize(n_original_q_points); + + if (this->update_each & update_volume_elements) + volume_elements.resize(n_original_q_points); if (dim > 1) { @@ -138,6 +152,38 @@ MappingManifold::InternalData::initialize_face( +template +void +MappingManifold::InternalData::store_vertices( + const typename Triangulation::cell_iterator &cell) const +{ + vertices.resize(GeometryInfo::vertices_per_cell); + for (const unsigned int i : GeometryInfo::vertex_indices()) + vertices[i] = cell->vertex(i); + this->cell = cell; +} + + + +template +void +MappingManifold::InternalData:: + compute_manifold_quadrature_weights(const Quadrature &quad) +{ + cell_manifold_quadrature_weights.resize( + quad.size(), std::vector(GeometryInfo::vertices_per_cell)); + for (unsigned int q = 0; q < quad.size(); ++q) + { + for (const unsigned int i : GeometryInfo::vertex_indices()) + { + cell_manifold_quadrature_weights[q][i] = + GeometryInfo::d_linear_shape_function(quad.point(q), i); + } + } +} + + + template MappingManifold::MappingManifold( const MappingManifold &) @@ -269,8 +315,7 @@ MappingManifold::get_data(const UpdateFlags update_flags, { std::unique_ptr::InternalDataBase> data_ptr = std::make_unique(); - auto &data = dynamic_cast(*data_ptr); - data.initialize(this->requires_update_flags(update_flags), q, q.size()); + data_ptr->reinit(this->requires_update_flags(update_flags), q); return data_ptr; } diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc index 5b77e81603..c51f40bf7e 100644 --- a/source/fe/mapping_q.cc +++ b/source/fe/mapping_q.cc @@ -78,10 +78,8 @@ MappingQ::InternalData::memory_consumption() const template void -MappingQ::InternalData::initialize( - const UpdateFlags update_flags, - const Quadrature &quadrature, - const unsigned int n_original_q_points) +MappingQ::InternalData::reinit(const UpdateFlags update_flags, + const Quadrature &quadrature) { // store the flags in the internal data object so we can access them // in fill_fe_*_values() @@ -90,7 +88,7 @@ MappingQ::InternalData::initialize( const unsigned int n_q_points = quadrature.size(); if (this->update_each & update_volume_elements) - volume_elements.resize(n_original_q_points); + volume_elements.resize(n_q_points); tensor_product_quadrature = quadrature.is_tensor_product(); @@ -163,7 +161,7 @@ MappingQ::InternalData::initialize_face( const Quadrature &quadrature, const unsigned int n_original_q_points) { - initialize(update_flags, quadrature, n_original_q_points); + reinit(update_flags, quadrature); quadrature_points = quadrature.get_points(); @@ -785,9 +783,7 @@ MappingQ::get_data(const UpdateFlags update_flags, { std::unique_ptr::InternalDataBase> data_ptr = std::make_unique(polynomial_degree); - auto &data = dynamic_cast(*data_ptr); - data.initialize(this->requires_update_flags(update_flags), q, q.size()); - + data_ptr->reinit(update_flags, q); return data_ptr; } -- 2.39.5