From 0e474824a8edd6942120c6a80da892cbac26e5fc Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 8 Sep 2015 10:24:18 -0500 Subject: [PATCH] Move transform_unit_to_real_cell from MappingQ1 to MappingQGeneric. --- include/deal.II/fe/mapping_q.h | 14 ++++++- include/deal.II/fe/mapping_q1.h | 25 ------------ include/deal.II/fe/mapping_q_generic.h | 32 +++++++++++++++ source/fe/mapping_q.cc | 33 +++++----------- source/fe/mapping_q1.cc | 54 -------------------------- source/fe/mapping_q_generic.cc | 43 ++++++++++++++++++++ 6 files changed, 96 insertions(+), 105 deletions(-) diff --git a/include/deal.II/fe/mapping_q.h b/include/deal.II/fe/mapping_q.h index afe0414df2..96bce39faf 100644 --- a/include/deal.II/fe/mapping_q.h +++ b/include/deal.II/fe/mapping_q.h @@ -450,9 +450,19 @@ protected: const FE_Q feq; /** - * Pointer to a Q1 mapping to be used on interior cells unless + * 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. + * 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. */ std_cxx11::unique_ptr > q1_mapping; diff --git a/include/deal.II/fe/mapping_q1.h b/include/deal.II/fe/mapping_q1.h index 4ec4d7b928..fae3b12478 100644 --- a/include/deal.II/fe/mapping_q1.h +++ b/include/deal.II/fe/mapping_q1.h @@ -69,12 +69,6 @@ public: * @{ */ - // for documentation, see the Mapping base class - virtual - Point - transform_unit_to_real_cell (const typename Triangulation::cell_iterator &cell, - const Point &p) const; - // for documentation, see the Mapping base class virtual Point @@ -146,25 +140,6 @@ protected: transform_real_to_unit_cell_initial_guess (const std::vector > &vertex, const Point &p) const; - - /** - * Transforms a point @p p on the unit cell to the point @p p_real on the - * real cell @p cell and returns @p p_real. - * - * This function is called by @p transform_unit_to_real_cell and multiple - * times (through the Newton iteration) by @p - * transform_real_to_unit_cell_internal. - * - * Takes a reference to an @p InternalData that must already include the - * shape values at point @p p and the mapping support points of the cell. - * - * This @p InternalData argument avoids multiple computations of the shape - * values at point @p p and especially multiple computations of the mapping - * support points. - */ - Point - transform_unit_to_real_cell_internal (const InternalData &mdata) const; - /** * Transforms the point @p p on the real cell to the corresponding point on * the unit cell @p cell by a Newton iteration. diff --git a/include/deal.II/fe/mapping_q_generic.h b/include/deal.II/fe/mapping_q_generic.h index cb5df7fd60..5eae7b0f9e 100644 --- a/include/deal.II/fe/mapping_q_generic.h +++ b/include/deal.II/fe/mapping_q_generic.h @@ -78,6 +78,20 @@ public: virtual bool preserves_vertex_locations () const; + /** + * @name Mapping points between reference and real cells + * @{ + */ + + // for documentation, see the Mapping base class + virtual + Point + transform_unit_to_real_cell (const typename Triangulation::cell_iterator &cell, + const Point &p) const; + + /** + * @} + */ /** * @name Functions to transform tensors from reference to real coordinates @@ -454,6 +468,24 @@ protected: compute_mapping_support_points (const typename Triangulation::cell_iterator &cell, std::vector > &a) const = 0; + /** + * Transforms a point @p p on the unit cell to the point @p p_real on the + * real cell @p cell and returns @p p_real. + * + * This function is called by @p transform_unit_to_real_cell and multiple + * times (through the Newton iteration) by @p + * transform_real_to_unit_cell_internal. + * + * Takes a reference to an @p InternalData that must already include the + * shape values at point @p p and the mapping support points of the cell. + * + * This @p InternalData argument avoids multiple computations of the shape + * values at point @p p and especially multiple computations of the mapping + * support points. + */ + Point + transform_unit_to_real_cell_internal (const InternalData &mdata) const; + /** * Make MappingQ a friend since it needs to call the * fill_fe_values() functions on its MappingQ1 sub-object. diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc index 4447272c2b..13443f4338 100644 --- a/source/fe/mapping_q.cc +++ b/source/fe/mapping_q.cc @@ -75,11 +75,9 @@ MappingQ::MappingQ (const unsigned int degree, || (dim != spacedim)), feq(degree), - q1_mapping (this->use_mapping_q_on_all_cells - ? - 0 - : - new MappingQ1()), + // 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 MappingQ1()), line_support_points(degree+1) { Assert(n_inner+n_outer==Utilities::fixed_power(degree+1), @@ -105,11 +103,9 @@ MappingQ::MappingQ (const MappingQ &mapping) 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), + // 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())), line_support_points(mapping.line_support_points) { Assert(n_inner+n_outer==Utilities::fixed_power(this->polynomial_degree+1), @@ -1001,19 +997,9 @@ 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()) - { - const Quadrature point_quadrature(p); - std_cxx11::unique_ptr mdata (get_data(update_quadrature_points, - point_quadrature)); - - compute_mapping_support_points(cell, mdata->mapping_support_points); - - return this->transform_unit_to_real_cell_internal(*mdata); - } + return this->MappingQGeneric::transform_unit_to_real_cell (cell, p); else - { - return q1_mapping->transform_unit_to_real_cell (cell, p); - } + return q1_mapping->transform_unit_to_real_cell (cell, p); } @@ -1036,8 +1022,7 @@ transform_real_to_unit_cell (const typename Triangulation::cell_it Point initial_p_unit; try { - initial_p_unit - = StaticMappingQ1::mapping.transform_real_to_unit_cell(cell, p); + initial_p_unit = q1_mapping->transform_real_to_unit_cell(cell, p); } catch (const typename Mapping::ExcTransformationFailed &) { diff --git a/source/fe/mapping_q1.cc b/source/fe/mapping_q1.cc index e3e1777a16..f52b098027 100644 --- a/source/fe/mapping_q1.cc +++ b/source/fe/mapping_q1.cc @@ -205,60 +205,6 @@ MappingQ1::compute_mapping_support_points( - -template -Point -MappingQ1::transform_unit_to_real_cell ( - const typename Triangulation::cell_iterator &cell, - const Point &p) const -{ - const Quadrature point_quadrature(p); - - //TODO: Use get_data() here once MappingQ is no longer derived from - //MappingQ1. this doesn't currently work because we here really need - //a Q1 InternalData, but MappingQGeneric produces one with the - //polynomial degree of the MappingQ - std_cxx11::unique_ptr mdata (new InternalData(1)); - mdata->initialize (this->requires_update_flags (update_quadrature_points), - point_quadrature, 1); - - // compute the mapping support - // points - compute_mapping_support_points(cell, mdata->mapping_support_points); - - // Mapping support points can be - // bigger than necessary. If this - // is the case, force them to be - // Q1. - if (mdata->mapping_support_points.size() > mdata->shape_values.size()) - mdata->mapping_support_points.resize - (GeometryInfo::vertices_per_cell); - - - return transform_unit_to_real_cell_internal(*mdata); -} - - - -template -Point -MappingQ1:: -transform_unit_to_real_cell_internal (const InternalData &data) const -{ - AssertDimension (data.shape_values.size(), - data.mapping_support_points.size()); - - // use now the InternalData to - // compute the point in real space. - Point p_real; - for (unsigned int i=0; i::get_degree() const +template +Point +MappingQGeneric:: +transform_unit_to_real_cell (const typename Triangulation::cell_iterator &cell, + const Point &p) const +{ + //TODO: This function is inefficient -- it might as well evaluate + //the shape functions directly, rather than first building the + //InternalData object and then working with it + + const Quadrature point_quadrature(p); + std_cxx11::unique_ptr mdata (new InternalData(polynomial_degree)); + mdata->initialize (this->requires_update_flags(update_quadrature_points), + point_quadrature, + 1); + + // compute the mapping support + // points + compute_mapping_support_points(cell, mdata->mapping_support_points); + return transform_unit_to_real_cell_internal(*mdata); +} + + + +template +Point +MappingQGeneric:: +transform_unit_to_real_cell_internal (const InternalData &data) const +{ + AssertDimension (data.shape_values.size(), + data.mapping_support_points.size()); + + // use now the InternalData to + // compute the point in real space. + Point p_real; + for (unsigned int i=0; i UpdateFlags MappingQGeneric::requires_update_flags (const UpdateFlags in) const -- 2.39.5