From e89706980cce0586d144c9335285eb5fb1534e8a Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Sat, 5 Sep 2015 08:07:33 -0500 Subject: [PATCH] Store the correct polynomial degree in MappingQGeneric. Currently, MappingQ derives from MappingQ1 which passes 1 down as polynomial degree. This worked 'by accident' because the polynomial degree is only used by MappingQGeneric in get_data() and friends, which are overloaded by MappingQ. However, the correct place to store the polynomial degree is clearly in MappingQGeneric. --- include/deal.II/fe/mapping_q.h | 12 ---- include/deal.II/fe/mapping_q1.h | 11 +++- include/deal.II/fe/mapping_q_generic.h | 6 ++ source/fe/mapping_q.cc | 87 ++++++++++++-------------- source/fe/mapping_q1.cc | 33 ++++++---- source/fe/mapping_q_generic.cc | 9 +++ 6 files changed, 82 insertions(+), 76 deletions(-) diff --git a/include/deal.II/fe/mapping_q.h b/include/deal.II/fe/mapping_q.h index 7f11705c25..b954602e3d 100644 --- a/include/deal.II/fe/mapping_q.h +++ b/include/deal.II/fe/mapping_q.h @@ -146,12 +146,6 @@ public: const typename Mapping::InternalDataBase &internal, VectorSlice > > output) const; - /** - * Return the degree of the mapping, i.e. the value which was passed to the - * constructor. - */ - unsigned int get_degree () const; - /** * Return a pointer to a copy of the present object. The caller of this copy * then assumes ownership of it. @@ -421,12 +415,6 @@ protected: int, << "laplace_vector not set for degree=" << arg1 << "."); - /** - * Degree @p p of the polynomials used as shape functions for the Qp mapping - * of cells at the boundary. - */ - const unsigned int degree; - /** * Number of inner mapping shape functions. */ diff --git a/include/deal.II/fe/mapping_q1.h b/include/deal.II/fe/mapping_q1.h index 92943025ac..4ec4d7b928 100644 --- a/include/deal.II/fe/mapping_q1.h +++ b/include/deal.II/fe/mapping_q1.h @@ -98,14 +98,19 @@ public: */ typedef typename MappingQGeneric::InternalData InternalData; -protected: - - /** * @} */ protected: + + /** + * Constructor. This constructor is for odd purposes: MappingQ is + * derived from this class (for historical reasons) and it needs a + * way to pass down the "true" polynomial degree of the mapping. + */ + MappingQ1 (const unsigned int degree); + /* Trick to templatize transform_real_to_unit_cell */ template Point diff --git a/include/deal.II/fe/mapping_q_generic.h b/include/deal.II/fe/mapping_q_generic.h index 5a1e97b2da..3791c7da33 100644 --- a/include/deal.II/fe/mapping_q_generic.h +++ b/include/deal.II/fe/mapping_q_generic.h @@ -63,6 +63,12 @@ public: */ MappingQGeneric (const unsigned int polynomial_degree); + /** + * 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 MappingQ1 preserves vertex locations. */ diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc index 6635f30f90..5eb71043b0 100644 --- a/source/fe/mapping_q.cc +++ b/source/fe/mapping_q.cc @@ -54,10 +54,10 @@ MappingQ::InternalData::memory_consumption () const template -MappingQ::MappingQ (const unsigned int p, +MappingQ::MappingQ (const unsigned int degree, const bool use_mapping_q_on_all_cells) : - degree(p), + MappingQ1(degree), n_inner(Utilities::fixed_power(degree-1)), n_outer((dim==1) ? 2 : ((dim==2) ? @@ -88,7 +88,7 @@ typename MappingQ::InternalData * MappingQ::get_data (const UpdateFlags update_flags, const Quadrature &quadrature) const { - InternalData *data = new InternalData(degree); + 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; @@ -120,7 +120,7 @@ typename Mapping::InternalDataBase * MappingQ::get_face_data (const UpdateFlags update_flags, const Quadrature& quadrature) const { - InternalData *data = new InternalData(degree); + 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 @@ -153,7 +153,7 @@ typename Mapping::InternalDataBase * MappingQ::get_subface_data (const UpdateFlags update_flags, const Quadrature& quadrature) const { - InternalData *data = new InternalData(degree); + 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 @@ -219,7 +219,7 @@ fill_fe_values (const typename Triangulation::cell_iterator &cell, const CellSimilarity::Similarity updated_cell_similarity = ((data.use_mapping_q1_on_current_cell == false) && - (get_degree() > 1) + (this->polynomial_degree > 1) ? CellSimilarity::invalid_next_cell : @@ -333,14 +333,14 @@ template void MappingQ::set_laplace_on_quad_vector(Table<2,double> &loqvs) const { - Assert(degree>1, ExcInternalError()); - const unsigned int n_inner_2d=(degree-1)*(degree-1); - const unsigned int n_outer_2d=4+4*(degree-1); + Assert(this->polynomial_degree>1, ExcInternalError()); + const unsigned int n_inner_2d=(this->polynomial_degree-1)*(this->polynomial_degree-1); + const unsigned int n_outer_2d=4+4*(this->polynomial_degree-1); // first check whether we have precomputed the values for some polynomial // degree; the sizes of arrays is n_inner_2d*n_outer_2d double const *loqv_ptr=0; - switch (degree) + switch (this->polynomial_degree) { // for degree==1, we shouldn't have to compute any support points, since // all of them are on the vertices @@ -380,7 +380,7 @@ MappingQ::set_laplace_on_quad_vector(Table<2,double> &loqvs) const compute_laplace_vector(loqvs); else if (dim == 3) { - MappingQ<2,2> mapping_2d(this->degree); + MappingQ<2,2> mapping_2d(this->polynomial_degree); loqvs = mapping_2d.laplace_on_quad_vector; } } @@ -389,7 +389,7 @@ MappingQ::set_laplace_on_quad_vector(Table<2,double> &loqvs) const // this for (unsigned int unit_point=0; unit_pointdegree, + loqvs[unit_point].end(),0.)-1)<1e-13*this->polynomial_degree, ExcInternalError()); } @@ -399,12 +399,12 @@ template <> void MappingQ<3>::set_laplace_on_hex_vector(Table<2,double> &lohvs) const { - Assert(degree>1, ExcInternalError()); + Assert(this->polynomial_degree>1, ExcInternalError()); // first check whether we have precomputed the values for some polynomial // degree double const *lohv_ptr=0; - if (degree==2) + if (this->polynomial_degree==2) { static const double loqv2[26] = {1/128., 1/128., 1/128., 1/128., 1/128., 1/128., 1/128., 1/128., @@ -432,7 +432,7 @@ MappingQ<3>::set_laplace_on_hex_vector(Table<2,double> &lohvs) const // this for (unsigned int unit_point=0; unit_pointdegree, + lohvs[unit_point].end(),0.) - 1)<1e-13*this->polynomial_degree, ExcInternalError()); } @@ -465,13 +465,13 @@ MappingQ::compute_laplace_vector(Table<2,double> &lvs) const // for degree==1, we shouldn't have to compute any support points, since all // of them are on the vertices - Assert(degree>1, ExcInternalError()); + Assert(this->polynomial_degree>1, ExcInternalError()); // compute the shape gradients at the quadrature points on the unit cell - const QGauss quadrature(degree+1); + const QGauss quadrature(this->polynomial_degree+1); const unsigned int n_q_points=quadrature.size(); - InternalData quadrature_data(degree); + InternalData quadrature_data(this->polynomial_degree); quadrature_data.shape_derivatives.resize(quadrature_data.n_shape_functions * n_q_points); quadrature_data.compute_shape_function_values(quadrature.get_points()); @@ -532,10 +532,10 @@ MappingQ::apply_laplace_vector(const Table<2,double> &lvs, // compute the quad laplace vector. this is mentioned in the constructor of // this class, although I don't understand the reason for not aborting there // any more [WB] - Assert(lvs.n_rows()!=0, ExcLaplaceVectorNotSet(degree)); + Assert(lvs.n_rows()!=0, ExcLaplaceVectorNotSet(this->polynomial_degree)); const unsigned int n_inner_apply=lvs.n_rows(); - Assert(n_inner_apply==n_inner || n_inner_apply==(degree-1)*(degree-1), + Assert(n_inner_apply==n_inner || n_inner_apply==(this->polynomial_degree-1)*(this->polynomial_degree-1), ExcInternalError()); const unsigned int n_outer_apply=lvs.n_cols(); Assert(a.size()==n_outer_apply, @@ -588,7 +588,7 @@ MappingQ::compute_support_points_laplace(const typename Triangulat for (unsigned int i=0; i::vertices_per_cell; ++i) a[i] = cell->vertex(i); - if (degree>1) + if (this->polynomial_degree>1) switch (dim) { case 1: @@ -625,7 +625,7 @@ MappingQ::add_line_support_points (const typename Triangulation > &a) const { // if we only need the midpoint, then ask for it. - if (degree==2) + if (this->polynomial_degree==2) { for (unsigned int line_no=0; line_no::lines_per_cell; ++line_no) { @@ -645,7 +645,7 @@ MappingQ::add_line_support_points (const typename Triangulation > line_points (degree-1); + std::vector > line_points (this->polynomial_degree-1); // loop over each of the lines, and if it is at the boundary, then first // get the boundary description and second compute the points on it for (unsigned int line_no=0; line_no::lines_per_cell; ++line_no) @@ -699,9 +699,9 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell, static const StraightBoundary<3> straight_boundary; // used if face quad at boundary or entirely in the interior of the domain - std::vector > quad_points ((degree-1)*(degree-1)); + std::vector > quad_points ((this->polynomial_degree-1)*(this->polynomial_degree-1)); // used if only one line of face quad is at boundary - std::vector > b(4*degree); + std::vector > b(4*this->polynomial_degree); // Used by the new Manifold interface. This vector collects the // vertices used to compute the intermediate points. @@ -769,12 +769,12 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell, // call of function apply_laplace_vector increases size of b // about 1. There resize b for the case the mentioned function // was already called. - b.resize(4*degree); + b.resize(4*this->polynomial_degree); // b is of size 4*degree, make sure that this is the right size - Assert(b.size()==vertices_per_face+lines_per_face*(degree-1), + Assert(b.size()==vertices_per_face+lines_per_face*(this->polynomial_degree-1), ExcDimensionMismatch(b.size(), - vertices_per_face+lines_per_face*(degree-1))); + vertices_per_face+lines_per_face*(this->polynomial_degree-1))); // sort the points into b. We used access from the cell (not // from the face) to fill b, so we can assume a standard face @@ -784,19 +784,19 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell, b[i]=a[GeometryInfo<3>::face_to_cell_vertices(face_no, i)]; for (unsigned int i=0; ipolynomial_degree-1; ++j) + b[vertices_per_face+i*(this->polynomial_degree-1)+j]= a[vertices_per_cell + GeometryInfo<3>::face_to_cell_lines( - face_no, i)*(degree-1)+j]; + face_no, i)*(this->polynomial_degree-1)+j]; // Now b includes the support points on the quad and we can // apply the laplace vector apply_laplace_vector(laplace_on_quad_vector, b); - Assert(b.size()==4*degree+(degree-1)*(degree-1), - ExcDimensionMismatch(b.size(), 4*degree+(degree-1)*(degree-1))); + Assert(b.size()==4*this->polynomial_degree+(this->polynomial_degree-1)*(this->polynomial_degree-1), + ExcDimensionMismatch(b.size(), 4*this->polynomial_degree+(this->polynomial_degree-1)*(this->polynomial_degree-1))); - for (unsigned int i=0; i<(degree-1)*(degree-1); ++i) - a.push_back(b[4*degree+i]); + for (unsigned int i=0; i<(this->polynomial_degree-1)*(this->polynomial_degree-1); ++i) + a.push_back(b[4*this->polynomial_degree+i]); } else { @@ -829,7 +829,7 @@ MappingQ<2,3>:: add_quad_support_points(const Triangulation<2,3>::cell_iterator &cell, std::vector > &a) const { - std::vector > quad_points ((degree-1)*(degree-1)); + std::vector > quad_points ((this->polynomial_degree-1)*(this->polynomial_degree-1)); get_intermediate_points_on_object (cell->get_manifold(), cell, quad_points); for (unsigned int i=0; i::cell_it -template -unsigned int -MappingQ::get_degree() const -{ - return degree; -} - - - template Mapping * MappingQ::clone () const { - return new MappingQ(degree); + return new MappingQ(this->polynomial_degree); } @@ -1117,7 +1108,7 @@ MappingQ::get_intermediate_points (const Manifold & { // If two points are passed, these are the two vertices, and // we can only compute degree-1 intermediate points. - AssertDimension(n, degree-1); + AssertDimension(n, this->polynomial_degree-1); for (unsigned int i=0; i::get_intermediate_points (const Manifold & // If four points are passed, these are the two vertices, and // we can only compute (degree-1)*(degree-1) intermediate // points. - AssertDimension(m, degree-1); + AssertDimension(m, this->polynomial_degree-1); for (unsigned int i=0; i::MappingQ1 () +template +MappingQ1::MappingQ1 (const unsigned int p) + : + MappingQGeneric (p) +{} + + + namespace internal { namespace MappingQ1 @@ -197,13 +205,14 @@ MappingQ1::transform_unit_to_real_cell ( const typename Triangulation::cell_iterator &cell, const Point &p) const { - // Use the get_data function to create an InternalData with data - // vectors of the right size and transformation shape values already - // computed at point p. const Quadrature point_quadrature(p); - std_cxx11::unique_ptr mdata (get_data(update_quadrature_points | update_jacobians, - point_quadrature)); + //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 | update_jacobians), point_quadrature, 1); // compute the mapping support // points @@ -476,20 +485,18 @@ transform_real_to_unit_cell (const typename Triangulation::cell_it // it shouldn't really make any difference... // initial_p_unit = GeometryInfo::project_to_unit_cell(initial_p_unit); - // Use the get_data function to - // create an InternalData with data - // vectors of the right size and - // transformation shape values and - // derivatives already computed at - // point initial_p_unit. const Quadrature point_quadrature(initial_p_unit); UpdateFlags update_flags = update_quadrature_points | update_jacobians; if (spacedim>dim) update_flags |= update_jacobian_grads; - std_cxx11::unique_ptr mdata(MappingQ1::get_data(update_flags, - point_quadrature)); + //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_flags), point_quadrature, 1); compute_mapping_support_points (cell, mdata->mapping_support_points); // The support points have to be at diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc index 3d8df89ddf..53534ffb3b 100644 --- a/source/fe/mapping_q_generic.cc +++ b/source/fe/mapping_q_generic.cc @@ -776,6 +776,15 @@ MappingQGeneric::MappingQGeneric (const unsigned int p) +template +unsigned int +MappingQGeneric::get_degree() const +{ + return polynomial_degree; +} + + + template UpdateFlags MappingQGeneric::requires_update_flags (const UpdateFlags in) const -- 2.39.5