From 2478b0eec5cc2f273440e0cb0d273f9d106ad4a7 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Fri, 11 Sep 2015 20:07:14 -0500 Subject: [PATCH] Move several functions from the public interface of MappingQ into the .cc file. --- include/deal.II/fe/mapping_q.h | 206 +++++---------- source/fe/mapping_q.cc | 444 ++++++++++++++++++--------------- 2 files changed, 305 insertions(+), 345 deletions(-) diff --git a/include/deal.II/fe/mapping_q.h b/include/deal.II/fe/mapping_q.h index 96bce39faf..19af543546 100644 --- a/include/deal.II/fe/mapping_q.h +++ b/include/deal.II/fe/mapping_q.h @@ -163,68 +163,6 @@ 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 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. - */ - 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. - * - * Needed by the @p compute_support_points_laplace function. For - * dim=1 and dim=2 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. - */ - virtual void - add_quad_support_points(const typename Triangulation::cell_iterator &cell, - std::vector > &a) const; - - -private: - /** - * Ask the manifold descriptor to return intermediate points on lines or - * faces. The function needs to return one or multiple points (depending on - * the number of elements in the output vector @p points that lie inside a - * line, quad or hex). Whether it is a line, quad or hex doesn't really - * matter to this function but it can be inferred from the number of input - * points in the @p surrounding_points vector. - */ - void get_intermediate_points(const Manifold &manifold, - const std::vector > &surrounding_points, - std::vector > &points) const; - - - /** - * Ask the manifold descriptor to return intermediate points on the object - * pointed to by the TriaIterator @p iter. This function tries to be - * backward compatible with respect to the differences between - * Boundary and Manifold, querying the first - * whenever the passed @p manifold can be upgraded to a - * Boundary. - */ - template - void get_intermediate_points_on_object(const Manifold &manifold, - const TriaIterator &iter, - std::vector > &points) const; - /** * @name Interface with FEValues * @{ @@ -332,90 +270,106 @@ protected: * This function is needed by the constructor of * MappingQ for dim= 2 and 3. * - * For degree<4 this function sets the @p laplace_on_quad_vector to + * 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 laplace_on_quad_vector please refer to + * For the definition of the @p support_point_weights_on_quad please refer to * equation (8) of the `mapping' report. */ void - set_laplace_on_quad_vector(Table<2,double> &loqvs) const; + 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 laplace_on_hex_vector to + * 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 laplace_on_hex_vector please refer to + * For the definition of the @p support_point_weights_on_hex please refer to * equation (8) of the `mapping' report. */ - void set_laplace_on_hex_vector(Table<2,double> &lohvs) const; + void set_support_point_weights_on_hex(Table<2,double> &lohvs) const; /** - * Computes the laplace_on_quad(hex)_vector. + * Compute the support_point_weights_on_quad(hex)_vector. * - * Called by the set_laplace_on_quad(hex)_vector functions if the + * Called by the set_support_point_weights_on_quad(hex)_vector functions if the * data is not yet hardcoded. * - * For the definition of the laplace_on_quad(hex)_vector please + * 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; /** - * Takes a laplace_on_hex(quad)_vector and applies it to the vector - * @p a to compute the inner support points as a linear combination of the - * exterior points. - * - * The vector @p a initially contains the locations of the @p n_outer - * points, the @p n_inner computed inner points are appended. - * - * See equation (7) of the `mapping' report. + * 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. */ - void apply_laplace_vector(const Table<2,double> &lvs, - std::vector > &a) const; + virtual + void + compute_mapping_support_points(const typename Triangulation::cell_iterator &cell, + std::vector > &a) const; + /** - * Computes the support points of the mapping. - */ - virtual void compute_mapping_support_points( - const typename Triangulation::cell_iterator &cell, - std::vector > &a) const; + * 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 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. + */ + virtual + void + add_line_support_points (const typename Triangulation::cell_iterator &cell, + std::vector > &a) const; /** - * Computes all support points of the mapping shape functions. The inner - * 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. + * 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. + * + * Needed by the @p compute_support_points_laplace function. For + * dim=1 and dim=2 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. */ - void compute_support_points_laplace( - const typename Triangulation::cell_iterator &cell, - std::vector > &a) const; + virtual + void + add_quad_support_points(const typename Triangulation::cell_iterator &cell, + std::vector > &a) const; /** - * Needed by the @p laplace_on_quad function (for dim==2). Filled + * Needed by the @p support_point_weights_on_quad function (for dim==2). Filled * by the constructor. * - * Sizes: laplace_on_quad_vector.size()= number of inner unit_support_points - * laplace_on_quad_vector[i].size()= number of outer unit_support_points, + * Sizes: support_point_weights_on_quad.size()= number of inner unit_support_points + * support_point_weights_on_quad[i].size()= number of outer unit_support_points, * i.e. unit_support_points on the boundary of the quad * * For the definition of this vector see equation (8) of the `mapping' * report. */ - Table<2,double> laplace_on_quad_vector; + Table<2,double> support_point_weights_on_quad; /** - * Needed by the @p laplace_on_hex function (for dim==3). Filled by + * Needed by the @p support_point_weights_on_hex function (for dim==3). Filled by * the constructor. * * For the definition of this vector see equation (8) of the `mapping' * report. */ - Table<2,double> laplace_on_hex_vector; + Table<2,double> support_point_weights_on_hex; /** * Exception. @@ -490,63 +444,15 @@ protected: #ifndef DOXYGEN template <> -void MappingQ<1>::set_laplace_on_quad_vector(Table<2,double> &) const; +void MappingQ<1>::set_support_point_weights_on_quad(Table<2,double> &) const; template <> -void MappingQ<3>::set_laplace_on_hex_vector(Table<2,double> &lohvs) const; +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; -template<> -void MappingQ<3>::add_quad_support_points(const Triangulation<3>::cell_iterator &cell, - std::vector > &a) const; - -// ---- Templated functions ---- // -template -template -void -MappingQ::get_intermediate_points_on_object(const Manifold &manifold, - const TriaIterator &iter, - std::vector > &points) const -{ - const unsigned int structdim = TriaIterator::AccessorType::structure_dimension; - - // Try backward compatibility option. - if (const Boundary *boundary - = dynamic_cast *>(&manifold)) - // This is actually a boundary. Call old methods. - { - switch (structdim) - { - case 1: - { - const typename Triangulation::line_iterator line = iter; - boundary->get_intermediate_points_on_line(line, points); - return; - } - case 2: - { - const typename Triangulation::quad_iterator quad = iter; - boundary->get_intermediate_points_on_quad(quad, points); - return; - } - default: - Assert(false, ExcInternalError()); - return; - } - } - else - { - std::vector > sp(GeometryInfo::vertices_per_cell); - for (unsigned int i=0; ivertex(i); - get_intermediate_points(manifold, sp, points); - } -} - - #endif // DOXYGEN DEAL_II_NAMESPACE_CLOSE diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc index 13443f4338..9db5d010b2 100644 --- a/source/fe/mapping_q.cc +++ b/source/fe/mapping_q.cc @@ -83,13 +83,13 @@ MappingQ::MappingQ (const unsigned int degree, Assert(n_inner+n_outer==Utilities::fixed_power(degree+1), ExcInternalError()); - // build laplace_on_quad_vector + // build support_point_weights_on_quad if (degree>1) { if (dim >= 2) - set_laplace_on_quad_vector(laplace_on_quad_vector); + set_support_point_weights_on_quad(support_point_weights_on_quad); if (dim >= 3) - set_laplace_on_hex_vector(laplace_on_hex_vector); + set_support_point_weights_on_hex(support_point_weights_on_hex); } } @@ -111,13 +111,13 @@ MappingQ::MappingQ (const MappingQ &mapping) Assert(n_inner+n_outer==Utilities::fixed_power(this->polynomial_degree+1), ExcInternalError()); - // build laplace_on_quad_vector + // build support_point_weights_on_quad if (this->polynomial_degree>1) { if (dim >= 2) - set_laplace_on_quad_vector(laplace_on_quad_vector); + set_support_point_weights_on_quad(support_point_weights_on_quad); if (dim >= 3) - set_laplace_on_hex_vector(laplace_on_hex_vector); + set_support_point_weights_on_hex(support_point_weights_on_hex); } } @@ -341,7 +341,7 @@ fill_fe_subface_values (const typename Triangulation::cell_iterato template <> void -MappingQ<1>::set_laplace_on_quad_vector(Table<2,double> &) const +MappingQ<1>::set_support_point_weights_on_quad(Table<2,double> &) const { Assert(false, ExcInternalError()); } @@ -350,7 +350,7 @@ MappingQ<1>::set_laplace_on_quad_vector(Table<2,double> &) const template void -MappingQ::set_laplace_on_quad_vector(Table<2,double> &loqvs) const +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); @@ -400,7 +400,7 @@ MappingQ::set_laplace_on_quad_vector(Table<2,double> &loqvs) const else if (dim == 3) { MappingQ<2,2> mapping_2d(this->polynomial_degree); - loqvs = mapping_2d.laplace_on_quad_vector; + loqvs = mapping_2d.support_point_weights_on_quad; } } @@ -416,7 +416,7 @@ MappingQ::set_laplace_on_quad_vector(Table<2,double> &loqvs) const template <> void -MappingQ<3>::set_laplace_on_hex_vector(Table<2,double> &lohvs) const +MappingQ<3>::set_support_point_weights_on_hex(Table<2,double> &lohvs) const { Assert(this->polynomial_degree>1, ExcInternalError()); @@ -459,7 +459,7 @@ MappingQ<3>::set_laplace_on_hex_vector(Table<2,double> &lohvs) const template void -MappingQ::set_laplace_on_hex_vector(Table<2,double> &) const +MappingQ::set_support_point_weights_on_hex(Table<2,double> &) const { Assert(false, ExcInternalError()); } @@ -541,107 +541,176 @@ MappingQ::compute_laplace_vector(Table<2,double> &lvs) const -template -void -MappingQ::apply_laplace_vector(const Table<2,double> &lvs, - std::vector > &a) const -{ - // check whether the data we need is really available. if you fail here and - // if lvs==laplace_on_quad_vector in the calling function, then we didn't - // 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(this->polynomial_degree)); - - const unsigned int n_inner_apply=lvs.n_rows(); - 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, - ExcDimensionMismatch(a.size(), n_outer_apply)); - - // compute each inner point as linear combination of the outer points. the - // weights are given by the lvs entries, the outer points are the first - // (existing) elements of a - for (unsigned int unit_point=0; unit_point p; - for (unsigned int k=0; k -void -MappingQ::compute_mapping_support_points( - const typename Triangulation::cell_iterator &cell, - std::vector > &a) const -{ - // if this is a cell for which we want to compute the full mapping, then get - // them from the following function - if (use_mapping_q_on_all_cells || cell->has_boundary_lines()) - compute_support_points_laplace(cell, a); - else - // otherwise: use a Q1 mapping for which the mapping shape function - // support points are simply the vertices of the cell - { - a.resize(GeometryInfo::vertices_per_cell); - - for (unsigned int i=0; i::vertices_per_cell; ++i) - a[i] = cell->vertex(i); - } -} - - -template -void -MappingQ::compute_support_points_laplace(const typename Triangulation::cell_iterator &cell, - std::vector > &a) const +namespace { - // in any case, we need the vertices first - a.resize(GeometryInfo::vertices_per_cell); - for (unsigned int i=0; i::vertices_per_cell; ++i) - a[i] = cell->vertex(i); - - if (this->polynomial_degree>1) - switch (dim) + /** + * Ask the manifold descriptor to return intermediate points on lines or + * faces. The function needs to return one or multiple points (depending on + * the number of elements in the output vector @p points that lie inside a + * line, quad or hex). Whether it is a line, quad or hex doesn't really + * matter to this function but it can be inferred from the number of input + * points in the @p surrounding_points vector. + */ + template + void + get_intermediate_points (const Manifold &manifold, + const QGaussLobatto<1> &line_support_points, + const std::vector > &surrounding_points, + std::vector > &points) + { + Assert(surrounding_points.size() >= 2, ExcMessage("At least 2 surrounding points are required")); + const unsigned int n=points.size(); + Assert(n>0, ExcMessage("You can't ask for 0 intermediate points.")); + std::vector w(surrounding_points.size()); + + switch (surrounding_points.size()) { - case 1: - add_line_support_points(cell, a); - break; case 2: - // in 2d, add the points on the four bounding lines to the exterior - // (outer) points - add_line_support_points (cell, a); - if (dim != spacedim) - add_quad_support_points(cell, a); - else - apply_laplace_vector (laplace_on_quad_vector,a); + { + // If two points are passed, these are the two vertices, and + // we can only compute degree-1 intermediate points. + for (unsigned int i=0; i quadrature(surrounding_points, w); + points[i] = manifold.get_new_point(quadrature); + } break; + } - case 3: + case 4: { - // in 3d also add the points located on the boundary faces - add_line_support_points (cell, a); - add_quad_support_points (cell, a); - apply_laplace_vector (laplace_on_hex_vector, a); + Assert(spacedim >= 2, ExcImpossibleInDim(spacedim)); + const unsigned m= + static_cast(std::sqrt(static_cast(n))); + // is n a square number + Assert(m*m==n, ExcInternalError()); + + // If four points are passed, these are the two vertices, and + // we can only compute (degree-1)*(degree-1) intermediate + // points. + for (unsigned int i=0; i quadrature(surrounding_points, w); + points[i*m+j]=manifold.get_new_point(quadrature); + } + } break; } - default: + + case 8: Assert(false, ExcNotImplemented()); break; + default: + Assert(false, ExcInternalError()); + break; + } + } + + + + + /** + * Ask the manifold descriptor to return intermediate points on the object + * pointed to by the TriaIterator @p iter. This function tries to be + * backward compatible with respect to the differences between + * Boundary and Manifold, querying the first + * whenever the passed @p manifold can be upgraded to a + * Boundary. + */ + template + void get_intermediate_points_on_object(const Manifold &manifold, + const QGaussLobatto<1> &line_support_points, + const TriaIterator &iter, + std::vector > &points) + { + const unsigned int structdim = TriaIterator::AccessorType::structure_dimension; + + // Try backward compatibility option. + if (const Boundary *boundary + = dynamic_cast *>(&manifold)) + // This is actually a boundary. Call old methods. + { + switch (structdim) + { + case 1: + { + const typename Triangulation::line_iterator line = iter; + boundary->get_intermediate_points_on_line(line, points); + return; + } + case 2: + { + const typename Triangulation::quad_iterator quad = iter; + boundary->get_intermediate_points_on_quad(quad, points); + return; + } + default: + Assert(false, ExcInternalError()); + return; + } } + else + { + std::vector > sp(GeometryInfo::vertices_per_cell); + for (unsigned int i=0; ivertex(i); + get_intermediate_points(manifold, line_support_points, sp, points); + } + } + + + /** + * Take a support_point_weights_on_hex(quad) and apply it to the vector + * @p a to compute the inner support points as a linear combination of the + * exterior points. + * + * The vector @p a initially contains the locations of the @p n_outer + * points, the @p n_inner computed inner points are appended. + * + * See equation (7) of the `mapping' report. + */ + template + void add_weighted_interior_points(const Table<2,double> &lvs, + std::vector > &a) + { + const unsigned int n_inner_apply=lvs.n_rows(); + const unsigned int n_outer_apply=lvs.n_cols(); + Assert(a.size()==n_outer_apply, + ExcDimensionMismatch(a.size(), n_outer_apply)); + + // compute each inner point as linear combination of the outer points. the + // weights are given by the lvs entries, the outer points are the first + // (existing) elements of a + for (unsigned int unit_point=0; unit_point p; + for (unsigned int k=0; k +template void -MappingQ::add_line_support_points (const typename Triangulation::cell_iterator &cell, - std::vector > &a) const +MappingQ:: +add_line_support_points (const typename Triangulation::cell_iterator &cell, + std::vector > &a) const { // if we only need the midpoint, then ask for it. if (this->polynomial_degree==2) @@ -655,10 +724,13 @@ MappingQ::add_line_support_points (const typename Triangulation &manifold = ( ( line->manifold_id() == numbers::invalid_manifold_id ) && - ( dim < spacedim ) ? cell->get_manifold() : + ( dim < spacedim ) + ? + cell->get_manifold() + : line->get_manifold() ); a.push_back(manifold.get_new_point_on_line(line)); - }; + } } else // otherwise call the more complicated functions and ask for inner points @@ -683,7 +755,7 @@ MappingQ::add_line_support_points (const typename Triangulationget_manifold() : line->get_manifold() ); - get_intermediate_points_on_object (manifold, line, line_points); + get_intermediate_points_on_object (manifold, line_support_points, line, line_points); if (dim==3) { @@ -698,16 +770,15 @@ MappingQ::add_line_support_points (const typename Triangulation +template <> void -MappingQ<3>:: +MappingQ<3,3>:: add_quad_support_points(const Triangulation<3>::cell_iterator &cell, std::vector > &a) const { @@ -718,9 +789,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 ((this->polynomial_degree-1)*(this->polynomial_degree-1)); + std::vector > quad_points ((polynomial_degree-1)*(polynomial_degree-1)); // used if only one line of face quad is at boundary - std::vector > b(4*this->polynomial_degree); + std::vector > b(4*polynomial_degree); // Used by the new Manifold interface. This vector collects the // vertices used to compute the intermediate points. @@ -758,12 +829,12 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell, // points on it if (face->at_boundary()) { - get_intermediate_points_on_object(face->get_manifold(), face, quad_points); + get_intermediate_points_on_object(face->get_manifold(), line_support_points, face, quad_points); // in 3D, the orientation, flip and rotation of the face might not // match what we expect here, namely the standard orientation. thus // reorder points accordingly. since a Mapping uses the same shape - // function as an FEQ, we can ask a FEQ to do the reordering for us. + // function as an FE_Q, we can ask a FE_Q to do the reordering for us. for (unsigned int i=0; i::cell_iterator &cell, // boundary, then collect points separately if (lines_at_boundary>0) { - // call of function apply_laplace_vector increases size of b + // call of function add_weighted_interior_points increases size of b // about 1. There resize b for the case the mentioned function // was already called. - b.resize(4*this->polynomial_degree); + b.resize(4*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*(this->polynomial_degree-1), + Assert(b.size()==vertices_per_face+lines_per_face*(polynomial_degree-1), ExcDimensionMismatch(b.size(), - vertices_per_face+lines_per_face*(this->polynomial_degree-1))); + vertices_per_face+lines_per_face*(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 @@ -803,19 +874,20 @@ 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]= + for (unsigned int j=0; j::face_to_cell_lines( - face_no, i)*(this->polynomial_degree-1)+j]; + face_no, i)*(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*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))); + add_weighted_interior_points (support_point_weights_on_quad, b); + AssertDimension (b.size(), + 4*this->polynomial_degree + + (this->polynomial_degree-1)*(this->polynomial_degree-1)); - for (unsigned int i=0; i<(this->polynomial_degree-1)*(this->polynomial_degree-1); ++i) - a.push_back(b[4*this->polynomial_degree+i]); + for (unsigned int i=0; i<(polynomial_degree-1)*(polynomial_degree-1); ++i) + a.push_back(b[4*polynomial_degree+i]); } else { @@ -824,7 +896,7 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell, vertices.resize(4); for (unsigned int i=0; i<4; ++i) vertices[i] = face->vertex(i); - get_intermediate_points (face->get_manifold(), vertices, quad_points); + get_intermediate_points (face->get_manifold(), line_support_points, vertices, quad_points); // in 3D, the orientation, flip and rotation of the face might // not match what we expect here, namely the standard // orientation. thus reorder points accordingly. since a Mapping @@ -842,27 +914,79 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell, -template<> +template <> void MappingQ<2,3>:: add_quad_support_points(const Triangulation<2,3>::cell_iterator &cell, std::vector > &a) const { - std::vector > quad_points ((this->polynomial_degree-1)*(this->polynomial_degree-1)); - get_intermediate_points_on_object (cell->get_manifold(), cell, quad_points); + std::vector > quad_points ((polynomial_degree-1)*(polynomial_degree-1)); + get_intermediate_points_on_object (cell->get_manifold(), line_support_points, + cell, quad_points); for (unsigned int i=0; i +template void MappingQ:: add_quad_support_points(const typename Triangulation::cell_iterator &, std::vector > &) const { - Assert (dim > 2, ExcImpossibleInDim(dim)); + Assert (false, ExcInternalError()); +} + + + + + +template +void +MappingQ:: +compute_mapping_support_points(const typename Triangulation::cell_iterator &cell, + std::vector > &a) const +{ + // get the vertices first + a.resize(GeometryInfo::vertices_per_cell); + for (unsigned int i=0; i::vertices_per_cell; ++i) + a[i] = cell->vertex(i); + + if (this->polynomial_degree>1) + switch (dim) + { + case 1: + add_line_support_points(cell, a); + break; + case 2: + // in 2d, add the points on the four bounding lines to the exterior + // (outer) points + add_line_support_points(cell, a); + + // then get the support points on the quad if we are on a + // manifold, otherwise compute them from the points around it + if (dim != spacedim) + add_quad_support_points(cell, a); + else + add_weighted_interior_points (support_point_weights_on_quad, a); + break; + + case 3: + { + // in 3d also add the points located on the boundary faces + add_line_support_points (cell, a); + add_quad_support_points (cell, a); + + // then compute the interior points + add_weighted_interior_points (support_point_weights_on_hex, a); + break; + } + + default: + Assert(false, ExcNotImplemented()); + break; + } } @@ -1093,76 +1217,6 @@ MappingQ::clone () const } -template -void -MappingQ::get_intermediate_points (const Manifold &manifold, - const std::vector > &surrounding_points, - std::vector > &points) const -{ - Assert(surrounding_points.size() >= 2, ExcMessage("At least 2 surrounding points are required")); - const unsigned int n=points.size(); - Assert(n>0, ExcMessage("You can't ask for 0 intermediate points.")); - std::vector w(surrounding_points.size()); - - switch (surrounding_points.size()) - { - case 2: - { - // If two points are passed, these are the two vertices, and - // we can only compute degree-1 intermediate points. - AssertDimension(n, this->polynomial_degree-1); - for (unsigned int i=0; i quadrature(surrounding_points, w); - points[i] = manifold.get_new_point(quadrature); - } - break; - } - - case 4: - { - Assert(spacedim >= 2, ExcImpossibleInDim(spacedim)); - const unsigned m= - static_cast(std::sqrt(static_cast(n))); - // is n a square number - Assert(m*m==n, ExcInternalError()); - - // If four points are passed, these are the two vertices, and - // we can only compute (degree-1)*(degree-1) intermediate - // points. - AssertDimension(m, this->polynomial_degree-1); - - for (unsigned int i=0; i quadrature(surrounding_points, w); - points[i*m+j]=manifold.get_new_point(quadrature); - } - } - break; - } - - case 8: - Assert(false, ExcNotImplemented()); - break; - default: - Assert(false, ExcInternalError()); - break; - } -} - - // explicit instantiations #include "mapping_q.inst" -- 2.39.5