From b1cb46bff51a9cf670b236045f8e98719841d6fd Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 28 Oct 2015 16:33:36 -0500 Subject: [PATCH] Break dependency of MappingQ from MappingQGeneric. As discussed in #1732, MappingQ doesn't have the is-a property with regard to MappingQGeneric. Nor does it have it with regard to its previous base class, MappingQ1, and this was already fixed previously (#1429). Rather, what it should be is a has-a relationship with regard to both the MappingQGeneric (used on cells on the boundary) and MappingQ1 (used for interior cells). This patch implements this. The part of the patch for MappingQ itself is relatively straightforward: have it store pointers to both Q1 and Qp mappings, and where it currently dispatches to the base class MappingQGeneric, just dispatch to the Qp object instead. We can then break the inheritance from MappingQ to MappingQGeneric, at the cost of reimplementing a couple of trivial functions of the Mapping interface that were previously provided by the MappingQGeneric class. The bigger problems appear in the MappingQEulerian and MappingC1 classes which previously overloaded functions declared in MappingQGeneric and that are now no longer available to the MappingQGeneric object we now keep a pointer to, rather than have as base class. This is worked around by having both MappingQEulerian and MappingC1 declare their own internal classes derived from MappingQGeneric that provide the now missing functions. --- include/deal.II/fe/mapping_c1.h | 87 ++++++----- include/deal.II/fe/mapping_q.h | 58 +++++++- include/deal.II/fe/mapping_q_eulerian.h | 45 +++++- source/fe/mapping_c1.cc | 43 ++++-- source/fe/mapping_q.cc | 185 +++++++++++++++--------- source/fe/mapping_q_eulerian.cc | 51 +++++-- source/fe/mapping_q_generic.cc | 5 +- 7 files changed, 328 insertions(+), 146 deletions(-) diff --git a/include/deal.II/fe/mapping_c1.h b/include/deal.II/fe/mapping_c1.h index 56f2fd984b..b2ac7e7682 100644 --- a/include/deal.II/fe/mapping_c1.h +++ b/include/deal.II/fe/mapping_c1.h @@ -33,13 +33,10 @@ DEAL_II_NAMESPACE_OPEN * class chooses them such that the discretized boundary is globally * continuously differentiable. * - * To use this class, make sure that the Boundary::@p get_normals_at_vertices + * To use this class, make sure that the Boundary::get_normals_at_vertices() * function is implemented for the user's boundary object. * - * For more information about the spacedim template parameter check - * the documentation of FiniteElement or the one of Triangulation. - * - * @author Wolfgang Bangerth, 2001 + * @author Wolfgang Bangerth, 2001, 2015 */ template class MappingC1 : public MappingQ @@ -59,37 +56,53 @@ public: Mapping *clone () const; protected: - /** - * 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. - * - * Needed by the compute_support_points_simple(laplace) functions. - * For dim=1 this function is empty. - * - * 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; /** - * 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. - * - * Needed by the @p compute_support_points_laplace function. For - * dim=1 and 2 this function is empty. - * - * 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. + * 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. */ - virtual void - add_quad_support_points(const typename Triangulation::cell_iterator &cell, - std::vector > &a) const; + 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. + * + * Needed by the compute_support_points_simple(laplace) functions. + * For dim=1 this function is empty. + * + * 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; + + /** + * 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. + * + * Needed by the @p compute_support_points_laplace function. For + * dim=1 and 2 this function is empty. + * + * 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; + }; }; /*@}*/ @@ -98,17 +111,17 @@ protected: #ifndef DOXYGEN -template <> void MappingC1<1>::add_line_support_points ( +template <> void MappingC1<1>::MappingC1Generic::add_line_support_points ( const Triangulation<1>::cell_iterator &, std::vector > &) const; -template <> void MappingC1<2>::add_line_support_points ( +template <> void MappingC1<2>::MappingC1Generic::add_line_support_points ( const Triangulation<2>::cell_iterator &cell, std::vector > &a) const; -template <> void MappingC1<1>::add_quad_support_points ( +template <> void MappingC1<1>::MappingC1Generic::add_quad_support_points ( const Triangulation<1>::cell_iterator &, std::vector > &) const; -template <> void MappingC1<2>::add_quad_support_points ( +template <> void MappingC1<2>::MappingC1Generic::add_quad_support_points ( const Triangulation<2>::cell_iterator &, std::vector > &) const; diff --git a/include/deal.II/fe/mapping_q.h b/include/deal.II/fe/mapping_q.h index 86447fc6fe..545d6fc37a 100644 --- a/include/deal.II/fe/mapping_q.h +++ b/include/deal.II/fe/mapping_q.h @@ -76,7 +76,7 @@ template class TensorProductPolynomials; * @author Ralf Hartmann, 2000, 2001, 2005; Guido Kanschat 2000, 2001, Wolfgang Bangerth, 2015 */ template -class MappingQ : public MappingQGeneric +class MappingQ : public Mapping { public: /** @@ -104,6 +104,19 @@ public: */ 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; + /** * Transforms 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. @@ -215,13 +228,13 @@ protected: * also store a pointer to an InternalData object that pertains * to a $Q_1$ mapping. */ - class InternalData : public MappingQGeneric::InternalData + class InternalData : public Mapping::InternalDataBase { public: /** * Constructor. */ - InternalData (const unsigned int polynomial_degree); + InternalData (); /** @@ -242,10 +255,21 @@ protected: * $Q_1$ mapping that is, by default, used on all interior cells. */ std_cxx11::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_cxx11::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; + // documentation can be found in Mapping::get_data() virtual InternalData * @@ -296,6 +320,12 @@ protected: 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. @@ -316,13 +346,29 @@ protected: * 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_cxx11::unique_ptr > q1_mapping; + std_cxx11::shared_ptr > q1_mapping; /** - * Declare other MappingQ classes friends. + * 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. */ - template friend class MappingQ; + std_cxx11::shared_ptr > qp_mapping; }; /*@}*/ diff --git a/include/deal.II/fe/mapping_q_eulerian.h b/include/deal.II/fe/mapping_q_eulerian.h index 7be9b15f48..9964d19c07 100644 --- a/include/deal.II/fe/mapping_q_eulerian.h +++ b/include/deal.II/fe/mapping_q_eulerian.h @@ -211,13 +211,46 @@ private: mutable Threads::Mutex fe_values_mutex; /** - * Compute the positions of the support points in the current configuration. - * See the documentation of MappingQGeneric::compute_mapping_support_points() - * for more information. + * 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. */ - virtual - std::vector > - compute_mapping_support_points(const typename Triangulation::cell_iterator &cell) const; + class MappingQEulerianGeneric : public MappingQGeneric + { + 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. + */ + virtual + std_cxx11::array, GeometryInfo::vertices_per_cell> + get_vertices (const typename Triangulation::cell_iterator &cell) const; + + /** + * 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; + + private: + /** + * Reference to the surrounding object off of which we live. + */ + const MappingQEulerian &mapping_q_eulerian; + }; + }; /*@}*/ diff --git a/source/fe/mapping_c1.cc b/source/fe/mapping_c1.cc index 9e69706fdd..d16141702d 100644 --- a/source/fe/mapping_c1.cc +++ b/source/fe/mapping_c1.cc @@ -23,20 +23,35 @@ DEAL_II_NAMESPACE_OPEN +template +MappingC1::MappingC1Generic::MappingC1Generic () + : + MappingQGeneric (3) +{} + + + template MappingC1::MappingC1 () : MappingQ (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.reset (new MappingC1::MappingC1Generic()); } template <> void -MappingC1<1>::add_line_support_points (const Triangulation<1>::cell_iterator &, - std::vector > &) const +MappingC1<1>::MappingC1Generic::add_line_support_points (const Triangulation<1>::cell_iterator &, + std::vector > &) const { const unsigned int dim = 1; (void)dim; @@ -47,8 +62,8 @@ MappingC1<1>::add_line_support_points (const Triangulation<1>::cell_iterator &, template <> void -MappingC1<2>::add_line_support_points (const Triangulation<2>::cell_iterator &cell, - std::vector > &a) const +MappingC1<2>::MappingC1Generic::add_line_support_points (const Triangulation<2>::cell_iterator &cell, + std::vector > &a) const { const unsigned int dim = 2; std::vector > line_points (2); @@ -139,16 +154,16 @@ MappingC1<2>::add_line_support_points (const Triangulation<2>::cell_iterator &ce static const StraightBoundary straight_boundary; straight_boundary.get_intermediate_points_on_line (line, line_points); a.insert (a.end(), line_points.begin(), line_points.end()); - }; - }; + } + } } template void -MappingC1::add_line_support_points (const typename Triangulation::cell_iterator &, - std::vector > &) const +MappingC1::MappingC1Generic::add_line_support_points (const typename Triangulation::cell_iterator &, + std::vector > &) const { Assert (false, ExcNotImplemented()); } @@ -157,8 +172,8 @@ MappingC1::add_line_support_points (const typename Triangulation void -MappingC1<1>::add_quad_support_points (const Triangulation<1>::cell_iterator &, - std::vector > &) const +MappingC1<1>::MappingC1Generic::add_quad_support_points (const Triangulation<1>::cell_iterator &, + std::vector > &) const { const unsigned int dim = 1; (void)dim; @@ -169,8 +184,8 @@ MappingC1<1>::add_quad_support_points (const Triangulation<1>::cell_iterator &, template <> void -MappingC1<2>::add_quad_support_points (const Triangulation<2>::cell_iterator &, - std::vector > &) const +MappingC1<2>::MappingC1Generic::add_quad_support_points (const Triangulation<2>::cell_iterator &, + std::vector > &) const { const unsigned int dim = 2; (void)dim; @@ -181,8 +196,8 @@ MappingC1<2>::add_quad_support_points (const Triangulation<2>::cell_iterator &, template void -MappingC1::add_quad_support_points (const typename Triangulation::cell_iterator &, - std::vector > &) const +MappingC1::MappingC1Generic::add_quad_support_points (const typename Triangulation::cell_iterator &, + std::vector > &) const { Assert (false, ExcNotImplemented()); } diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc index 7236647f74..bd50a012d2 100644 --- a/source/fe/mapping_q.cc +++ b/source/fe/mapping_q.cc @@ -33,9 +33,8 @@ DEAL_II_NAMESPACE_OPEN template -MappingQ::InternalData::InternalData (const unsigned int polynomial_degree) +MappingQ::InternalData::InternalData () : - MappingQGeneric::InternalData(polynomial_degree), use_mapping_q1_on_current_cell(false) {} @@ -45,9 +44,10 @@ template std::size_t MappingQ::InternalData::memory_consumption () const { - return (MappingQGeneric::InternalData::memory_consumption () + + return (Mapping::InternalDataBase::memory_consumption () + MemoryConsumption::memory_consumption (use_mapping_q1_on_current_cell) + - MemoryConsumption::memory_consumption (mapping_q1_data)); + MemoryConsumption::memory_consumption (mapping_q1_data) + + MemoryConsumption::memory_consumption (mapping_qp_data)); } @@ -56,7 +56,7 @@ template MappingQ::MappingQ (const unsigned int degree, const bool use_mapping_q_on_all_cells) : - MappingQGeneric(degree), + 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 @@ -70,7 +70,15 @@ MappingQ::MappingQ (const unsigned int degree, (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 (new MappingQGeneric(1)) + q1_mapping (new MappingQGeneric(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_cxx11::shared_ptr >(new MappingQGeneric(degree)) + : + q1_mapping) {} @@ -78,34 +86,71 @@ MappingQ::MappingQ (const unsigned int degree, template MappingQ::MappingQ (const MappingQ &mapping) : - MappingQGeneric(mapping.get_degree()), + polynomial_degree (mapping.polynomial_degree), use_mapping_q_on_all_cells (mapping.use_mapping_q_on_all_cells), // clone the Q1 mapping for use on interior cells (if necessary) // or to create a good initial guess in transform_real_to_unit_cell() - q1_mapping (dynamic_cast*>(mapping.q1_mapping->clone())) + q1_mapping (dynamic_cast*>(mapping.q1_mapping->clone())), + // 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_cxx11::shared_ptr >(dynamic_cast*>(mapping.qp_mapping->clone())) + : + q1_mapping) {} +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 typename MappingQ::InternalData * MappingQ::get_data (const UpdateFlags update_flags, const Quadrature &quadrature) const { - InternalData *data = new InternalData(this->polynomial_degree); - - // fill the data of both the Q_p and the Q_1 objects in parallel - Threads::Task<> - task = Threads::new_task (std_cxx11::function - (std_cxx11::bind(&MappingQGeneric::InternalData::initialize, - data, - update_flags, - std_cxx11::cref(quadrature), - quadrature.size()))); + InternalData *data = new InternalData; + + // build the Q1 and Qp internal data objects in parallel + Threads::Task::InternalData *> + 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.reset (q1_mapping->get_data (update_flags, quadrature)); - task.join (); + // wait for the task above to finish and use returned value + data->mapping_qp_data.reset (do_get_data.return_value()); return data; } @@ -116,21 +161,20 @@ typename MappingQ::InternalData * MappingQ::get_face_data (const UpdateFlags update_flags, const Quadrature& quadrature) const { - InternalData *data = new InternalData(this->polynomial_degree); - const Quadrature q (QProjector::project_to_all_faces(quadrature)); - - // fill the data of both the Q_p and the Q_1 objects in parallel - Threads::Task<> - task = Threads::new_task (std_cxx11::function - (std_cxx11::bind(&MappingQGeneric::InternalData::initialize_face, - data, - update_flags, - std_cxx11::cref(q), - quadrature.size()))); + InternalData *data = new InternalData; + + // build the Q1 and Qp internal data objects in parallel + Threads::Task::InternalData *> + do_get_data = Threads::new_task (&MappingQGeneric::get_face_data, + *qp_mapping, + update_flags, + quadrature); + if (!use_mapping_q_on_all_cells) data->mapping_q1_data.reset (q1_mapping->get_face_data (update_flags, quadrature)); - task.join (); + // wait for the task above to finish and use returned value + data->mapping_qp_data.reset (do_get_data.return_value()); return data; } @@ -141,21 +185,20 @@ typename MappingQ::InternalData * MappingQ::get_subface_data (const UpdateFlags update_flags, const Quadrature& quadrature) const { - InternalData *data = new InternalData(this->polynomial_degree); - const Quadrature q (QProjector::project_to_all_subfaces(quadrature)); - - // fill the data of both the Q_p and the Q_1 objects in parallel - Threads::Task<> - task = Threads::new_task (std_cxx11::function - (std_cxx11::bind(&MappingQGeneric::InternalData::initialize_face, - data, - update_flags, - std_cxx11::cref(q), - quadrature.size()))); + InternalData *data = new InternalData; + + // build the Q1 and Qp internal data objects in parallel + Threads::Task::InternalData *> + 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.reset (q1_mapping->get_subface_data (update_flags, quadrature)); - task.join (); + // wait for the task above to finish and use returned value + data->mapping_qp_data.reset (do_get_data.return_value()); return data; } @@ -206,11 +249,11 @@ fill_fe_values (const typename Triangulation::cell_iterator &cell, *data.mapping_q1_data, output_data); else - MappingQGeneric::fill_fe_values(cell, - updated_cell_similarity, - quadrature, - data, - output_data); + qp_mapping->fill_fe_values(cell, + updated_cell_similarity, + quadrature, + *data.mapping_qp_data, + output_data); return updated_cell_similarity; } @@ -250,11 +293,11 @@ fill_fe_face_values (const typename Triangulation::cell_iterator & *data.mapping_q1_data, output_data); else - MappingQGeneric::fill_fe_face_values(cell, - face_no, - quadrature, - data, - output_data); + qp_mapping->fill_fe_face_values(cell, + face_no, + quadrature, + *data.mapping_qp_data, + output_data); } @@ -293,12 +336,12 @@ fill_fe_subface_values (const typename Triangulation::cell_iterato *data.mapping_q1_data, output_data); else - MappingQGeneric::fill_fe_subface_values(cell, - face_no, - subface_no, - quadrature, - data, - output_data); + qp_mapping->fill_fe_subface_values(cell, + face_no, + subface_no, + quadrature, + *data.mapping_qp_data, + output_data); } @@ -318,10 +361,10 @@ transform (const VectorSlice > > input, // check whether we should in fact work on the Q1 portion of it if (data->use_mapping_q1_on_current_cell) - return q1_mapping->transform (input, mapping_type, *data->mapping_q1_data, output); + q1_mapping->transform (input, mapping_type, *data->mapping_q1_data, output); else // otherwise use the full mapping - MappingQGeneric::transform(input, mapping_type, mapping_data, output); + qp_mapping->transform(input, mapping_type, *data->mapping_qp_data, output); } @@ -342,10 +385,10 @@ transform (const VectorSlice // check whether we should in fact work on the Q1 portion of it if (data->use_mapping_q1_on_current_cell) - return q1_mapping->transform (input, mapping_type, *data->mapping_q1_data, output); + q1_mapping->transform (input, mapping_type, *data->mapping_q1_data, output); else // otherwise use the full mapping - MappingQGeneric::transform(input, mapping_type, mapping_data, output); + qp_mapping->transform(input, mapping_type, *data->mapping_qp_data, output); } @@ -364,10 +407,10 @@ transform (const VectorSlice > > input, // check whether we should in fact work on the Q1 portion of it if (data->use_mapping_q1_on_current_cell) - return q1_mapping->transform (input, mapping_type, *data->mapping_q1_data, output); + q1_mapping->transform (input, mapping_type, *data->mapping_q1_data, output); else // otherwise use the full mapping - MappingQGeneric::transform(input, mapping_type, mapping_data, output); + qp_mapping->transform(input, mapping_type, *data->mapping_qp_data, output); } @@ -388,10 +431,10 @@ transform (const VectorSlice // check whether we should in fact work on the Q1 portion of it if (data->use_mapping_q1_on_current_cell) - return q1_mapping->transform (input, mapping_type, *data->mapping_q1_data, output); + q1_mapping->transform (input, mapping_type, *data->mapping_q1_data, output); else // otherwise use the full mapping - MappingQGeneric::transform(input, mapping_type, mapping_data, output); + qp_mapping->transform(input, mapping_type, *data->mapping_qp_data, output); } @@ -410,10 +453,10 @@ transform (const VectorSlice > > input, // check whether we should in fact work on the Q1 portion of it if (data->use_mapping_q1_on_current_cell) - return q1_mapping->transform (input, mapping_type, *data->mapping_q1_data, output); + q1_mapping->transform (input, mapping_type, *data->mapping_q1_data, output); else // otherwise use the full mapping - MappingQGeneric::transform(input, mapping_type, mapping_data, output); + qp_mapping->transform(input, mapping_type, *data->mapping_qp_data, output); } @@ -427,7 +470,7 @@ transform_unit_to_real_cell (const typename Triangulation::cell_it // 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 this->MappingQGeneric::transform_unit_to_real_cell (cell, p); + return qp_mapping->transform_unit_to_real_cell (cell, p); else return q1_mapping->transform_unit_to_real_cell (cell, p); } @@ -445,7 +488,7 @@ transform_real_to_unit_cell (const typename Triangulation::cell_it use_mapping_q_on_all_cells || (dim!=spacedim) ) - return MappingQGeneric::transform_real_to_unit_cell(cell, p); + return qp_mapping->transform_real_to_unit_cell(cell, p); else return q1_mapping->transform_real_to_unit_cell(cell, p); } diff --git a/source/fe/mapping_q_eulerian.cc b/source/fe/mapping_q_eulerian.cc index 56d4cbad0d..60213ecb6b 100644 --- a/source/fe/mapping_q_eulerian.cc +++ b/source/fe/mapping_q_eulerian.cc @@ -34,6 +34,15 @@ 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) +{} + template MappingQEulerian:: MappingQEulerian (const unsigned int degree, @@ -53,6 +62,9 @@ MappingQEulerian (const unsigned int degree, // current vector this->q1_mapping.reset (new MappingQ1Eulerian(euler_vector, euler_dof_handler)); + + // also reset the qp mapping pointer with our own class + this->qp_mapping.reset (new MappingQEulerianGeneric(degree,*this)); } @@ -76,6 +88,9 @@ MappingQEulerian (const unsigned int degree, // current vector this->q1_mapping.reset (new MappingQ1Eulerian(euler_vector, euler_dof_handler)); + + // also reset the qp mapping pointer with our own class + this->qp_mapping.reset (new MappingQEulerianGeneric(degree,*this)); } @@ -133,7 +148,8 @@ get_vertices (const typename Triangulation::cell_iterator &cell) const { // get the vertices as the first 2^dim mapping support points - const std::vector > a = compute_mapping_support_points(cell); + const std::vector > a + = dynamic_cast(*this->qp_mapping).compute_mapping_support_points(cell); std_cxx11::array, GeometryInfo::vertices_per_cell> vertex_locations; std::copy (a.begin(), @@ -145,15 +161,27 @@ get_vertices +template +std_cxx11::array, GeometryInfo::vertices_per_cell> +MappingQEulerian::MappingQEulerianGeneric:: +get_vertices +(const typename Triangulation::cell_iterator &cell) const +{ + return mapping_q_eulerian.get_vertices (cell); +} + + + + template std::vector > -MappingQEulerian:: +MappingQEulerian::MappingQEulerianGeneric:: compute_mapping_support_points (const typename Triangulation::cell_iterator &cell) const { // first, basic assertion with respect to vector size, - const types::global_dof_index n_dofs = euler_dof_handler->n_dofs(); - const types::global_dof_index vector_size = euler_vector->size(); + const types::global_dof_index n_dofs = mapping_q_eulerian.euler_dof_handler->n_dofs(); + const types::global_dof_index vector_size = mapping_q_eulerian.euler_vector->size(); (void)n_dofs; (void)vector_size; @@ -161,7 +189,8 @@ compute_mapping_support_points (const typename Triangulation::cell // 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, euler_dof_handler); + typename DoFHandler::cell_iterator dof_cell(*cell, + mapping_q_eulerian.euler_dof_handler); Assert (dof_cell->active() == true, ExcInactiveCell()); @@ -177,8 +206,8 @@ compute_mapping_support_points (const typename Triangulation::cell // 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 = euler_dof_handler->get_fe().n_components(); + const unsigned int n_support_pts = mapping_q_eulerian.support_quadrature.size(); + const unsigned int n_components = mapping_q_eulerian.euler_dof_handler->get_fe().n_components(); Assert (n_components >= spacedim, ExcDimensionMismatch(n_components, spacedim) ); @@ -189,16 +218,16 @@ compute_mapping_support_points (const typename Triangulation::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 - Threads::Mutex::ScopedLock lock(fe_values_mutex); - fe_values.reinit(dof_cell); - fe_values.get_function_values(*euler_vector, shift_vector); + Threads::Mutex::ScopedLock lock(mapping_q_eulerian.fe_values_mutex); + mapping_q_eulerian.fe_values.reinit(dof_cell); + mapping_q_eulerian.fe_values.get_function_values(*mapping_q_eulerian.euler_vector, shift_vector); // and finally compute the positions of the support points in the deformed // configuration. std::vector > a(n_support_pts); for (unsigned int q=0; q::MappingQGeneric (const unsigned int p) fe_q(dim == 3 ? new FE_Q(this->polynomial_degree) : 0), support_point_weights_on_quad (compute_support_point_weights_on_quad(this->polynomial_degree)), support_point_weights_on_hex (compute_support_point_weights_on_hex(this->polynomial_degree)) -{} +{ + Assert (p >= 1, ExcMessage ("It only makes sense to create polynomial mappings " + "with a polynomial degree greater or equal to one.")); +} -- 2.39.5