From: Martin Kronbichler Date: Wed, 9 Jun 2021 18:08:07 +0000 (+0200) Subject: Make the class MappingQ equivalent to MappingQGeneric X-Git-Tag: v9.4.0-rc1~1155^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e73d5268f18f49c7c87be0ca00268f82a2962b06;p=dealii.git Make the class MappingQ equivalent to MappingQGeneric --- diff --git a/include/deal.II/fe/mapping_c1.h b/include/deal.II/fe/mapping_c1.h index b5f98850dd..621d4547d1 100644 --- a/include/deal.II/fe/mapping_c1.h +++ b/include/deal.II/fe/mapping_c1.h @@ -28,14 +28,14 @@ DEAL_II_NAMESPACE_OPEN /** * Mapping class that uses C1 (continuously differentiable) cubic mappings of - * the boundary. This class is built atop of MappingQ by simply determining - * the interpolation points for a cubic mapping of the boundary differently: - * MappingQ chooses them such that they interpolate the boundary, while this - * class chooses them such that the discretized boundary is globally - * continuously differentiable. + * the boundary. This class is built atop of MappingQGeneric by simply + * determining the interpolation points for a cubic mapping of the boundary + * differently: MappingQ chooses them such that they interpolate the boundary, + * while this class chooses them such that the discretized boundary is + * globally continuously differentiable. */ template -class MappingC1 : public MappingQ +class MappingC1 : public MappingQGeneric { public: /** @@ -51,48 +51,33 @@ public: virtual std::unique_ptr> clone() const override; -protected: /** - * A class derived from MappingQGeneric that provides the generic mapping - * with support points on boundary objects so that the corresponding Q3 - * mapping ends up being C1. + * For dim=2,3. Append the support points of all shape functions + * located on bounding lines to the vector @p a. Points located on the + * line but on vertices are not included. + * + * This function chooses the respective points not such that they are + * interpolating the boundary (as does the base class), but rather such + * that the resulting cubic mapping is a continuous one. */ - class MappingC1Generic : public MappingQGeneric - { - public: - /** - * Constructor. - */ - MappingC1Generic(); - - /** - * For dim=2,3. Append the support points of all shape functions - * located on bounding lines to the vector @p a. Points located on the - * line but on vertices are not included. - * - * This function chooses the respective points not such that they are - * interpolating the boundary (as does the base class), but rather such - * that the resulting cubic mapping is a continuous one. - */ - virtual void - add_line_support_points( - const typename Triangulation::cell_iterator &cell, - std::vector> & a) const override; - - /** - * For dim=3. Append the support points of all shape functions - * located on bounding faces (quads in 3d) to the vector @p a. Points - * located on the line but on vertices are not included. - * - * This function chooses the respective points not such that they are - * interpolating the boundary (as does the base class), but rather such - * that the resulting cubic mapping is a continuous one. - */ - virtual void - add_quad_support_points( - const typename Triangulation::cell_iterator &cell, - std::vector> & a) const override; - }; + virtual void + add_line_support_points( + const typename Triangulation::cell_iterator &cell, + std::vector> & a) const override; + + /** + * For dim=3. Append the support points of all shape functions + * located on bounding faces (quads in 3d) to the vector @p a. Points + * located on the line but on vertices are not included. + * + * This function chooses the respective points not such that they are + * interpolating the boundary (as does the base class), but rather such + * that the resulting cubic mapping is a continuous one. + */ + virtual void + add_quad_support_points( + const typename Triangulation::cell_iterator &cell, + std::vector> & a) const override; }; /*@}*/ @@ -103,25 +88,22 @@ protected: template <> void -MappingC1<1>::MappingC1Generic::add_line_support_points( - const Triangulation<1>::cell_iterator &, - std::vector> &) const; +MappingC1<1>::add_line_support_points(const Triangulation<1>::cell_iterator &, + std::vector> &) const; template <> void -MappingC1<2>::MappingC1Generic::add_line_support_points( +MappingC1<2>::add_line_support_points( const Triangulation<2>::cell_iterator &cell, std::vector> & a) const; template <> void -MappingC1<1>::MappingC1Generic::add_quad_support_points( - const Triangulation<1>::cell_iterator &, - std::vector> &) const; +MappingC1<1>::add_quad_support_points(const Triangulation<1>::cell_iterator &, + std::vector> &) const; template <> void -MappingC1<2>::MappingC1Generic::add_quad_support_points( - const Triangulation<2>::cell_iterator &, - std::vector> &) const; +MappingC1<2>::add_quad_support_points(const Triangulation<2>::cell_iterator &, + std::vector> &) const; #endif // DOXYGEN diff --git a/include/deal.II/fe/mapping_q.h b/include/deal.II/fe/mapping_q.h index abe2f7b43b..5c8262aaa4 100644 --- a/include/deal.II/fe/mapping_q.h +++ b/include/deal.II/fe/mapping_q.h @@ -27,354 +27,32 @@ DEAL_II_NAMESPACE_OPEN /*@{*/ /** - * A class that implements a polynomial mapping $Q_p$ of degree $p$ on cells - * at the boundary of the domain (or, if requested in the constructor, for all - * cells) and linear mappings for interior cells. - * - * The class is in fact poorly named since (unless explicitly specified during - * the construction of the object, see below), it does not actually use - * mappings of degree $p$ everywhere, but only on cells at the - * boundary. This is in contrast to the MappingQGeneric class which indeed - * does use a polynomial mapping $Q_p$ of degree $p$ everywhere. The point of - * the current class is that in many situations, curved domains are only - * provided with information about how exactly edges at the boundary are - * shaped, but we do not know anything about internal edges. Thus, in the - * absence of other information, we can only assume that internal edges are - * straight lines, and in that case internal cells may as well be treated is - * bilinear quadrilaterals or trilinear hexahedra. (An example of how such - * meshes look is shown in step-1 already, but it is also discussed in the - * "Results" section of step-6.) Because bi-/trilinear mappings are - * significantly cheaper to compute than higher order mappings, it is - * advantageous in such situations to use the higher order mapping only on - * cells at the boundary of the domain. This class implements exactly this - * behavior. - * - * There are a number of special cases worth considering: - * - If you want to use a higher order mapping for all cells, you can - * achieve this by setting the second argument to the constructor to true. - * This only makes sense if you can actually provide information about how - * interior edges and faces of the mesh should be curved. This is typically - * done by associating a Manifold with interior cells and edges. A simple - * example of this is discussed in the "Results" section of step-6; a full - * discussion of manifolds is provided in step-53. - * - If you pass true as the second argument to this class, then it - * is in fact completely equivalent to generating a MappingQGeneric object - * right away. - * - This class is also entirely equivalent to MappingQGeneric if the - * polynomial degree provided is one. This is because in that case, no - * distinction between the mapping used on cells in the interior and on the - * boundary of the domain can be made. - * - If you are working on meshes embedded in higher space dimensions, - * i.e., if dim!=spacedim, then every cell is considered to be at the boundary - * of the domain and consequently a higher order mapping is used for all - * cells; again this class is then equivalent to using MappingQGeneric right - * away. - * - *

Behavior along curved boundaries and with different manifolds

- * - * For the behavior of the mapping and convergence rates in case of mixing - * different manifolds, please consult the respective section of - * MappingQGeneric. + * A class that implements a polynomial mapping $Q_p$ of degree $p$ on all + * cells. This class is completely equivalent to the MappingQGeneric class. */ template -class MappingQ : public Mapping +class MappingQ : public MappingQGeneric { public: /** * Constructor. @p polynomial_degree denotes the polynomial degree of the * polynomials that are used to map cells boundary. - * - * The second argument determines whether the higher order mapping should - * also be used on interior cells. If its value is false (the - * default), then a lower order mapping is used in the interior. This is - * sufficient for most cases where higher order mappings are only used to - * better approximate the boundary. In that case, cells bounded by straight - * lines are acceptable in the interior. However, there are cases where one - * would also like to use a higher order mapping in the interior. The - * MappingQEulerian class is one such case. - * - * The value of @p use_mapping_q_on_all_cells is ignored if @p dim is not - * equal to @p spacedim, i.e., if we are considering meshes on surfaces - * embedded into higher dimensional spaces. - */ - MappingQ(const unsigned int polynomial_degree, - const bool use_mapping_q_on_all_cells = false); - - /** - * Copy constructor. - */ - MappingQ(const MappingQ &mapping); - - /** - * Return the degree of the mapping, i.e. the value which was passed to the - * constructor. - */ - unsigned int - get_degree() const; - - /** - * Always returns @p true because the default implementation of functions in - * this class preserves vertex locations. - */ - virtual bool - preserves_vertex_locations() const override; - - // for documentation, see the Mapping base class - virtual BoundingBox - get_bounding_box(const typename Triangulation::cell_iterator - &cell) const override; - - virtual bool - is_compatible_with(const ReferenceCell &reference_cell) const override; - - /** - * Transform the point @p p on the unit cell to the point @p p_real on the - * real cell @p cell and returns @p p_real. - */ - virtual Point - transform_unit_to_real_cell( - const typename Triangulation::cell_iterator &cell, - const Point &p) const override; - - /** - * Transform the point @p p on the real cell to the point @p p_unit on the - * unit cell @p cell and returns @p p_unit. - * - * Uses Newton iteration and the @p transform_unit_to_real_cell function. - * - * In the codimension one case, this function returns the normal projection - * of the real point @p p on the curve or surface identified by the @p cell. - * - * @note Polynomial mappings from the reference (unit) cell coordinates to - * the coordinate system of a real cell are not always invertible if the - * point for which the inverse mapping is to be computed lies outside the - * cell's boundaries. In such cases, the current function may fail to - * compute a point on the reference cell whose image under the mapping - * equals the given point @p p. If this is the case then this function - * throws an exception of type Mapping::ExcTransformationFailed . Whether - * the given point @p p lies outside the cell can therefore be determined by - * checking whether the return reference coordinates lie inside of outside - * the reference cell (e.g., using GeometryInfo::is_inside_unit_cell) or - * whether the exception mentioned above has been thrown. - */ - virtual Point - transform_real_to_unit_cell( - const typename Triangulation::cell_iterator &cell, - const Point &p) const override; - - // for documentation, see the Mapping base class - virtual void - transform(const ArrayView> & input, - const MappingKind kind, - const typename Mapping::InternalDataBase &internal, - const ArrayView> &output) const override; - - // for documentation, see the Mapping base class - virtual void - transform(const ArrayView> &input, - const MappingKind kind, - const typename Mapping::InternalDataBase &internal, - const ArrayView> &output) const override; - - // for documentation, see the Mapping base class - virtual void - transform(const ArrayView> & input, - const MappingKind kind, - const typename Mapping::InternalDataBase &internal, - const ArrayView> &output) const override; - - // for documentation, see the Mapping base class - virtual void - transform(const ArrayView> &input, - const MappingKind kind, - const typename Mapping::InternalDataBase &internal, - const ArrayView> &output) const override; - - // for documentation, see the Mapping base class - virtual void - transform(const ArrayView> & input, - const MappingKind kind, - const typename Mapping::InternalDataBase &internal, - const ArrayView> &output) const override; - - /** - * Return a pointer to a copy of the present object. The caller of this copy - * then assumes ownership of it. - */ - - virtual std::unique_ptr> - clone() const override; - - - /** - * @name Interface with FEValues - * @{ */ + MappingQ(const unsigned int polynomial_degree); -protected: /** - * Storage for internal data of this mapping. See Mapping::InternalDataBase - * for an extensive description. - * - * This includes data that is computed once when the object is created (in - * get_data()) as well as data the class wants to store from between the - * call to fill_fe_values(), fill_fe_face_values(), or - * fill_fe_subface_values() until possible later calls from the finite - * element to functions such as transform(). The latter class of member - * variables are marked as 'mutable'. - * - * The current class uses essentially the same fields for storage as the - * MappingQGeneric class. Consequently, it inherits from - * MappingQGeneric::InternalData, rather than from - * Mapping::InternalDataBase. The principal difference to - * MappingQGeneric::InternalData is that MappingQ switches between $Q_1$ and - * $Q_p$ mappings depending on the cell we are on, so the internal data - * object needs to also store a pointer to an InternalData object that - * pertains to a $Q_1$ mapping. + * The second argument is here for backward compatibility with previous + * versions of deal.II, but it does not have any effect on the workings of + * this class. */ - class InternalData : public Mapping::InternalDataBase - { - public: - /** - * Constructor. - */ - InternalData(); - - - /** - * Return an estimate (in bytes) for the memory consumption of this object. - */ - virtual std::size_t - memory_consumption() const override; - - /** - * Flag that is set by the fill_fe_[[sub]face]_values function. - * - * If this flag is @p true we are on an interior cell and the @p - * mapping_q1_data is used. - */ - mutable bool use_mapping_q1_on_current_cell; - - /** - * A pointer to a structure to store the information for the pure $Q_1$ - * mapping that is, by default, used on all interior cells. - */ - std::unique_ptr::InternalData> - mapping_q1_data; - - /** - * A pointer to a structure to store the information for the full $Q_p$ - * mapping that is, by default, used on all boundary cells. - */ - std::unique_ptr::InternalData> - mapping_qp_data; - }; - -protected: - // documentation can be found in Mapping::requires_update_flags() - virtual UpdateFlags - requires_update_flags(const UpdateFlags update_flags) const override; - - // documentation can be found in Mapping::get_data() - virtual std::unique_ptr::InternalDataBase> - get_data(const UpdateFlags, const Quadrature &quadrature) const override; - - using Mapping::get_face_data; - - // documentation can be found in Mapping::get_face_data() - virtual std::unique_ptr::InternalDataBase> - get_face_data(const UpdateFlags flags, - const hp::QCollection &quadrature) const override; - - // documentation can be found in Mapping::get_subface_data() - virtual std::unique_ptr::InternalDataBase> - get_subface_data(const UpdateFlags flags, - const Quadrature &quadrature) const override; - - // documentation can be found in Mapping::fill_fe_values() - virtual CellSimilarity::Similarity - fill_fe_values( - const typename Triangulation::cell_iterator &cell, - const CellSimilarity::Similarity cell_similarity, - const Quadrature & quadrature, - const typename Mapping::InternalDataBase & internal_data, - internal::FEValuesImplementation::MappingRelatedData - &output_data) const override; - - using Mapping::fill_fe_face_values; - - // documentation can be found in Mapping::fill_fe_face_values() - virtual void - fill_fe_face_values( - const typename Triangulation::cell_iterator &cell, - const unsigned int face_no, - const hp::QCollection & quadrature, - const typename Mapping::InternalDataBase & internal_data, - internal::FEValuesImplementation::MappingRelatedData - &output_data) const override; - - // documentation can be found in Mapping::fill_fe_subface_values() - virtual void - fill_fe_subface_values( - const typename Triangulation::cell_iterator &cell, - const unsigned int face_no, - const unsigned int subface_no, - const Quadrature & quadrature, - const typename Mapping::InternalDataBase & internal_data, - internal::FEValuesImplementation::MappingRelatedData - &output_data) const override; - - /** - * @} - */ - -protected: - /** - * The polynomial degree of the cells to be used on all cells at the - * boundary of the domain, or everywhere if so specified. - */ - const unsigned int polynomial_degree; - - /** - * If this flag is set @p true then @p MappingQ is used on all cells, not - * only on boundary cells. - */ - const bool use_mapping_q_on_all_cells; - - /** - * Pointer to a Q1 mapping. This mapping is used on interior cells unless - * use_mapping_q_on_all_cells was set in the call to the constructor. The - * mapping is also used on any cell in the transform_real_to_unit_cell() to - * compute a cheap initial guess for the position of the point before we - * employ the more expensive Newton iteration using the full mapping. - * - * @note MappingQEulerian resets this pointer to an object of type - * MappingQ1Eulerian to ensure that the Q1 mapping also knows about the - * proper shifts and transformations of the Eulerian displacements. This - * also means that we really need to store our own Q1 mapping here, rather - * than simply resorting to StaticMappingQ1::mapping. - * - * @note If the polynomial degree used for the current object is one, then - * the qp_mapping and q1_mapping variables point to the same underlying - * object. - */ - std::shared_ptr> q1_mapping; + DEAL_II_DEPRECATED_EARLY + MappingQ(const unsigned int polynomial_degree, + const bool use_mapping_q_on_all_cells); /** - * Pointer to a Q_p mapping. This mapping is used on boundary cells unless - * use_mapping_q_on_all_cells was set in the call to the constructor (in - * which case it is used for all cells). - * - * @note MappingQEulerian and MappingC1 reset this pointer to an object of - * their own implementation to ensure that the Q_p mapping also knows about - * the proper shifts and transformations of the Eulerian displacements - * (Eulerian case) and proper choice of support points (C1 case). - * - * @note If the polynomial degree used for the current object is one, then - * the qp_mapping and q1_mapping variables point to the same underlying - * object. + * Copy constructor. */ - std::shared_ptr> qp_mapping; + MappingQ(const MappingQ &mapping); }; /*@}*/ diff --git a/include/deal.II/fe/mapping_q_eulerian.h b/include/deal.II/fe/mapping_q_eulerian.h index 1e19e037d2..bd1a4686e6 100644 --- a/include/deal.II/fe/mapping_q_eulerian.h +++ b/include/deal.II/fe/mapping_q_eulerian.h @@ -90,7 +90,7 @@ class Vector; * of the vector can be specified as template parameter VectorType. */ template , int spacedim = dim> -class MappingQEulerian : public MappingQ +class MappingQEulerian : public MappingQGeneric { public: /** @@ -141,9 +141,10 @@ public: preserves_vertex_locations() const override; // for documentation, see the Mapping base class - virtual BoundingBox - get_bounding_box(const typename Triangulation::cell_iterator - &cell) const override; + virtual std::vector> + compute_mapping_support_points( + const typename Triangulation::cell_iterator &cell) + const override; /** * Exception which is thrown when the mapping is being evaluated at @@ -182,7 +183,6 @@ protected: MappingQEulerian> euler_dof_handler; - private: /** * Multigrid level at which the mapping is to be used. @@ -190,89 +190,37 @@ private: const unsigned int level; /** - * A class derived from MappingQGeneric that provides the generic mapping - * with support points on boundary objects so that the corresponding Q3 - * mapping ends up being C1. + * Special quadrature rule used to define the support points in the + * reference configuration. */ - class MappingQEulerianGeneric : public MappingQGeneric + class SupportQuadrature : public Quadrature { public: /** - * Constructor. - */ - MappingQEulerianGeneric( - const unsigned int degree, - const MappingQEulerian &mapping_q_eulerian); - - /** - * Return the mapped vertices of the cell. For the current class, this - * function does not use the support points from the geometry of the - * current cell but instead evaluates an externally given displacement - * field in addition to the geometry of the cell. + * Constructor, with an argument defining the desired polynomial degree. */ - virtual boost::container::small_vector, - GeometryInfo::vertices_per_cell> - get_vertices(const typename Triangulation::cell_iterator - &cell) const override; - - /** - * Compute the positions of the support points in the current - * configuration. See the documentation of - * MappingQGeneric::compute_mapping_support_points() for more information. - */ - virtual std::vector> - compute_mapping_support_points( - const typename Triangulation::cell_iterator &cell) - const override; - - /** - * Always return @p false because MappingQEulerianGeneric does not in general - * preserve vertex locations (unless the translation vector happens to - * provide for zero displacements at vertex locations). - */ - virtual bool - preserves_vertex_locations() const override; - - private: - /** - * Reference to the surrounding object off of which we live. - */ - const MappingQEulerian &mapping_q_eulerian; - - - /** - * Special quadrature rule used to define the support points in the - * reference configuration. - */ - class SupportQuadrature : public Quadrature - { - public: - /** - * Constructor, with an argument defining the desired polynomial degree. - */ - SupportQuadrature(const unsigned int map_degree); - }; + SupportQuadrature(const unsigned int map_degree); + }; - /** - * A member variable holding the quadrature points in the right order. - */ - const SupportQuadrature support_quadrature; + /** + * A member variable holding the quadrature points in the right order. + */ + const SupportQuadrature support_quadrature; - /** - * FEValues object used to query the given finite element field at the - * support points in the reference configuration. - * - * The variable is marked as mutable since we have to call - * FEValues::reinit from compute_mapping_support_points, a function that - * is 'const'. - */ - mutable FEValues fe_values; + /** + * FEValues object used to query the given finite element field at the + * support points in the reference configuration. + * + * The variable is marked as mutable since we have to call + * FEValues::reinit from compute_mapping_support_points, a function that + * is 'const'. + */ + mutable FEValues fe_values; - /** - * A variable to guard access to the fe_values variable. - */ - mutable Threads::Mutex fe_values_mutex; - }; + /** + * A variable to guard access to the fe_values variable. + */ + mutable Threads::Mutex fe_values_mutex; }; /*@}*/ diff --git a/source/fe/mapping_c1.cc b/source/fe/mapping_c1.cc index de998fc4d7..d5165e3887 100644 --- a/source/fe/mapping_c1.cc +++ b/source/fe/mapping_c1.cc @@ -27,35 +27,19 @@ DEAL_II_NAMESPACE_OPEN -template -MappingC1::MappingC1Generic::MappingC1Generic() - : MappingQGeneric(3) -{} - - - template MappingC1::MappingC1() - : MappingQ(3) + : MappingQGeneric(3) { Assert(dim > 1, ExcImpossibleInDim(dim)); - - // replace the mapping_qp objects of the base class by something - // that knows about generating data points based on the geometry - // - // we only need to replace the Qp mapping because that's the one that's - // used on boundary cells where it matters - this->qp_mapping = - std::make_shared::MappingC1Generic>(); } template <> void -MappingC1<1>::MappingC1Generic::add_line_support_points( - const Triangulation<1>::cell_iterator &, - std::vector> &) const +MappingC1<1>::add_line_support_points(const Triangulation<1>::cell_iterator &, + std::vector> &) const { const unsigned int dim = 1; (void)dim; @@ -66,7 +50,7 @@ MappingC1<1>::MappingC1Generic::add_line_support_points( template <> void -MappingC1<2>::MappingC1Generic::add_line_support_points( +MappingC1<2>::add_line_support_points( const Triangulation<2>::cell_iterator &cell, std::vector> & a) const { @@ -172,7 +156,7 @@ MappingC1<2>::MappingC1Generic::add_line_support_points( template void -MappingC1::MappingC1Generic::add_line_support_points( +MappingC1::add_line_support_points( const typename Triangulation::cell_iterator &, std::vector> &) const { @@ -183,9 +167,8 @@ MappingC1::MappingC1Generic::add_line_support_points( template <> void -MappingC1<1>::MappingC1Generic::add_quad_support_points( - const Triangulation<1>::cell_iterator &, - std::vector> &) const +MappingC1<1>::add_quad_support_points(const Triangulation<1>::cell_iterator &, + std::vector> &) const { const unsigned int dim = 1; (void)dim; @@ -196,9 +179,8 @@ MappingC1<1>::MappingC1Generic::add_quad_support_points( template <> void -MappingC1<2>::MappingC1Generic::add_quad_support_points( - const Triangulation<2>::cell_iterator &, - std::vector> &) const +MappingC1<2>::add_quad_support_points(const Triangulation<2>::cell_iterator &, + std::vector> &) const { const unsigned int dim = 2; (void)dim; @@ -209,7 +191,7 @@ MappingC1<2>::MappingC1Generic::add_quad_support_points( template void -MappingC1::MappingC1Generic::add_quad_support_points( +MappingC1::add_quad_support_points( const typename Triangulation::cell_iterator &, std::vector> &) const { diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc index 06da719bb1..aaf181dd8f 100644 --- a/source/fe/mapping_q.cc +++ b/source/fe/mapping_q.cc @@ -38,542 +38,23 @@ DEAL_II_NAMESPACE_OPEN template -MappingQ::InternalData::InternalData() - : use_mapping_q1_on_current_cell(false) +MappingQ::MappingQ(const unsigned int degree) + : MappingQGeneric(degree) {} template -std::size_t -MappingQ::InternalData::memory_consumption() const -{ - return ( - Mapping::InternalDataBase::memory_consumption() + - MemoryConsumption::memory_consumption(use_mapping_q1_on_current_cell) + - MemoryConsumption::memory_consumption(mapping_q1_data) + - MemoryConsumption::memory_consumption(mapping_qp_data)); -} - - - -template -MappingQ::MappingQ(const unsigned int degree, - const bool use_mapping_q_on_all_cells) - : polynomial_degree(degree) - , - - // see whether we want to use *this* mapping objects on *all* cells, - // or defer to an explicit Q1 mapping on interior cells. if - // degree==1, then we are already that Q1 mapping, so we don't need - // it; if dim!=spacedim, there is also no need for anything because - // we're most likely on a curved manifold - use_mapping_q_on_all_cells(degree == 1 || use_mapping_q_on_all_cells || - (dim != spacedim)) - , - // create a Q1 mapping for use on interior cells (if necessary) - // or to create a good initial guess in transform_real_to_unit_cell() - q1_mapping(std::make_shared>(1)) - , - - // create a Q_p mapping; if p=1, simply share the Q_1 mapping already - // created via the shared_ptr objects - qp_mapping(this->polynomial_degree > 1 ? - std::make_shared>(degree) : - q1_mapping) +MappingQ::MappingQ(const unsigned int degree, const bool) + : MappingQGeneric(degree) {} template MappingQ::MappingQ(const MappingQ &mapping) - : polynomial_degree(mapping.polynomial_degree) - , use_mapping_q_on_all_cells(mapping.use_mapping_q_on_all_cells) -{ - // Note that we really do have to use clone() here, since mapping.q1_mapping - // may be MappingQ1Eulerian and mapping.qp_mapping may be MappingQEulerian. - std::shared_ptr> other_q1_map = - mapping.q1_mapping->clone(); - q1_mapping = std::dynamic_pointer_cast>( - other_q1_map); - Assert(q1_mapping != nullptr, ExcInternalError()); - Assert(q1_mapping->get_degree() == 1, ExcInternalError()); - - // Same as the other constructor: if possible reuse the Q1 mapping - if (this->polynomial_degree == 1) - { - qp_mapping = q1_mapping; - } - else - { - std::shared_ptr> other_qp_map = - mapping.qp_mapping->clone(); - qp_mapping = - std::dynamic_pointer_cast>( - other_qp_map); - Assert(qp_mapping != nullptr, ExcInternalError()); - } -} - - - -template -unsigned int -MappingQ::get_degree() const -{ - return polynomial_degree; -} - - - -template -inline bool -MappingQ::preserves_vertex_locations() const -{ - return true; -} - - - -template -UpdateFlags -MappingQ::requires_update_flags(const UpdateFlags in) const -{ - return (q1_mapping->requires_update_flags(in) | - qp_mapping->requires_update_flags(in)); -} - - - -template -std::unique_ptr::InternalDataBase> -MappingQ::get_data(const UpdateFlags update_flags, - const Quadrature &quadrature) const -{ - std::unique_ptr::InternalDataBase> data_ptr = - std::make_unique(); - auto &data = dynamic_cast(*data_ptr); - - // build the Q1 and Qp internal data objects in parallel - Threads::Task< - std::unique_ptr::InternalDataBase>> - do_get_data = Threads::new_task(&MappingQGeneric::get_data, - *qp_mapping, - update_flags, - quadrature); - - if (!use_mapping_q_on_all_cells) - data.mapping_q1_data = Utilities::dynamic_unique_cast< - typename MappingQGeneric::InternalData>( - std::move(q1_mapping->get_data(update_flags, quadrature))); - - // wait for the task above to finish and use returned value - data.mapping_qp_data = Utilities::dynamic_unique_cast< - typename MappingQGeneric::InternalData>( - std::move(do_get_data.return_value())); - return data_ptr; -} - - - -template -std::unique_ptr::InternalDataBase> -MappingQ::get_face_data( - const UpdateFlags update_flags, - const hp::QCollection &quadrature) const -{ - std::unique_ptr::InternalDataBase> data_ptr = - std::make_unique(); - auto &data = dynamic_cast(*data_ptr); - - std::unique_ptr::InternalDataBase> ( - MappingQGeneric::*mapping_get_face_data)( - const UpdateFlags, const hp::QCollection &) const = - &MappingQGeneric::get_face_data; - - // build the Q1 and Qp internal data objects in parallel - Threads::Task< - std::unique_ptr::InternalDataBase>> - do_get_data = Threads::new_task(mapping_get_face_data, - *qp_mapping, - update_flags, - quadrature); - - if (!use_mapping_q_on_all_cells) - data.mapping_q1_data = Utilities::dynamic_unique_cast< - typename MappingQGeneric::InternalData>( - std::move(q1_mapping->get_face_data(update_flags, quadrature))); - - // wait for the task above to finish and use returned value - data.mapping_qp_data = Utilities::dynamic_unique_cast< - typename MappingQGeneric::InternalData>( - std::move(do_get_data.return_value())); - return data_ptr; -} - - - -template -std::unique_ptr::InternalDataBase> -MappingQ::get_subface_data( - const UpdateFlags update_flags, - const Quadrature &quadrature) const -{ - std::unique_ptr::InternalDataBase> data_ptr = - std::make_unique(); - auto &data = dynamic_cast(*data_ptr); - - // build the Q1 and Qp internal data objects in parallel - Threads::Task< - std::unique_ptr::InternalDataBase>> - do_get_data = - Threads::new_task(&MappingQGeneric::get_subface_data, - *qp_mapping, - update_flags, - quadrature); - - if (!use_mapping_q_on_all_cells) - data.mapping_q1_data = Utilities::dynamic_unique_cast< - typename MappingQGeneric::InternalData>( - std::move(q1_mapping->get_subface_data(update_flags, quadrature))); - - // wait for the task above to finish and use returned value - data.mapping_qp_data = Utilities::dynamic_unique_cast< - typename MappingQGeneric::InternalData>( - std::move(do_get_data.return_value())); - return data_ptr; -} - - -// Note that the CellSimilarity flag is modifiable, since MappingQ can need to -// recalculate data even when cells are similar. -template -CellSimilarity::Similarity -MappingQ::fill_fe_values( - const typename Triangulation::cell_iterator &cell, - const CellSimilarity::Similarity cell_similarity, - const Quadrature & quadrature, - const typename Mapping::InternalDataBase & internal_data, - internal::FEValuesImplementation::MappingRelatedData - &output_data) const -{ - // 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); - - // check whether this cell needs the full mapping or can be treated by a - // reduced Q1 mapping, e.g. if the cell is in the interior of the domain - data.use_mapping_q1_on_current_cell = - !(use_mapping_q_on_all_cells || cell->has_boundary_lines()); - - - // call the base class. we need to ensure that the flag indicating whether - // we can use some similarity has to be modified - for a general MappingQ, - // the data needs to be recomputed anyway since then the mapping changes the - // data. this needs to be known also for later operations, so modify the - // variable here. this also affects the calculation of the next cell -- if - // we use Q1 data on the next cell, the data will still be invalid. - const CellSimilarity::Similarity updated_cell_similarity = - ((data.use_mapping_q1_on_current_cell == false) && - (this->polynomial_degree > 1) ? - CellSimilarity::invalid_next_cell : - cell_similarity); - - // depending on the results above, decide whether the Q1 mapping or - // the Qp mapping needs to handle this cell - if (data.use_mapping_q1_on_current_cell) - q1_mapping->fill_fe_values(cell, - updated_cell_similarity, - quadrature, - *data.mapping_q1_data, - output_data); - else - qp_mapping->fill_fe_values(cell, - updated_cell_similarity, - quadrature, - *data.mapping_qp_data, - output_data); - - return updated_cell_similarity; -} - - - -template -void -MappingQ::fill_fe_face_values( - const typename Triangulation::cell_iterator &cell, - const unsigned int face_no, - const hp::QCollection & quadrature, - const typename Mapping::InternalDataBase & internal_data, - internal::FEValuesImplementation::MappingRelatedData - &output_data) const -{ - // 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); - - // check whether this cell needs the full mapping or can be treated by a - // reduced Q1 mapping, e.g. if the cell is entirely in the interior of the - // domain. note that it is not sufficient to ask whether the present _face_ - // is in the interior, as the mapping on the face depends on the mapping of - // the cell, which in turn depends on the fact whether _any_ of the faces of - // this cell is at the boundary, not only the present face - data.use_mapping_q1_on_current_cell = - !(use_mapping_q_on_all_cells || cell->has_boundary_lines()); - - // depending on the results above, decide whether the Q1 mapping or - // the Qp mapping needs to handle this cell - if (data.use_mapping_q1_on_current_cell) - q1_mapping->fill_fe_face_values( - cell, face_no, quadrature, *data.mapping_q1_data, output_data); - else - qp_mapping->fill_fe_face_values( - cell, face_no, quadrature, *data.mapping_qp_data, output_data); -} - - -template -void -MappingQ::fill_fe_subface_values( - const typename Triangulation::cell_iterator &cell, - const unsigned int face_no, - const unsigned int subface_no, - const Quadrature & quadrature, - const typename Mapping::InternalDataBase & internal_data, - internal::FEValuesImplementation::MappingRelatedData - &output_data) const -{ - // 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); - - // check whether this cell needs the full mapping or can be treated by a - // reduced Q1 mapping, e.g. if the cell is entirely in the interior of the - // domain. note that it is not sufficient to ask whether the present _face_ - // is in the interior, as the mapping on the face depends on the mapping of - // the cell, which in turn depends on the fact whether _any_ of the faces of - // this cell is at the boundary, not only the present face - data.use_mapping_q1_on_current_cell = - !(use_mapping_q_on_all_cells || cell->has_boundary_lines()); - - // depending on the results above, decide whether the Q1 mapping or - // the Qp mapping needs to handle this cell - if (data.use_mapping_q1_on_current_cell) - q1_mapping->fill_fe_subface_values(cell, - face_no, - subface_no, - quadrature, - *data.mapping_q1_data, - output_data); - else - qp_mapping->fill_fe_subface_values(cell, - face_no, - subface_no, - quadrature, - *data.mapping_qp_data, - output_data); -} - - - -template -void -MappingQ::transform( - const ArrayView> & input, - const MappingKind mapping_kind, - const typename Mapping::InternalDataBase &mapping_data, - const ArrayView> & output) const -{ - AssertDimension(input.size(), output.size()); - - const InternalData *data = dynamic_cast(&mapping_data); - Assert(data != nullptr, ExcInternalError()); - - // check whether we should in fact work on the Q1 portion of it - if (data->use_mapping_q1_on_current_cell) - q1_mapping->transform(input, mapping_kind, *data->mapping_q1_data, output); - else - // otherwise use the full mapping - qp_mapping->transform(input, mapping_kind, *data->mapping_qp_data, output); -} - - - -template -void -MappingQ::transform( - const ArrayView> &input, - const MappingKind mapping_kind, - const typename Mapping::InternalDataBase &mapping_data, - const ArrayView> & output) const -{ - AssertDimension(input.size(), output.size()); - Assert((dynamic_cast::InternalData *>( - &mapping_data) != nullptr), - ExcInternalError()); - const InternalData *data = dynamic_cast(&mapping_data); - - // check whether we should in fact work on the Q1 portion of it - if (data->use_mapping_q1_on_current_cell) - q1_mapping->transform(input, mapping_kind, *data->mapping_q1_data, output); - else - // otherwise use the full mapping - qp_mapping->transform(input, mapping_kind, *data->mapping_qp_data, output); -} - - - -template -void -MappingQ::transform( - const ArrayView> & input, - const MappingKind mapping_kind, - const typename Mapping::InternalDataBase &mapping_data, - const ArrayView> & output) const -{ - AssertDimension(input.size(), output.size()); - Assert((dynamic_cast::InternalData *>( - &mapping_data) != nullptr), - ExcInternalError()); - const InternalData *data = dynamic_cast(&mapping_data); - - // check whether we should in fact work on the Q1 portion of it - if (data->use_mapping_q1_on_current_cell) - q1_mapping->transform(input, mapping_kind, *data->mapping_q1_data, output); - else - // otherwise use the full mapping - qp_mapping->transform(input, mapping_kind, *data->mapping_qp_data, output); -} - - - -template -void -MappingQ::transform( - const ArrayView> &input, - const MappingKind mapping_kind, - const typename Mapping::InternalDataBase &mapping_data, - const ArrayView> & output) const -{ - AssertDimension(input.size(), output.size()); - Assert((dynamic_cast::InternalData *>( - &mapping_data) != nullptr), - ExcInternalError()); - const InternalData *data = dynamic_cast(&mapping_data); - - // check whether we should in fact work on the Q1 portion of it - if (data->use_mapping_q1_on_current_cell) - q1_mapping->transform(input, mapping_kind, *data->mapping_q1_data, output); - else - // otherwise use the full mapping - qp_mapping->transform(input, mapping_kind, *data->mapping_qp_data, output); -} - - - -template -void -MappingQ::transform( - const ArrayView> & input, - const MappingKind mapping_kind, - const typename Mapping::InternalDataBase &mapping_data, - const ArrayView> & output) const -{ - AssertDimension(input.size(), output.size()); - Assert((dynamic_cast::InternalData *>( - &mapping_data) != nullptr), - ExcInternalError()); - const InternalData *data = dynamic_cast(&mapping_data); - - // check whether we should in fact work on the Q1 portion of it - if (data->use_mapping_q1_on_current_cell) - q1_mapping->transform(input, mapping_kind, *data->mapping_q1_data, output); - else - // otherwise use the full mapping - qp_mapping->transform(input, mapping_kind, *data->mapping_qp_data, output); -} - - - -template -Point -MappingQ::transform_unit_to_real_cell( - const typename Triangulation::cell_iterator &cell, - const Point & p) const -{ - // first see, whether we want to use a linear or a higher order - // mapping, then either use our own facilities or that of the Q1 - // mapping we store - if (use_mapping_q_on_all_cells || cell->has_boundary_lines()) - return qp_mapping->transform_unit_to_real_cell(cell, p); - else - return q1_mapping->transform_unit_to_real_cell(cell, p); -} - - - -template -Point -MappingQ::transform_real_to_unit_cell( - const typename Triangulation::cell_iterator &cell, - const Point & p) const -{ - if (cell->has_boundary_lines() || use_mapping_q_on_all_cells || - (dim != spacedim)) - return qp_mapping->transform_real_to_unit_cell(cell, p); - else - return q1_mapping->transform_real_to_unit_cell(cell, p); -} - - - -template -std::unique_ptr> -MappingQ::clone() const -{ - return std::make_unique>( - this->polynomial_degree, this->use_mapping_q_on_all_cells); -} - - - -template -BoundingBox -MappingQ::get_bounding_box( - const typename Triangulation::cell_iterator &cell) const -{ - if (cell->has_boundary_lines() || use_mapping_q_on_all_cells || - (dim != spacedim)) - return BoundingBox( - qp_mapping->compute_mapping_support_points(cell)); - else - return BoundingBox(q1_mapping->get_vertices(cell)); -} - - - -template -bool -MappingQ::is_compatible_with( - const ReferenceCell &reference_cell) const -{ - Assert(dim == reference_cell.get_dimension(), - ExcMessage("The dimension of your mapping (" + - Utilities::to_string(dim) + - ") and the reference cell cell_type (" + - Utilities::to_string(reference_cell.get_dimension()) + - " ) do not agree.")); - - return reference_cell.is_hyper_cube(); -} - + : MappingQGeneric(mapping) +{} // explicit instantiations diff --git a/source/fe/mapping_q_eulerian.cc b/source/fe/mapping_q_eulerian.cc index 98675809de..9102152678 100644 --- a/source/fe/mapping_q_eulerian.cc +++ b/source/fe/mapping_q_eulerian.cc @@ -44,20 +44,6 @@ DEAL_II_NAMESPACE_OPEN // .... MAPPING Q EULERIAN CONSTRUCTOR -template -MappingQEulerian::MappingQEulerianGeneric:: - MappingQEulerianGeneric( - const unsigned int degree, - const MappingQEulerian &mapping_q_eulerian) - : MappingQGeneric(degree) - , mapping_q_eulerian(mapping_q_eulerian) - , support_quadrature(degree) - , fe_values(mapping_q_eulerian.euler_dof_handler->get_fe(), - support_quadrature, - update_values | update_quadrature_points) -{} - - template MappingQEulerian::MappingQEulerian( @@ -65,21 +51,15 @@ MappingQEulerian::MappingQEulerian( const DoFHandler &euler_dof_handler, const VectorType & euler_vector, const unsigned int level) - : MappingQ(degree, true) + : MappingQGeneric(degree) , euler_vector(&euler_vector) , euler_dof_handler(&euler_dof_handler) , level(level) -{ - // reset the q1 mapping we use for interior cells (and previously - // set by the MappingQ constructor) to a MappingQ1Eulerian with the - // current vector - this->q1_mapping = - std::make_shared>( - euler_dof_handler, euler_vector); - - // also reset the qp mapping pointer with our own class - this->qp_mapping = std::make_shared(degree, *this); -} + , support_quadrature(degree) + , fe_values(euler_dof_handler.get_fe(), + support_quadrature, + update_values | update_quadrature_points) +{} @@ -96,8 +76,8 @@ MappingQEulerian::clone() const // .... SUPPORT QUADRATURE CONSTRUCTOR template -MappingQEulerian::MappingQEulerianGeneric:: - SupportQuadrature::SupportQuadrature(const unsigned int map_degree) +MappingQEulerian::SupportQuadrature:: + SupportQuadrature(const unsigned int map_degree) : Quadrature(Utilities::fixed_power(map_degree + 1)) { // first we determine the support points on the unit cell in lexicographic @@ -127,9 +107,7 @@ MappingQEulerian::get_vertices( const typename Triangulation::cell_iterator &cell) const { // get the vertices as the first 2^dim mapping support points - const std::vector> a = - dynamic_cast(*this->qp_mapping) - .compute_mapping_support_points(cell); + const std::vector> a = compute_mapping_support_points(cell); boost::container::small_vector, GeometryInfo::vertices_per_cell> @@ -141,43 +119,16 @@ MappingQEulerian::get_vertices( -template -boost::container::small_vector, - GeometryInfo::vertices_per_cell> -MappingQEulerian::MappingQEulerianGeneric:: - get_vertices( - const typename Triangulation::cell_iterator &cell) const -{ - return mapping_q_eulerian.get_vertices(cell); -} - - - -template -bool -MappingQEulerian::MappingQEulerianGeneric:: - preserves_vertex_locations() const -{ - return false; -} - - - template std::vector> -MappingQEulerian::MappingQEulerianGeneric:: - compute_mapping_support_points( - const typename Triangulation::cell_iterator &cell) const +MappingQEulerian::compute_mapping_support_points( + const typename Triangulation::cell_iterator &cell) const { - const bool mg_vector = - mapping_q_eulerian.level != numbers::invalid_unsigned_int; + const bool mg_vector = level != numbers::invalid_unsigned_int; const types::global_dof_index n_dofs = - mg_vector ? - mapping_q_eulerian.euler_dof_handler->n_dofs(mapping_q_eulerian.level) : - mapping_q_eulerian.euler_dof_handler->n_dofs(); - const types::global_dof_index vector_size = - mapping_q_eulerian.euler_vector->size(); + mg_vector ? euler_dof_handler->n_dofs(level) : euler_dof_handler->n_dofs(); + const types::global_dof_index vector_size = euler_vector->size(); (void)n_dofs; (void)vector_size; @@ -186,8 +137,8 @@ MappingQEulerian::MappingQEulerianGeneric:: // we then transform our tria iterator into a dof iterator so we can access // data not associated with triangulations - typename DoFHandler::cell_iterator dof_cell( - *cell, mapping_q_eulerian.euler_dof_handler); + typename DoFHandler::cell_iterator dof_cell(*cell, + euler_dof_handler); Assert(mg_vector || dof_cell->is_active() == true, ExcInactiveCell()); @@ -203,8 +154,7 @@ MappingQEulerian::MappingQEulerianGeneric:: // that. This implies that the user should order components appropriately, // or create a separate dof handler for the displacements. const unsigned int n_support_pts = support_quadrature.size(); - const unsigned int n_components = - mapping_q_eulerian.euler_dof_handler->get_fe(0).n_components(); + const unsigned int n_components = euler_dof_handler->get_fe(0).n_components(); Assert(n_components >= spacedim, ExcDimensionMismatch(n_components, spacedim)); @@ -213,7 +163,7 @@ MappingQEulerian::MappingQEulerianGeneric:: n_support_pts, Vector(n_components)); std::vector dof_indices( - mapping_q_eulerian.euler_dof_handler->get_fe(0).n_dofs_per_cell()); + euler_dof_handler->get_fe(0).n_dofs_per_cell()); // fill shift vector for each support point using an fe_values object. make // sure that the fe_values variable isn't used simultaneously from different // threads @@ -222,13 +172,10 @@ MappingQEulerian::MappingQEulerianGeneric:: if (mg_vector) { dof_cell->get_mg_dof_indices(dof_indices); - fe_values.get_function_values(*mapping_q_eulerian.euler_vector, - dof_indices, - shift_vector); + fe_values.get_function_values(*euler_vector, dof_indices, shift_vector); } else - fe_values.get_function_values(*mapping_q_eulerian.euler_vector, - shift_vector); + fe_values.get_function_values(*euler_vector, shift_vector); // and finally compute the positions of the support points in the deformed // configuration. @@ -258,11 +205,12 @@ MappingQEulerian::fill_fe_values( // call the function of the base class, but ignoring // any potentially detected cell similarity between // the current and the previous cell - MappingQ::fill_fe_values(cell, - CellSimilarity::invalid_next_cell, - quadrature, - internal_data, - output_data); + MappingQGeneric::fill_fe_values( + cell, + CellSimilarity::invalid_next_cell, + quadrature, + internal_data, + output_data); // also return the updated flag since any detected // similarity wasn't based on the mapped field, but // the original vertices which are meaningless @@ -270,19 +218,6 @@ MappingQEulerian::fill_fe_values( } - -template -BoundingBox -MappingQEulerian::get_bounding_box( - const typename Triangulation::cell_iterator &cell) const -{ - return BoundingBox( - dynamic_cast(*this->qp_mapping) - .compute_mapping_support_points(cell)); -} - - - // explicit instantiations #include "mapping_q_eulerian.inst"