From 306a739e62608bdf3d0aaad0b6e2a4ecd936905e Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Sun, 6 Sep 2015 16:47:01 -0500 Subject: [PATCH] Let MappingQ have a MappingQ1 member. This mostly breaks the derivation of MappingQ from MappingQ1. There are a few hold-outs left, though. --- include/deal.II/fe/mapping_q.h | 21 ++- include/deal.II/fe/mapping_q_generic.h | 15 +- source/fe/mapping_c1.cc | 2 +- source/fe/mapping_q.cc | 216 +++++++++++++------------ source/fe/mapping_q_generic.cc | 4 +- 5 files changed, 148 insertions(+), 110 deletions(-) diff --git a/include/deal.II/fe/mapping_q.h b/include/deal.II/fe/mapping_q.h index b954602e3d..b2c49db3e4 100644 --- a/include/deal.II/fe/mapping_q.h +++ b/include/deal.II/fe/mapping_q.h @@ -62,16 +62,25 @@ public: * * The second argument determines whether the higher order mapping should * also be used on interior cells. If its value is false (the - * default), the a lower-order mapping is used in the interior. This is + * 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 ignore 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); + /** * 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. @@ -277,13 +286,13 @@ protected: // documentation can be found in Mapping::get_face_data() virtual - typename Mapping::InternalDataBase * + InternalData * get_face_data (const UpdateFlags flags, const Quadrature& quadrature) const; // documentation can be found in Mapping::get_subface_data() virtual - typename Mapping::InternalDataBase * + InternalData * get_subface_data (const UpdateFlags flags, const Quadrature& quadrature) const; @@ -440,6 +449,12 @@ protected: */ const FE_Q feq; + /** + * Pointer to a Q1 mapping to be used on interior cells unless + * use_mapping_q_on_all_cells was set in the call to the + * constructor. + */ + std_cxx11::unique_ptr > q1_mapping; /* * The default line support points. These are used when computing diff --git a/include/deal.II/fe/mapping_q_generic.h b/include/deal.II/fe/mapping_q_generic.h index 3791c7da33..cb5df7fd60 100644 --- a/include/deal.II/fe/mapping_q_generic.h +++ b/include/deal.II/fe/mapping_q_generic.h @@ -29,6 +29,9 @@ DEAL_II_NAMESPACE_OPEN +template class MappingQ; + + /*!@addtogroup mapping */ /*@{*/ @@ -374,7 +377,6 @@ public: mutable std::vector volume_elements; }; -protected: // documentation can be found in Mapping::requires_update_flags() virtual @@ -389,13 +391,13 @@ protected: // documentation can be found in Mapping::get_face_data() virtual - typename Mapping::InternalDataBase * + InternalData * get_face_data (const UpdateFlags flags, const Quadrature& quadrature) const; // documentation can be found in Mapping::get_subface_data() virtual - typename Mapping::InternalDataBase * + InternalData * get_subface_data (const UpdateFlags flags, const Quadrature& quadrature) const; @@ -429,6 +431,7 @@ protected: * @} */ +protected: /** * The degree of the polynomials used as shape functions for the mapping @@ -450,6 +453,12 @@ protected: void compute_mapping_support_points (const typename Triangulation::cell_iterator &cell, std::vector > &a) const = 0; + + /** + * Make MappingQ a friend since it needs to call the + * fill_fe_values() functions on its MappingQ1 sub-object. + */ + template friend class MappingQ; }; diff --git a/source/fe/mapping_c1.cc b/source/fe/mapping_c1.cc index 52e53a9b33..9e69706fdd 100644 --- a/source/fe/mapping_c1.cc +++ b/source/fe/mapping_c1.cc @@ -193,7 +193,7 @@ template Mapping * MappingC1::clone () const { - return new MappingC1(*this); + return new MappingC1(); } diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc index 5eb71043b0..ef90792511 100644 --- a/source/fe/mapping_q.cc +++ b/source/fe/mapping_q.cc @@ -36,7 +36,7 @@ DEAL_II_NAMESPACE_OPEN template MappingQ::InternalData::InternalData (const unsigned int polynomial_degree) : - MappingQ1::InternalData(polynomial_degree), + MappingQGeneric::InternalData(polynomial_degree), use_mapping_q1_on_current_cell(false) {} @@ -66,6 +66,11 @@ MappingQ::MappingQ (const unsigned int degree, use_mapping_q_on_all_cells (use_mapping_q_on_all_cells || (dim != spacedim)), feq(degree), + q1_mapping (this->use_mapping_q_on_all_cells + ? + 0 + : + new MappingQ1()), line_support_points(degree+1) { Assert(n_inner+n_outer==Utilities::fixed_power(degree+1), @@ -83,6 +88,36 @@ MappingQ::MappingQ (const unsigned int degree, +template +MappingQ::MappingQ (const MappingQ &mapping) + : + MappingQ1(mapping.get_degree()), + n_inner(mapping.n_inner), + n_outer(mapping.n_outer), + use_mapping_q_on_all_cells (mapping.use_mapping_q_on_all_cells), + feq(mapping.get_degree()), + q1_mapping (mapping.q1_mapping != 0 + ? + dynamic_cast*>(mapping.q1_mapping->clone()) + : + 0), + line_support_points(mapping.line_support_points) +{ + Assert(n_inner+n_outer==Utilities::fixed_power(this->polynomial_degree+1), + ExcInternalError()); + + // build laplace_on_quad_vector + if (this->polynomial_degree>1) + { + if (dim >= 2) + set_laplace_on_quad_vector(laplace_on_quad_vector); + if (dim >= 3) + set_laplace_on_hex_vector(laplace_on_hex_vector); + } +} + + + template typename MappingQ::InternalData * MappingQ::get_data (const UpdateFlags update_flags, @@ -91,32 +126,24 @@ MappingQ::get_data (const UpdateFlags update_flags, InternalData *data = new InternalData(this->polynomial_degree); // fill the data of both the Q_p and the Q_1 objects in parallel - Threads::TaskGroup<> tasks; - tasks += Threads::new_task (std_cxx11::function - (std_cxx11::bind(&MappingQ1::InternalData::initialize, - data, - update_flags, - std_cxx11::cref(quadrature), - quadrature.size()))); + Threads::Task<> + task = Threads::new_task (std_cxx11::function + (std_cxx11::bind(&MappingQ1::InternalData::initialize, + data, + update_flags, + std_cxx11::cref(quadrature), + quadrature.size()))); if (!use_mapping_q_on_all_cells) - { - data->mapping_q1_data.reset (new typename MappingQ1::InternalData(1)); - tasks += Threads::new_task (std_cxx11::function - (std_cxx11::bind(&MappingQ1::InternalData::initialize, - data->mapping_q1_data.get(), - update_flags, - std_cxx11::cref(quadrature), - quadrature.size()))); - } + data->mapping_q1_data.reset (q1_mapping->get_data (update_flags, quadrature)); - tasks.join_all (); + task.join (); return data; } template -typename Mapping::InternalDataBase * +typename MappingQ::InternalData * MappingQ::get_face_data (const UpdateFlags update_flags, const Quadrature& quadrature) const { @@ -124,32 +151,24 @@ MappingQ::get_face_data (const UpdateFlags update_flags, 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::TaskGroup<> tasks; - tasks += Threads::new_task (std_cxx11::function - (std_cxx11::bind(&MappingQ1::InternalData::initialize_face, - data, - update_flags, - std_cxx11::cref(q), - quadrature.size()))); + Threads::Task<> + task = Threads::new_task (std_cxx11::function + (std_cxx11::bind(&MappingQ1::InternalData::initialize_face, + data, + update_flags, + std_cxx11::cref(q), + quadrature.size()))); if (!use_mapping_q_on_all_cells) - { - data->mapping_q1_data.reset (new typename MappingQ1::InternalData(1)); - tasks += Threads::new_task (std_cxx11::function - (std_cxx11::bind(&MappingQ1::InternalData::initialize_face, - data->mapping_q1_data.get(), - update_flags, - std_cxx11::cref(q), - quadrature.size()))); - } + data->mapping_q1_data.reset (q1_mapping->get_face_data (update_flags, quadrature)); - tasks.join_all (); + task.join (); return data; } template -typename Mapping::InternalDataBase * +typename MappingQ::InternalData * MappingQ::get_subface_data (const UpdateFlags update_flags, const Quadrature& quadrature) const { @@ -157,25 +176,17 @@ MappingQ::get_subface_data (const UpdateFlags update_flags, 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::TaskGroup<> tasks; - tasks += Threads::new_task (std_cxx11::function - (std_cxx11::bind(&MappingQ1::InternalData::initialize_face, - data, - update_flags, - std_cxx11::cref(q), - quadrature.size()))); + Threads::Task<> + task = Threads::new_task (std_cxx11::function + (std_cxx11::bind(&MappingQ1::InternalData::initialize_face, + data, + update_flags, + std_cxx11::cref(q), + quadrature.size()))); if (!use_mapping_q_on_all_cells) - { - data->mapping_q1_data.reset (new typename MappingQ1::InternalData(1)); - tasks += Threads::new_task (std_cxx11::function - (std_cxx11::bind(&MappingQ1::InternalData::initialize_face, - data->mapping_q1_data.get(), - update_flags, - std_cxx11::cref(q), - quadrature.size()))); - } + data->mapping_q1_data.reset (q1_mapping->get_subface_data (update_flags, quadrature)); - tasks.join_all (); + task.join (); return data; } @@ -201,14 +212,6 @@ fill_fe_values (const typename Triangulation::cell_iterator &cell, data.use_mapping_q1_on_current_cell = !(use_mapping_q_on_all_cells || cell->has_boundary_lines()); - // depending on this result, use this or the other data object for the - // mapping. - const typename MappingQ1::InternalData * - p_data = (data.use_mapping_q1_on_current_cell - ? - data.mapping_q1_data.get() - : - &data); // 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, @@ -224,11 +227,21 @@ fill_fe_values (const typename Triangulation::cell_iterator &cell, CellSimilarity::invalid_next_cell : cell_similarity); - MappingQ1::fill_fe_values(cell, - updated_cell_similarity, - quadrature, - *p_data, - output_data); + + // 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 + MappingQ1::fill_fe_values(cell, + updated_cell_similarity, + quadrature, + data, + output_data); return updated_cell_similarity; } @@ -259,20 +272,20 @@ fill_fe_face_values (const typename Triangulation::cell_iterator & data.use_mapping_q1_on_current_cell = !(use_mapping_q_on_all_cells || cell->has_boundary_lines()); - // depending on this result, use this or the other data object for the - // mapping and call the function in the base class with this - // data object to do the work - const typename MappingQ1::InternalData *p_data - = (data.use_mapping_q1_on_current_cell - ? - data.mapping_q1_data.get() - : - &data); - MappingQ1::fill_fe_face_values(cell, - face_no, - quadrature, - *p_data, - output_data); + // 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 + MappingQ1::fill_fe_face_values(cell, + face_no, + quadrature, + data, + output_data); } @@ -301,21 +314,22 @@ fill_fe_subface_values (const typename Triangulation::cell_iterato data.use_mapping_q1_on_current_cell = !(use_mapping_q_on_all_cells || cell->has_boundary_lines()); - // depending on this result, use this or the other data object for the - // mapping and call the function in the base class with this - // data object to do the work - const typename MappingQ1::InternalData *p_data - = (data.use_mapping_q1_on_current_cell - ? - data.mapping_q1_data.get() - : - &data); - MappingQ1::fill_fe_subface_values(cell, - face_no, - subface_no, - quadrature, - *p_data, - output_data); + // 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 + MappingQ1::fill_fe_subface_values(cell, + face_no, + subface_no, + quadrature, + data, + output_data); } @@ -865,7 +879,7 @@ transform (const VectorSlice > > input, // whether we should in fact work on the Q1 portion of it if (const InternalData *data = dynamic_cast(&mapping_data)) if (data->use_mapping_q1_on_current_cell) - return MappingQ1::transform (input, mapping_type, *data->mapping_q1_data, output); + return q1_mapping->transform (input, mapping_type, *data->mapping_q1_data, output); // otherwise just stick with what we already had MappingQ1::transform(input, mapping_type, mapping_data, output); @@ -890,7 +904,7 @@ transform (const VectorSlice // whether we should in fact work on the Q1 portion of it if (const InternalData *data = dynamic_cast(&mapping_data)) if (data->use_mapping_q1_on_current_cell) - return MappingQ1::transform (input, mapping_type, *data->mapping_q1_data, output); + return q1_mapping->transform (input, mapping_type, *data->mapping_q1_data, output); // otherwise just stick with what we already had MappingQ1::transform(input, mapping_type, mapping_data, output); @@ -913,7 +927,7 @@ transform (const VectorSlice > > input, // whether we should in fact work on the Q1 portion of it if (const InternalData *data = dynamic_cast(&mapping_data)) if (data->use_mapping_q1_on_current_cell) - return MappingQ1::transform (input, mapping_type, *data->mapping_q1_data, output); + return q1_mapping->transform (input, mapping_type, *data->mapping_q1_data, output); // otherwise just stick with what we already had MappingQ1::transform(input, mapping_type, mapping_data, output); @@ -938,7 +952,7 @@ transform (const VectorSlice // whether we should in fact work on the Q1 portion of it if (const InternalData *data = dynamic_cast(&mapping_data)) if (data->use_mapping_q1_on_current_cell) - return MappingQ1::transform (input, mapping_type, *data->mapping_q1_data, output); + return q1_mapping->transform (input, mapping_type, *data->mapping_q1_data, output); // otherwise just stick with what we already had MappingQ1::transform(input, mapping_type, mapping_data, output); @@ -961,7 +975,7 @@ transform (const VectorSlice > > input, // whether we should in fact work on the Q1 portion of it if (const InternalData *data = dynamic_cast(&mapping_data)) if (data->use_mapping_q1_on_current_cell) - return MappingQ1::transform (input, mapping_type, *data->mapping_q1_data, output); + return q1_mapping->transform (input, mapping_type, *data->mapping_q1_data, output); // otherwise just stick with what we already had MappingQ1::transform(input, mapping_type, mapping_data, output); diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc index 5fe49778e6..06a42dd06f 100644 --- a/source/fe/mapping_q_generic.cc +++ b/source/fe/mapping_q_generic.cc @@ -859,7 +859,7 @@ MappingQGeneric::get_data (const UpdateFlags update_flags, template -typename Mapping::InternalDataBase * +typename MappingQGeneric::InternalData * MappingQGeneric::get_face_data (const UpdateFlags update_flags, const Quadrature &quadrature) const { @@ -874,7 +874,7 @@ MappingQGeneric::get_face_data (const UpdateFlags update_fl template -typename Mapping::InternalDataBase * +typename MappingQGeneric::InternalData * MappingQGeneric::get_subface_data (const UpdateFlags update_flags, const Quadrature& quadrature) const { -- 2.39.5