From 4db2a411eafd5991cf08e09a71ce6d74c401ce79 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Sat, 12 Sep 2015 07:07:07 -0500 Subject: [PATCH] Make several more functions part of the internal implementation. This allows us to remove them from the public interface. --- include/deal.II/fe/mapping_q.h | 109 +++----- source/fe/mapping_q.cc | 439 +++++++++++++++++---------------- 2 files changed, 258 insertions(+), 290 deletions(-) diff --git a/include/deal.II/fe/mapping_q.h b/include/deal.II/fe/mapping_q.h index 19af543546..72e8678c09 100644 --- a/include/deal.II/fe/mapping_q.h +++ b/include/deal.II/fe/mapping_q.h @@ -266,48 +266,22 @@ protected: protected: - /** - * This function is needed by the constructor of - * MappingQ for dim= 2 and 3. - * - * For degree<4 this function sets the @p support_point_weights_on_quad to - * the hardcoded data. For degree>=4 and MappingQ<2> this vector is - * computed. - * - * For the definition of the @p support_point_weights_on_quad please refer to - * equation (8) of the `mapping' report. - */ - void - set_support_point_weights_on_quad(Table<2,double> &loqvs) const; - - /** - * This function is needed by the constructor of MappingQ<3>. - * - * For degree==2 this function sets the @p support_point_weights_on_hex to - * the hardcoded data. For degree>2 this vector is computed. - * - * For the definition of the @p support_point_weights_on_hex please refer to - * equation (8) of the `mapping' report. - */ - void set_support_point_weights_on_hex(Table<2,double> &lohvs) const; - - /** - * Compute the support_point_weights_on_quad(hex)_vector. - * - * Called by the set_support_point_weights_on_quad(hex)_vector functions if the - * data is not yet hardcoded. - * - * For the definition of the support_point_weights_on_quad(hex)_vector please - * refer to equation (8) of the `mapping' report. - */ - void compute_laplace_vector(Table<2,double> &lvs) const; - /** * Compute the support points of the mapping. Interior support * points (ie. support points in quads for 2d, in hexes for 3d) are * computed using the solution of a Laplace equation with the * position of the outer support points as boundary values, in order * to make the transformation as smooth as possible. + * + * The function works its way from the vertices (which it takes from + * the given cell) via the support points on the line (for which it + * calls the add_line_support_points() function) and the support + * points on the quad faces (in 3d, for which it calls the + * add_quad_support_points() function). It then adds interior + * support points that are either computed by interpolation from the + * surrounding points using weights computed by solving a Laplace + * equation, or if dimdim=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 not on vertices are not included. - * - * Needed by the @p compute_support_points_laplace function . For - * dim=1 this function is empty. - * - * This function is made virtual in order to allow derived classes to choose - * shape function support points differently than the present class, which - * chooses the points as interpolation points on the boundary. - */ + * For dim=2,3. Append the support points of all shape + * functions located on bounding lines of the given cell to the + * vector @p a. Points located on the vertices of a line are not + * included. + * + * Needed by the @p compute_support_points() function. For + * dim=1 this function is empty. The function uses the + * underlying manifold object of the line (or, if none is set, of + * the cell) for the location of the requested points. + * + * This function is made virtual in order to allow derived classes + * to choose shape function support points differently than the + * present class, which chooses the points as interpolation points + * on the boundary. + */ 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 quad but not on vertices are not included. + * For dim=3. Append the support points of all shape + * functions located on bounding faces (quads in 3d) of the given + * cell to the vector @p a. Points located on the vertices or lines + * of a quad are not included. * - * Needed by the @p compute_support_points_laplace function. For - * dim=1 and dim=2 this function is empty. + * Needed by the @p compute_support_points() function. For + * dim=1 and dim=2 this function is empty. The + * function uses the underlying manifold object of the quad (or, if + * none is set, of the cell) for the location of the requested + * points. * - * This function is made virtual in order to allow derived classes to choose - * shape function support points differently than the present class, which - * chooses the points as interpolation points on the boundary. + * This function is made virtual in order to allow derived classes + * to choose shape function support points differently than the + * present class, which chooses the points as interpolation points + * on the boundary. */ virtual void @@ -439,22 +422,6 @@ protected: /*@}*/ -/* -------------- declaration of explicit specializations ------------- */ - -#ifndef DOXYGEN - -template <> -void MappingQ<1>::set_support_point_weights_on_quad(Table<2,double> &) const; - -template <> -void MappingQ<3>::set_support_point_weights_on_hex(Table<2,double> &lohvs) const; - -template <> -void MappingQ<1>::compute_laplace_vector(Table<2,double> &) const; - - -#endif // DOXYGEN - DEAL_II_NAMESPACE_CLOSE #endif diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc index 9db5d010b2..62b02f8197 100644 --- a/source/fe/mapping_q.cc +++ b/source/fe/mapping_q.cc @@ -53,6 +53,222 @@ MappingQ::InternalData::memory_consumption () const +namespace +{ + /** + * Compute the support_point_weights_on_quad(hex) arrays. + * + * Called by the compute_support_point_weights_on_quad(hex) functions if the + * data is not yet hardcoded. + * + * For the definition of the support_point_weights_on_quad(hex) please + * refer to equation (8) of the `mapping' report. + */ + template + Table<2,double> + compute_laplace_vector(const unsigned int polynomial_degree) + { + Table<2,double> lvs; + + Assert(lvs.n_rows()==0, ExcInternalError()); + Assert(dim==2 || dim==3, ExcNotImplemented()); + + // for degree==1, we shouldn't have to compute any support points, since all + // of them are on the vertices + Assert(polynomial_degree>1, ExcInternalError()); + + const unsigned int n_inner = Utilities::fixed_power(polynomial_degree-1); + const unsigned int n_outer = (dim==1) ? 2 : + ((dim==2) ? + 4+4*(polynomial_degree-1) : + 8+12*(polynomial_degree-1)+6*(polynomial_degree-1)*(polynomial_degree-1)); + + + // compute the shape gradients at the quadrature points on the unit cell + const QGauss quadrature(polynomial_degree+1); + const unsigned int n_q_points=quadrature.size(); + + typename MappingQGeneric::InternalData quadrature_data(polynomial_degree); + quadrature_data.shape_derivatives.resize(quadrature_data.n_shape_functions * + n_q_points); + quadrature_data.compute_shape_function_values(quadrature.get_points()); + + // Compute the stiffness matrix of the inner dofs + FullMatrix S(n_inner); + for (unsigned int point=0; point T(n_inner, n_outer); + for (unsigned int point=0; point S_1(n_inner); + S_1.invert(S); + + FullMatrix S_1_T(n_inner, n_outer); + + // S:=S_1*T + S_1.mmult(S_1_T,T); + + // Resize and initialize the lvs + lvs.reinit (n_inner, n_outer); + for (unsigned int i=0; iMappingQ for dim= 2 and 3. + * + * For degree<4 this function sets the @p support_point_weights_on_quad to + * the hardcoded data. For degree>=4 and MappingQ<2> this vector is + * computed. + * + * For the definition of the @p support_point_weights_on_quad please refer to + * equation (8) of the `mapping' report. + */ + template + Table<2,double> + compute_support_point_weights_on_quad(const unsigned int polynomial_degree) + { + Table<2,double> loqvs; + + // in 1d, there are no quads, so return an empty object + if (dim == 1) + return loqvs; + + // we are asked to compute weights for interior support points, but + // there are no interior points if degree==1 + if (polynomial_degree == 1) + return loqvs; + + const unsigned int n_inner_2d=(polynomial_degree-1)*(polynomial_degree-1); + const unsigned int n_outer_2d=4+4*(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 + if (polynomial_degree == 2) + { + // (checked these values against the output of compute_laplace_vector + // again, and found they're indeed right -- just in case someone wonders + // where they come from -- WB) + static const double loqv2[1*8] + = {1/16., 1/16., 1/16., 1/16., 3/16., 3/16., 3/16., 3/16.}; + Assert (sizeof(loqv2)/sizeof(loqv2[0]) == + n_inner_2d * n_outer_2d, + ExcInternalError()); + + // copy and return + loqvs.reinit(n_inner_2d, n_outer_2d); + for (unsigned int unit_point=0; unit_point(polynomial_degree); + } + + // the sum of weights of the points at the outer rim should be one. check + // this + for (unsigned int unit_point=0; unit_pointMappingQ<3>. + * + * For degree==2 this function sets the @p support_point_weights_on_hex to + * the hardcoded data. For degree>2 this vector is computed. + * + * For the definition of the @p support_point_weights_on_hex please refer to + * equation (8) of the `mapping' report. + */ + template + Table<2,double> + compute_support_point_weights_on_hex(const unsigned int polynomial_degree) + { + Table<2,double> lohvs; + + // in 1d and 2d, there are no hexes, so return an empty object + if (dim < 3) + return lohvs; + + // we are asked to compute weights for interior support points, but + // there are no interior points if degree==1 + if (polynomial_degree == 1) + return lohvs; + + const unsigned int n_inner = Utilities::fixed_power(polynomial_degree-1); + const unsigned int n_outer = 8+12*(polynomial_degree-1)+6*(polynomial_degree-1)*(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 + if (polynomial_degree == 2) + { + static const double lohv2[26] + = {1/128., 1/128., 1/128., 1/128., 1/128., 1/128., 1/128., 1/128., + 7/192., 7/192., 7/192., 7/192., 7/192., 7/192., 7/192., 7/192., + 7/192., 7/192., 7/192., 7/192., + 1/12., 1/12., 1/12., 1/12., 1/12., 1/12. + }; + + // copy and return + lohvs.reinit(n_inner, n_outer); + for (unsigned int unit_point=0; unit_point(polynomial_degree); + } + + // the sum of weights of the points at the outer rim should be one. check + // this + for (unsigned int unit_point=0; unit_point MappingQ::MappingQ (const unsigned int degree, const bool use_mapping_q_on_all_cells) @@ -83,14 +299,8 @@ MappingQ::MappingQ (const unsigned int degree, Assert(n_inner+n_outer==Utilities::fixed_power(degree+1), ExcInternalError()); - // build support_point_weights_on_quad - if (degree>1) - { - if (dim >= 2) - set_support_point_weights_on_quad(support_point_weights_on_quad); - if (dim >= 3) - set_support_point_weights_on_hex(support_point_weights_on_hex); - } + 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); } @@ -111,14 +321,8 @@ MappingQ::MappingQ (const MappingQ &mapping) Assert(n_inner+n_outer==Utilities::fixed_power(this->polynomial_degree+1), ExcInternalError()); - // build support_point_weights_on_quad - if (this->polynomial_degree>1) - { - if (dim >= 2) - set_support_point_weights_on_quad(support_point_weights_on_quad); - if (dim >= 3) - set_support_point_weights_on_hex(support_point_weights_on_hex); - } + 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); } @@ -338,209 +542,6 @@ fill_fe_subface_values (const typename Triangulation::cell_iterato } - -template <> -void -MappingQ<1>::set_support_point_weights_on_quad(Table<2,double> &) const -{ - Assert(false, ExcInternalError()); -} - - - -template -void -MappingQ::set_support_point_weights_on_quad(Table<2,double> &loqvs) const -{ - 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 (this->polynomial_degree) - { - // for degree==1, we shouldn't have to compute any support points, since - // all of them are on the vertices - - case 2: - { - // (checked these values against the output of compute_laplace_vector - // again, and found they're indeed right -- just in case someone wonders - // where they come from -- WB) - static const double loqv2[1*8] - = {1/16., 1/16., 1/16., 1/16., 3/16., 3/16., 3/16., 3/16.}; - loqv_ptr=&loqv2[0]; - Assert (sizeof(loqv2)/sizeof(loqv2[0]) == - n_inner_2d * n_outer_2d, - ExcInternalError()); - - break; - } - - // no other cases implemented, so simply fall through - default: - break; - } - - if (loqv_ptr!=0) - { - // precomputed. copy values to the loqvs array - loqvs.reinit(n_inner_2d, n_outer_2d); - for (unsigned int unit_point=0; unit_point mapping_2d(this->polynomial_degree); - loqvs = mapping_2d.support_point_weights_on_quad; - } - } - - // the sum of weights of the points at the outer rim should be one. check - // this - for (unsigned int unit_point=0; unit_pointpolynomial_degree, - ExcInternalError()); -} - - - -template <> -void -MappingQ<3>::set_support_point_weights_on_hex(Table<2,double> &lohvs) const -{ - Assert(this->polynomial_degree>1, ExcInternalError()); - - // first check whether we have precomputed the values for some polynomial - // degree - double const *lohv_ptr=0; - 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., - 7/192., 7/192., 7/192., 7/192., 7/192., 7/192., 7/192., 7/192., - 7/192., 7/192., 7/192., 7/192., - 1/12., 1/12., 1/12., 1/12., 1/12., 1/12. - }; - - lohv_ptr=&loqv2[0]; - } - - if (lohv_ptr!=0) - { - // precomputed. copy values to the lohvs array - lohvs.reinit(n_inner, n_outer); - for (unsigned int unit_point=0; unit_pointpolynomial_degree, - ExcInternalError()); -} - - - -template -void -MappingQ::set_support_point_weights_on_hex(Table<2,double> &) const -{ - Assert(false, ExcInternalError()); -} - - - -template <> -void -MappingQ<1>::compute_laplace_vector(Table<2,double> &) const -{ - Assert(false, ExcInternalError()); -} - - - -template -void -MappingQ::compute_laplace_vector(Table<2,double> &lvs) const -{ - Assert(lvs.n_rows()==0, ExcInternalError()); - Assert(dim==2 || dim==3, ExcNotImplemented()); - - // for degree==1, we shouldn't have to compute any support points, since all - // of them are on the vertices - Assert(this->polynomial_degree>1, ExcInternalError()); - - // compute the shape gradients at the quadrature points on the unit cell - const QGauss quadrature(this->polynomial_degree+1); - const unsigned int n_q_points=quadrature.size(); - - 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()); - - // Compute the stiffness matrix of the inner dofs - FullMatrix S(n_inner); - for (unsigned int point=0; point T(n_inner, n_outer); - for (unsigned int point=0; point S_1(n_inner); - S_1.invert(S); - - FullMatrix S_1_T(n_inner, n_outer); - - // S:=S_1*T - S_1.mmult(S_1_T,T); - - // Resize and initialize the lvs - lvs.reinit (n_inner, n_outer); - for (unsigned int i=0; i