From 743feb32f262f89cb3d92d24b8dfc05eb0f122fe Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Sun, 13 Sep 2015 13:49:23 -0500 Subject: [PATCH] Move compute_mapping_support_points() and friends to MappingQGeneric. --- include/deal.II/fe/mapping_q.h | 115 +---- include/deal.II/fe/mapping_q_generic.h | 112 +++- source/fe/mapping_q.cc | 672 +----------------------- source/fe/mapping_q_generic.cc | 683 ++++++++++++++++++++++++- 4 files changed, 800 insertions(+), 782 deletions(-) diff --git a/include/deal.II/fe/mapping_q.h b/include/deal.II/fe/mapping_q.h index 72e8678c09..81f16ea5bf 100644 --- a/include/deal.II/fe/mapping_q.h +++ b/include/deal.II/fe/mapping_q.h @@ -18,11 +18,8 @@ #include -#include #include #include -#include -#include DEAL_II_NAMESPACE_OPEN @@ -266,94 +263,6 @@ protected: protected: - /** - * 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 dim::cell_iterator &cell, - std::vector > &a) const; - - - /** - * 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) 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() 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. - */ - virtual - void - add_quad_support_points(const typename Triangulation::cell_iterator &cell, - std::vector > &a) const; - - /** - * Needed by the @p support_point_weights_on_quad function (for dim==2). Filled - * by the constructor. - * - * 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> support_point_weights_on_quad; - - /** - * 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> support_point_weights_on_hex; - /** * Exception. */ @@ -377,15 +286,6 @@ protected: */ const bool use_mapping_q_on_all_cells; - /** - * An FE_Q object which is only needed in 3D, since it knows how to reorder - * shape functions/DoFs on non-standard faces. This is used to reorder - * support points in the same way. We could make this a pointer to prevent - * construction in 1D and 2D, but since memory and time requirements are not - * particularly high this seems unnecessary at the moment. - */ - const FE_Q feq; - /** * 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 @@ -403,16 +303,11 @@ protected: */ std_cxx11::unique_ptr > q1_mapping; - /* - * The default line support points. These are used when computing - * the location in real space of the support points on lines and - * quads, which are asked to the Manifold class. - * - * The number of quadrature points depends on the degree of this - * class, and it matches the number of degrees of freedom of an - * FE_Q<1>(this->degree). - */ - QGaussLobatto<1> line_support_points; + //TODO: Remove again -- all the function does is bypass the inherited function from MappingQ1 and go back to the one in the MappingQGeneric base class + virtual + void + compute_mapping_support_points (const typename Triangulation::cell_iterator &cell, + std::vector > &a) const; /** * Declare other MappingQ classes friends. diff --git a/include/deal.II/fe/mapping_q_generic.h b/include/deal.II/fe/mapping_q_generic.h index f7b4ae5afd..46fa2b4b85 100644 --- a/include/deal.II/fe/mapping_q_generic.h +++ b/include/deal.II/fe/mapping_q_generic.h @@ -20,10 +20,12 @@ #include #include #include +#include #include #include #include #include +#include #include @@ -66,6 +68,11 @@ public: */ MappingQGeneric (const unsigned int polynomial_degree); + /** + * Copy constructor. + */ + MappingQGeneric (const MappingQGeneric &mapping); + /** * Return the degree of the mapping, i.e. the value which was passed to the * constructor. @@ -453,6 +460,48 @@ protected: */ const unsigned int polynomial_degree; + /* + * The default line support points. These are used when computing + * the location in real space of the support points on lines and + * quads, which are asked to the Manifold class. + * + * The number of quadrature points depends on the degree of this + * class, and it matches the number of degrees of freedom of an + * FE_Q<1>(this->degree). + */ + QGaussLobatto<1> line_support_points; + + /** + * An FE_Q object which is only needed in 3D, since it knows how to reorder + * shape functions/DoFs on non-standard faces. This is used to reorder + * support points in the same way. + */ + const std_cxx11::unique_ptr > fe_q; + + /** + * A table of weights by which we multiply the locations of the + * support points on the perimeter of a quad to get the location of + * interior 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> support_point_weights_on_quad; + + /** + * A table of weights by which we multiply the locations of the + * support points on the perimeter of a hex to get the location of + * interior support points. + * + * For the definition of this vector see equation (8) of the `mapping' + * report. + */ + Table<2,double> support_point_weights_on_hex; + /** * An interface that derived classes have to implement and that * computes the locations of support points for the mapping. For @@ -462,11 +511,72 @@ protected: * the support points from the geometry of the current cell but * instead evaluating an externally given displacement field in * addition to the geometry of the cell. + * + * The default implementation of this function is appropriate for + * most cases. It takes the locations of support points on the + * boundary of the cell from the underlying manifold. Interior + * support points (ie. support points in quads for 2d, in hexes for + * 3d) are then 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 dim::cell_iterator &cell, - std::vector > &a) const = 0; + std::vector > &a) const; + + /** + * 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) 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() 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. + */ + virtual + void + add_quad_support_points(const typename Triangulation::cell_iterator &cell, + std::vector > &a) const; /** * Make MappingQ a friend since it needs to call the diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc index 62b02f8197..6ac2a071d5 100644 --- a/source/fe/mapping_q.cc +++ b/source/fe/mapping_q.cc @@ -22,7 +22,6 @@ #include #include #include -#include #include #include #include @@ -53,222 +52,6 @@ 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) @@ -290,17 +73,12 @@ MappingQ::MappingQ (const unsigned int degree, use_mapping_q_on_all_cells || (dim != spacedim)), - feq(degree), // 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) + q1_mapping (new MappingQ1()) { Assert(n_inner+n_outer==Utilities::fixed_power(degree+1), ExcInternalError()); - - 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); } @@ -312,17 +90,12 @@ MappingQ::MappingQ (const MappingQ &mapping) 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()), // 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) + q1_mapping (dynamic_cast*>(mapping.q1_mapping->clone())) { Assert(n_inner+n_outer==Utilities::fixed_power(this->polynomial_degree+1), ExcInternalError()); - - 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); } @@ -542,452 +315,13 @@ fill_fe_subface_values (const typename Triangulation::cell_iterato } -namespace -{ - /** - * 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 2: - { - // 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 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. - 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; - } - } - - - - - /** - * 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 -void -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) - { - for (unsigned int line_no=0; line_no::lines_per_cell; ++line_no) - { - const typename Triangulation::line_iterator line = - (dim == 1 ? - static_cast::line_iterator>(cell) : - cell->line(line_no)); - - const Manifold &manifold = - ( ( line->manifold_id() == numbers::invalid_manifold_id ) && - ( 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 - // from the boundary description - { - 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) - { - const typename Triangulation::line_iterator - line = (dim == 1 - ? - static_cast::line_iterator>(cell) - : - cell->line(line_no)); - - const Manifold &manifold = - ( ( line->manifold_id() == numbers::invalid_manifold_id ) && - ( dim < spacedim ) - ? - cell->get_manifold() : - line->get_manifold() ); - - get_intermediate_points_on_object (manifold, line_support_points, line, line_points); - - if (dim==3) - { - // in 3D, lines might be in wrong orientation. if so, reverse - // the vector - if (cell->line_orientation(line_no)) - a.insert (a.end(), line_points.begin(), line_points.end()); - else - a.insert (a.end(), line_points.rbegin(), line_points.rend()); - } - else - // in 2D, lines always have the correct orientation. simply append - // all points - a.insert (a.end(), line_points.begin(), line_points.end()); - } - } -} - - - -template <> -void -MappingQ<3,3>:: -add_quad_support_points(const Triangulation<3>::cell_iterator &cell, - std::vector > &a) const -{ - const unsigned int faces_per_cell = GeometryInfo<3>::faces_per_cell, - vertices_per_face = GeometryInfo<3>::vertices_per_face, - lines_per_face = GeometryInfo<3>::lines_per_face, - vertices_per_cell = GeometryInfo<3>::vertices_per_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 ((polynomial_degree-1)*(polynomial_degree-1)); - // used if only one line of face quad is at boundary - std::vector > b(4*polynomial_degree); - - // Used by the new Manifold interface. This vector collects the - // vertices used to compute the intermediate points. - std::vector > vertices(4); - - // loop over all faces and collect points on them - for (unsigned int face_no=0; face_no::face_iterator face = cell->face(face_no); - - // select the correct mappings for the present face - const bool face_orientation = cell->face_orientation(face_no), - face_flip = cell->face_flip (face_no), - face_rotation = cell->face_rotation (face_no); - -#ifdef DEBUG - // some sanity checks up front - for (unsigned int i=0; ivertex_index(i)==cell->vertex_index( - GeometryInfo<3>::face_to_cell_vertices(face_no, i, - face_orientation, - face_flip, - face_rotation)), - ExcInternalError()); - - // indices of the lines that bound a face are given by GeometryInfo<3>:: - // face_to_cell_lines - for (unsigned int i=0; iline(i)==cell->line(GeometryInfo<3>::face_to_cell_lines( - face_no, i, face_orientation, face_flip, face_rotation)), - ExcInternalError()); -#endif - - // if face at boundary, then ask boundary object to return intermediate - // points on it - if (face->at_boundary()) - { - 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 FE_Q, we can ask a FE_Q to do the reordering for us. - for (unsigned int i=0; iline(i)->at_boundary()) - ++lines_at_boundary; - - Assert(lines_at_boundary<=lines_per_face, ExcInternalError()); - - // if at least one of the lines bounding this quad is at the - // boundary, then collect points separately - if (lines_at_boundary>0) - { - // 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*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*(polynomial_degree-1), - ExcDimensionMismatch(b.size(), - 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 - // orientation. Doing so, the calculated points will be in - // standard orientation as well. - for (unsigned int i=0; i::face_to_cell_vertices(face_no, i)]; - - for (unsigned int i=0; i::face_to_cell_lines( - face_no, i)*(polynomial_degree-1)+j]; - - // Now b includes the support points on the quad and we can - // apply the laplace vector - 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<(polynomial_degree-1)*(polynomial_degree-1); ++i) - a.push_back(b[4*polynomial_degree+i]); - } - else - { - // face is entirely in the interior. get intermediate - // points from the relevant manifold object. - vertices.resize(4); - for (unsigned int i=0; i<4; ++i) - vertices[i] = face->vertex(i); - 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 - // uses the same shape function as an FEQ, we can ask a FEQ to - // do the reordering for us. - for (unsigned int i=0; i -void -MappingQ<2,3>:: -add_quad_support_points(const Triangulation<2,3>::cell_iterator &cell, - std::vector > &a) const -{ - 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 -void -MappingQ:: -add_quad_support_points(const typename Triangulation::cell_iterator &, - std::vector > &) const -{ - 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; - } + MappingQGeneric::compute_mapping_support_points (cell, a); } diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc index 155e6a23a9..adaac27fe7 100644 --- a/source/fe/mapping_q_generic.cc +++ b/source/fe/mapping_q_generic.cc @@ -25,6 +25,7 @@ #include #include #include +#include #include #include #include @@ -767,11 +768,243 @@ compute_shape_function_values (const std::vector > &unit_points) } +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 MappingQGeneric::MappingQGeneric (const unsigned int p) : - polynomial_degree(p) + polynomial_degree(p), + line_support_points(this->polynomial_degree+1), + fe_q(dim == 3 ? new FE_Q(this->polynomial_degree) : 0), + 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)) +{} + + + +template +MappingQGeneric::MappingQGeneric (const MappingQGeneric &mapping) + : + polynomial_degree(mapping.polynomial_degree), + line_support_points(mapping.line_support_points), + fe_q(dim == 3 ? new FE_Q(*mapping.fe_q) : 0), + support_point_weights_on_quad (mapping.support_point_weights_on_quad), + support_point_weights_on_hex (mapping.support_point_weights_on_hex) {} @@ -2281,7 +2514,6 @@ transform (const VectorSlice > > input, const typename Mapping::InternalDataBase &mapping_data, VectorSlice > > output) const { - switch (mapping_type) { case mapping_piola_hessian: @@ -2292,7 +2524,454 @@ transform (const VectorSlice > > input, default: Assert(false, ExcNotImplemented()); } +} + + + +namespace +{ + /** + * 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 2: + { + // 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 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. + 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; + } + } + + + + + /** + * 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 +void +MappingQGeneric:: +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) + { + for (unsigned int line_no=0; line_no::lines_per_cell; ++line_no) + { + const typename Triangulation::line_iterator line = + (dim == 1 ? + static_cast::line_iterator>(cell) : + cell->line(line_no)); + + const Manifold &manifold = + ( ( line->manifold_id() == numbers::invalid_manifold_id ) && + ( 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 + // from the boundary description + { + 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) + { + const typename Triangulation::line_iterator + line = (dim == 1 + ? + static_cast::line_iterator>(cell) + : + cell->line(line_no)); + + const Manifold &manifold = + ( ( line->manifold_id() == numbers::invalid_manifold_id ) && + ( dim < spacedim ) + ? + cell->get_manifold() : + line->get_manifold() ); + + get_intermediate_points_on_object (manifold, line_support_points, line, line_points); + + if (dim==3) + { + // in 3D, lines might be in wrong orientation. if so, reverse + // the vector + if (cell->line_orientation(line_no)) + a.insert (a.end(), line_points.begin(), line_points.end()); + else + a.insert (a.end(), line_points.rbegin(), line_points.rend()); + } + else + // in 2D, lines always have the correct orientation. simply append + // all points + a.insert (a.end(), line_points.begin(), line_points.end()); + } + } +} + + + +template <> +void +MappingQGeneric<3,3>:: +add_quad_support_points(const Triangulation<3>::cell_iterator &cell, + std::vector > &a) const +{ + const unsigned int faces_per_cell = GeometryInfo<3>::faces_per_cell, + vertices_per_face = GeometryInfo<3>::vertices_per_face, + lines_per_face = GeometryInfo<3>::lines_per_face, + vertices_per_cell = GeometryInfo<3>::vertices_per_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 ((polynomial_degree-1)*(polynomial_degree-1)); + // used if only one line of face quad is at boundary + std::vector > b(4*polynomial_degree); + + // Used by the new Manifold interface. This vector collects the + // vertices used to compute the intermediate points. + std::vector > vertices(4); + + // loop over all faces and collect points on them + for (unsigned int face_no=0; face_no::face_iterator face = cell->face(face_no); + + // select the correct mappings for the present face + const bool face_orientation = cell->face_orientation(face_no), + face_flip = cell->face_flip (face_no), + face_rotation = cell->face_rotation (face_no); + +#ifdef DEBUG + // some sanity checks up front + for (unsigned int i=0; ivertex_index(i)==cell->vertex_index( + GeometryInfo<3>::face_to_cell_vertices(face_no, i, + face_orientation, + face_flip, + face_rotation)), + ExcInternalError()); + + // indices of the lines that bound a face are given by GeometryInfo<3>:: + // face_to_cell_lines + for (unsigned int i=0; iline(i)==cell->line(GeometryInfo<3>::face_to_cell_lines( + face_no, i, face_orientation, face_flip, face_rotation)), + ExcInternalError()); +#endif + + // if face at boundary, then ask boundary object to return intermediate + // points on it + if (face->at_boundary()) + { + 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 FE_Q, we can ask a FE_Q to do the reordering for us. + for (unsigned int i=0; iadjust_quad_dof_index_for_face_orientation(i, + face_orientation, + face_flip, + face_rotation)]); + } + else + { + // face is not at boundary, but maybe some of its lines are. count + // them + unsigned int lines_at_boundary=0; + for (unsigned int i=0; iline(i)->at_boundary()) + ++lines_at_boundary; + + Assert(lines_at_boundary<=lines_per_face, ExcInternalError()); + + // if at least one of the lines bounding this quad is at the + // boundary, then collect points separately + if (lines_at_boundary>0) + { + // 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*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*(polynomial_degree-1), + ExcDimensionMismatch(b.size(), + 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 + // orientation. Doing so, the calculated points will be in + // standard orientation as well. + for (unsigned int i=0; i::face_to_cell_vertices(face_no, i)]; + + for (unsigned int i=0; i::face_to_cell_lines( + face_no, i)*(polynomial_degree-1)+j]; + + // Now b includes the support points on the quad and we can + // apply the laplace vector + 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<(polynomial_degree-1)*(polynomial_degree-1); ++i) + a.push_back(b[4*polynomial_degree+i]); + } + else + { + // face is entirely in the interior. get intermediate + // points from the relevant manifold object. + vertices.resize(4); + for (unsigned int i=0; i<4; ++i) + vertices[i] = face->vertex(i); + 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 + // uses the same shape function as an FE_Q, we can ask a FE_Q to + // do the reordering for us. + for (unsigned int i=0; iadjust_quad_dof_index_for_face_orientation(i, + face_orientation, + face_flip, + face_rotation)]); + } + } + } +} + + + +template <> +void +MappingQGeneric<2,3>:: +add_quad_support_points(const Triangulation<2,3>::cell_iterator &cell, + std::vector > &a) const +{ + 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 +void +MappingQGeneric:: +add_quad_support_points(const typename Triangulation::cell_iterator &, + std::vector > &) const +{ + Assert (false, ExcInternalError()); +} + + + +template +void +MappingQGeneric:: +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; + } } -- 2.39.5