From 92a06f6229a1c242171d6a24f3eb2c7cdd2bbf6f Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 17 Mar 2017 15:21:31 +0100 Subject: [PATCH] Unify interpolation from surrounding points. Query multiple points from mapping. --- include/deal.II/fe/mapping_q.h | 64 ++-- include/deal.II/fe/mapping_q_generic.h | 143 ++++--- source/fe/mapping_q.cc | 4 +- source/fe/mapping_q_generic.cc | 495 +++++++++++++------------ 4 files changed, 369 insertions(+), 337 deletions(-) diff --git a/include/deal.II/fe/mapping_q.h b/include/deal.II/fe/mapping_q.h index 2694e8521d..dd6946c3b6 100644 --- a/include/deal.II/fe/mapping_q.h +++ b/include/deal.II/fe/mapping_q.h @@ -72,13 +72,14 @@ DEAL_II_NAMESPACE_OPEN *

Behavior along curved boundaries and with different manifolds

* * As described above, one often only knows a manifold description of a - * surface but not the interior of the computational domain. Thus, a + * surface but not the interior of the computational domain. In such a case, a * StraightBoundary object will be assigned to the interior entities that - * describes a usual planar coordinate system where the additional points of - * this class are placed exactly according to a bi-/trilinear mapping. When - * combined with a non-flat manifold on the boundary, for example a circle, - * the two manifold descriptions are in general incompatible. For example, a - * StraightBoundary would put an interior point located at some small distance + * describes a usual planar coordinate system where the additional points for + * the higher order mapping are placed exactly according to a bi-/trilinear + * mapping. When combined with a non-flat manifold on the boundary, for + * example a circle, the two manifold descriptions are in general + * incompatible. For example, a StraightBoundary defined solely through the + * cell's vertices would put an interior point located at some small distance * epsilon away from the boundary along a flat line and thus in general * outside the concave part of a circle. If the polynomial degree of * MappingQGeneric is sufficiently high, the transformation from the reference @@ -86,28 +87,31 @@ DEAL_II_NAMESPACE_OPEN * boundary. * * In order to avoid this situation, this class applies a smoothing on cells - * adjacent to the boundary by using so-called Laplace smoothing. In the - * algorithm computing new points, all the entities of the cells are passed - * through hierarchically, starting from the lines to the quads and finally - * hexes. The elements higher up in the hierarchy that sit on a - * StraightBoundary will then get their points interpolated from all the - * surrounding points and not just the corner points. If only a line is - * associated to a curved boundary but the adjacent quad is on a flat - * manifold, the points inside the quad will be computed according to the - * deformed line and thus always result in a well-defined transformation. If - * this smoothing is undesired, the optional argument @p - * smooth_support_points_on_flat_manifold can be used to disable the - * smoothing, placing the additional points strictly according to the - * manifold. + * adjacent to the boundary by using so-called Laplace smoothing by + * default. In the algorithm that computing additional points, the + * compute_mapping_support_points() method, all the entities of the cells are + * passed through hierarchically, starting from the lines to the quads and + * finally hexes. The elements higher up in the hierarchy that sit on the + * boundary will then get their points interpolated from all the surrounding + * points and not just the corner points. If only a line is assigned a curved + * boundary but the adjacent quad is on a flat manifold, the points inside the + * quad will be computed according to the deformed line and thus always result + * in a well-defined transformation. This smoothing can be disabled by setting + * the optional argument @p smooth_support_points to false, placing the + * additional points strictly according to the manifold. This is usually the + * most efficient choice in case different manifolds are present that are + * compatible with each other. * - * While this smoothing approach works well for low and medium convergence - * orders up to approximately three to four, this mechanism has an inherent - * shortcoming because it switches from a curved manifold to a flat manifold - * within one layer of elements. This will cause the Jacobian transformation - * to have jumps between the first and second element layer that can reduce - * the order of convergence. For example, the convergence rates for solving - * the Laplacian on a circle where only the boundary is deformed and the above - * mesh smoothing algorithm is applied will typically not exceed 3.5, even for + * While the smoothing approach works well for filling holes or avoiding + * inversions with low and medium convergence orders up to approximately three + * to four, there is nonetheless an inherent shortcoming because of a + * discontinuous mapping that switches from a curved manifold to a flat + * manifold within one layer of elements. This will cause the Jacobian + * transformation to have jumps between the first and second element layer + * that can reduce the order of convergence. For example, the convergence + * rates for solving the Laplacian on a circle where only the boundary is + * deformed and the above mesh smoothing algorithm is applied will typically + * not exceed 3.5 (or 3 in the elements adjacent to the boundary), even for * fourth or fifth degree polynomials. In such a case, the curved manifold * needs to be switched to a flat manifold in a smooth way that does not * depend on the mesh size and eventuell covers a whole layer of cells. @@ -136,13 +140,13 @@ public: * equal to @p spacedim, i.e., if we are considering meshes on surfaces * embedded into higher dimensional spaces. * - * The optional parameter @p smooth_support_points_on_flat_manifold controls - * whether smoothing adjancent to StraightBoundary objects according to the + * The optional parameter @p smooth_support_points controls whether + * smoothing on objects where different manifolds meet according to the * general class description should be enabled (default) or not. */ MappingQ (const unsigned int polynomial_degree, const bool use_mapping_q_on_all_cells = false, - const bool smooth_support_points_on_flat_manifold = true); + const bool smooth_support_points = true); /** * Copy constructor. diff --git a/include/deal.II/fe/mapping_q_generic.h b/include/deal.II/fe/mapping_q_generic.h index 1832e37786..4c604c5c45 100644 --- a/include/deal.II/fe/mapping_q_generic.h +++ b/include/deal.II/fe/mapping_q_generic.h @@ -84,13 +84,14 @@ template class MappingQ; *

Behavior along curved boundaries and with different manifolds

* * As described above, one often only knows a manifold description of a - * surface but not the interior of the computational domain. Thus, a + * surface but not the interior of the computational domain. In such a case, a * StraightBoundary object will be assigned to the interior entities that - * describes a usual planar coordinate system where the additional points of - * this class are placed exactly according to a bi-/trilinear mapping. When - * combined with a non-flat manifold on the boundary, for example a circle, - * the two manifold descriptions are in general incompatible. For example, a - * StraightBoundary would put an interior point located at some small distance + * describes a usual planar coordinate system where the additional points for + * the higher order mapping are placed exactly according to a bi-/trilinear + * mapping. When combined with a non-flat manifold on the boundary, for + * example a circle, the two manifold descriptions are in general + * incompatible. For example, a StraightBoundary defined solely through the + * cell's vertices would put an interior point located at some small distance * epsilon away from the boundary along a flat line and thus in general * outside the concave part of a circle. If the polynomial degree of * MappingQGeneric is sufficiently high, the transformation from the reference @@ -98,33 +99,36 @@ template class MappingQ; * boundary. * * In order to avoid this situation, this class applies a smoothing on cells - * adjacent to the boundary by using so-called Laplace smoothing. In the - * algorithm computing new points, all the entities of the cells are passed - * through hierarchically, starting from the lines to the quads and finally - * hexes. The elements higher up in the hierarchy that sit on a - * StraightBoundary will then get their points interpolated from all the - * surrounding points and not just the corner points. If only a line is - * associated to a curved boundary but the adjacent quad is on a flat - * manifold, the points inside the quad will be computed according to the - * deformed line and thus always result in a well-defined transformation. If - * this smoothing is undesired, the optional argument @p - * smooth_support_points_on_flat_manifold can be used to disable the - * smoothing, placing the additional points strictly according to the - * manifold. + * adjacent to the boundary by using so-called Laplace smoothing by + * default. In the algorithm that computing additional points, the + * compute_mapping_support_points() method, all the entities of the cells are + * passed through hierarchically, starting from the lines to the quads and + * finally hexes. The elements higher up in the hierarchy that sit on the + * boundary will then get their points interpolated from all the surrounding + * points and not just the corner points. If only a line is assigned a curved + * boundary but the adjacent quad is on a flat manifold, the points inside the + * quad will be computed according to the deformed line and thus always result + * in a well-defined transformation. This smoothing can be disabled by setting + * the optional argument @p smooth_support_points to false, placing the + * additional points strictly according to the manifold. This is usually the + * most efficient choice in case different manifolds are present that are + * compatible with each other. * - * While this smoothing approach works well for low and medium convergence - * orders up to approximately three to four, this mechanism has an inherent - * shortcoming because it switches from a curved manifold to a flat manifold - * within one layer of elements. This will cause the Jacobian transformation - * to have jumps between the first and second element layer that can reduce - * the order of convergence. For example, the convergence rates for solving - * the Laplacian on a circle where only the boundary is deformed and the above - * mesh smoothing algorithm is applied will typically not exceed 3.5, even for + * While the smoothing approach works well for filling holes or avoiding + * inversions with low and medium convergence orders up to approximately three + * to four, there is nonetheless an inherent shortcoming because of a + * discontinuous mapping that switches from a curved manifold to a flat + * manifold within one layer of elements. This will cause the Jacobian + * transformation to have jumps between the first and second element layer + * that can reduce the order of convergence. For example, the convergence + * rates for solving the Laplacian on a circle where only the boundary is + * deformed and the above mesh smoothing algorithm is applied will typically + * not exceed 3.5 (or 3 in the elements adjacent to the boundary), even for * fourth or fifth degree polynomials. In such a case, the curved manifold * needs to be switched to a flat manifold in a smooth way that does not * depend on the mesh size and eventuell covers a whole layer of cells. * - * @author Wolfgang Bangerth, 2015 + * @author Wolfgang Bangerth, 2015, Martin Kronbichler, 2017 */ template class MappingQGeneric : public Mapping @@ -135,12 +139,12 @@ public: * polynomials that are used to map cells from the reference to the real * cell. * - * The optional parameter @p smooth_support_points_on_flat_manifold controls - * whether smoothing adjancent to StraightBoundary objects according to the - * general class description should be enabled (default) or not. + * The optional parameter @p smooth_support_points controls whether + * smoothing according to the general class description should be enabled + * (default) or not. */ MappingQGeneric (const unsigned int polynomial_degree, - const bool smooth_support_points_on_flat_manifold = true); + const bool smooth_support_points = true); /** * Copy constructor. @@ -437,16 +441,6 @@ public: */ const unsigned int n_shape_functions; - /* - * The default line support points. Is used in when the shape function - * values are computed. - * - * 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; - /** * Tensors of covariant transformation at each of the quadrature points. * The matrix stored is the Jacobian * G^{-1}, where G = Jacobian^{t} * @@ -556,7 +550,7 @@ protected: * Stores whether we want to smooth the placement of interior points on flat * manifolds. */ - const bool smooth_support_points_on_flat_manifold; + const bool smooth_support_points; /* * The default line support points. These are used when computing @@ -577,29 +571,62 @@ protected: 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. + * A vector of tables of weights by which we multiply the locations of the + * support points on the perimeter of an object (line, quad, hex) 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 + * Access into this table is by @p [structdim-1], i.e., use 0 to access the + * support point weights on a line (i.e., the interior points of the + * GaussLobatto quadrature), use 1 to access the support point weights from + * to perimeter to the interior of a quad, and use 2 to access the support + * point weights from the perimeter to the interior of a hex. * - * For the definition of this vector see equation (8) of the `mapping' + * The table itself contains as many columns as there are surrounding points + * to a particular object (2 for a vertex, 4 + 4*(degree-1) for + * a quad, 8 + 12*(degree-1) + 6*(degree-1)*(degree-1) for a + * hex) and as many rows as there are strictly interior points. + * + * For the definition of this table see equation (8) of the `mapping' * report. */ - Table<2,double> support_point_weights_on_quad; + std::vector > support_point_weights_perimeter_to_interior; /** - * 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. + * A vector of tables of weights by which we multiply the locations of the + * vertex points of an object (line, quad, hex) to get the location of + * interior support points. * - * For the definition of this vector see equation (8) of the `mapping' - * report. + * As opposed to @p support_point_weights_perimeter_to_interior, this table + * takes only the vertex points into account and not intermediate points + * inside the objects. Thus, the content of this table simply corresponds to + * the evaluation of the linear shape functions in the dimension of the + * object. + * + * Access into the vector of tables is by @p [structdim-1], i.e., use 0 to + * access the support point weights on a line (i.e., the interior points of + * the GaussLobatto quadrature), use 1 to access the support point weights + * from to perimeter to the interior of a quad, and use 2 to access the + * support point weights from the perimeter to the interior of a hex. + * + * The table itself contains as many columns as there are vertices on a + * particular object (2 for a vertex, 4 for a quad, 8 for a hex) and as many + * rows as there are strictly interior points. + */ + std::vector > support_point_weights_interior; + + /** + * A table of weights by which we multiply the locations of the vertex + * points of the cell to get the location of all additional support points, + * both on lines, quads, and hexes (as appropriate). This data structure is + * used when we fill all support points at once, which is the case if the + * same manifold is attached to all sub-entities of a cell. This way, we can + * avoid some of the overhead in transforming data for mappings. + * + * The table has as many rows as there are vertices to the cell (2 in 1D, 4 + * in 2D, 8 in 3D), and as many rows as there are additional support points + * in the mapping, i.e., (degree+1)^dim - 2^dim. */ - Table<2,double> support_point_weights_on_hex; + Table<2,double> support_point_weights_cell; /** * Return the locations of support points for the mapping. For example, for diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc index 87ed327bb8..4b0b849ac7 100644 --- a/source/fe/mapping_q.cc +++ b/source/fe/mapping_q.cc @@ -55,7 +55,7 @@ MappingQ::InternalData::memory_consumption () const template MappingQ::MappingQ (const unsigned int degree, const bool use_mapping_q_on_all_cells, - const bool smooth_support_points_on_flat_manifold) + const bool smooth_support_points) : polynomial_degree (degree), @@ -77,7 +77,7 @@ MappingQ::MappingQ (const unsigned int degree, // created via the shared_ptr objects qp_mapping (this->polynomial_degree>1 ? - std_cxx11::shared_ptr >(new MappingQGeneric(degree, smooth_support_points_on_flat_manifold)) + std_cxx11::shared_ptr >(new MappingQGeneric(degree, smooth_support_points)) : q1_mapping) {} diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc index 27a1fbb99f..4392dd8551 100644 --- a/source/fe/mapping_q_generic.cc +++ b/source/fe/mapping_q_generic.cc @@ -353,8 +353,9 @@ namespace internal // Construct the tensor product polynomials used as shape functions for the // Qp mapping of cells at the boundary. + QGaussLobatto<1> line_support_points(data.polynomial_degree+1); const TensorProductPolynomials - tensor_pols (Polynomials::generate_complete_Lagrange_basis(data.line_support_points.get_points())); + tensor_pols (Polynomials::generate_complete_Lagrange_basis(line_support_points.get_points())); Assert (n_shape_functions==tensor_pols.n(), ExcInternalError()); @@ -752,8 +753,7 @@ template MappingQGeneric::InternalData::InternalData (const unsigned int polynomial_degree) : polynomial_degree (polynomial_degree), - n_shape_functions (Utilities::fixed_power(polynomial_degree+1)), - line_support_points (polynomial_degree + 1) + n_shape_functions (Utilities::fixed_power(polynomial_degree+1)) {} @@ -1086,16 +1086,11 @@ namespace * 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) @@ -1150,22 +1145,17 @@ namespace * 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_inner = Utilities::fixed_power<3>(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 @@ -1188,7 +1178,7 @@ namespace else { // not precomputed, then do so now - lohvs = compute_laplace_vector(polynomial_degree); + lohvs = compute_laplace_vector<3>(polynomial_degree); } // the sum of weights of the points at the outer rim should be one. check @@ -1200,6 +1190,115 @@ namespace return lohvs; } + + /** + * This function collects the output of + * compute_support_point_weights_on_{quad,hex} in a single data structure. + */ + std::vector > + compute_support_point_weights_perimeter_to_interior(const unsigned int polynomial_degree, + const unsigned int dim) + { + Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim)); + std::vector > output(dim); + if (polynomial_degree <= 1) + return output; + + // fill the 1D interior weights + QGaussLobatto<1> quadrature(polynomial_degree+1); + output[0].reinit(polynomial_degree-1, GeometryInfo<1>::vertices_per_cell); + for (unsigned int q=0; q::vertices_per_cell; ++i) + output[0](q,i) = GeometryInfo<1>::d_linear_shape_function(quadrature.point(q+1), + i); + + if (dim > 1) + output[1] = compute_support_point_weights_on_quad(polynomial_degree); + + if (dim > 2) + output[2] = compute_support_point_weights_on_hex(polynomial_degree); + + return output; + } + + /** + * Collects all interior points for the various dimensions. + */ + std::vector > + compute_support_point_weights_interior(const unsigned int polynomial_degree, + const unsigned int dim) + { + Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim)); + std::vector > output(dim); + if (polynomial_degree <= 1) + return output; + + // fill the 1D interior weights + QGaussLobatto<1> quadrature(polynomial_degree+1); + output[0].reinit(polynomial_degree-1, GeometryInfo<1>::vertices_per_cell); + for (unsigned int q=0; q::vertices_per_cell; ++i) + output[0](q,i) = GeometryInfo<1>::d_linear_shape_function(quadrature.point(q+1), + i); + + // fill the 2D interior points + if (dim > 1) + { + output[1].reinit(Utilities::fixed_power<2>(polynomial_degree-1), + GeometryInfo<2>::vertices_per_cell); + for (unsigned int q=0, q2=0; q2 point(quadrature.point(q1+1)[0], + quadrature.point(q2+1)[0]); + for (unsigned int i=0; i::vertices_per_cell; ++i) + output[1](q,i) = GeometryInfo<2>::d_linear_shape_function(point, i); + } + } + + // fill the 3D interior points + if (dim > 2) + { + output[2].reinit(Utilities::fixed_power<3>(polynomial_degree-1), + GeometryInfo<3>::vertices_per_cell); + for (unsigned int q=0, q3=0; q3 point(quadrature.point(q1+1)[0], + quadrature.point(q2+1)[0], + quadrature.point(q3+1)[0]); + for (unsigned int i=0; i::vertices_per_cell; ++i) + output[2](q,i) = GeometryInfo<3>::d_linear_shape_function(point, i); + } + } + + return output; + } + + /** + * Collects all interior points for the various dimensions. + */ + template + Table<2,double> + compute_support_point_weights_cell(const unsigned int polynomial_degree) + { + Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim)); + if (polynomial_degree <= 1) + return Table<2,double>(); + + FE_Q fe(polynomial_degree); + std::vector > unit_support_points = fe.get_unit_support_points(); + + Table<2,double> output(fe.dofs_per_cell - GeometryInfo::vertices_per_cell, + GeometryInfo::vertices_per_cell); + for (unsigned int q=0; q::vertices_per_cell; ++i) + output(q,i) = GeometryInfo::d_linear_shape_function(unit_support_points[q+GeometryInfo::vertices_per_cell], + i); + + return output; + } } @@ -1207,14 +1306,15 @@ namespace template MappingQGeneric::MappingQGeneric (const unsigned int p, - const bool smooth_support_points_on_flat_manifold) + const bool smooth_support_points) : polynomial_degree(p), - smooth_support_points_on_flat_manifold(smooth_support_points_on_flat_manifold), + smooth_support_points(smooth_support_points), 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)) + support_point_weights_perimeter_to_interior (compute_support_point_weights_perimeter_to_interior(this->polynomial_degree, dim)), + support_point_weights_interior (compute_support_point_weights_interior(this->polynomial_degree, dim)), + support_point_weights_cell (compute_support_point_weights_cell(this->polynomial_degree)) { Assert (p >= 1, ExcMessage ("It only makes sense to create polynomial mappings " "with a polynomial degree greater or equal to one.")); @@ -1226,11 +1326,12 @@ template MappingQGeneric::MappingQGeneric (const MappingQGeneric &mapping) : polynomial_degree(mapping.polynomial_degree), - smooth_support_points_on_flat_manifold(mapping.smooth_support_points_on_flat_manifold), + smooth_support_points(mapping.smooth_support_points), 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) + support_point_weights_perimeter_to_interior (mapping.support_point_weights_perimeter_to_interior), + support_point_weights_interior (mapping.support_point_weights_interior), + support_point_weights_cell (mapping.support_point_weights_cell) {} @@ -3513,108 +3614,6 @@ transform (const ArrayView > &input, 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, so - // we can only compute degree-1 intermediate points. - for (unsigned int i=0; i= 2, ExcImpossibleInDim(spacedim)); - const unsigned m= - static_cast(std::sqrt(static_cast(n))); - // is n a square number - Assert(m*m==n, ExcInternalError()); - - // If m points are passed, two of them are the outer vertices, and we - // can only compute (degree-1)*(degree-1) intermediate points. - Point<2> p; - for (unsigned int i=0, c=0; i::d_linear_shape_function(p, l); - - points[c]=manifold.get_new_point(surrounding_points, w); - } - } - break; - } - - case 8: - { - Assert(spacedim >= 3, ExcImpossibleInDim(spacedim)); - unsigned int m=1; - for ( ; m < n; ++m) - if (m*m*m == n) - break; - // is n a cube number - Assert(m*m*m==n, ExcInternalError()); - - Point<3> p; - for (unsigned int k=0, c=0; k::d_linear_shape_function(p, l); - - points[c]=manifold.get_new_point(surrounding_points, w); - } - } - } - 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 @@ -3624,82 +3623,70 @@ namespace * Boundary. */ template - void get_intermediate_points_on_object(const Manifold &manifold, - const QGaussLobatto<1> &line_support_points, - const TriaIterator &iter, - std::vector > &points) + void add_intermediate_points(const Manifold &manifold, + const Table<2,double> &weight_table, + 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)) + const Boundary *boundary + = dynamic_cast *>(&manifold); + if (structdim < spacedim && boundary != NULL) // This is actually a boundary. Call old methods. { + std::vector > new_points(weight_table.size(0)); switch (structdim) { case 1: { const typename Triangulation::line_iterator line = iter; - boundary->get_intermediate_points_on_line(line, points); - return; + boundary->get_intermediate_points_on_line(line, new_points); + break; } case 2: { const typename Triangulation::quad_iterator quad = iter; - boundary->get_intermediate_points_on_quad(quad, points); - return; + boundary->get_intermediate_points_on_quad(quad, new_points); + break; } default: Assert(false, ExcInternalError()); - return; + break; } + for (unsigned int i=0; i > sp(GeometryInfo::vertices_per_cell); for (unsigned int i=0; ivertex(i); - get_intermediate_points(manifold, line_support_points, sp, points); + manifold.add_new_points(sp, weight_table, 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) + + template + bool check_identical_manifolds_of_quads(const Iterator &) { - 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::cell_iterator &cell) + { + for (unsigned int f=0; f::faces_per_cell; ++f) + if (&cell->face(f)->get_manifold() != &cell->get_manifold()) + return false; + return true; } } + template void MappingQGeneric:: @@ -3749,7 +3736,9 @@ add_line_support_points (const typename Triangulation::cell_iterat cell->get_manifold() : line->get_manifold() ); - get_intermediate_points_on_object (manifold, line_support_points, line, line_points); + line_points.resize(0); + add_intermediate_points (manifold, support_point_weights_interior[0], + line, line_points); if (dim==3) { @@ -3774,12 +3763,12 @@ template <> void MappingQGeneric<3,3>:: add_quad_support_points(const Triangulation<3,3>::cell_iterator &cell, - std::vector > &a) const + std::vector > &a) const { const unsigned int faces_per_cell = GeometryInfo<3>::faces_per_cell; // used if face quad at boundary or entirely in the interior of the domain - std::vector > quad_points ((polynomial_degree-1)*(polynomial_degree-1)); + std::vector > quad_points; std::vector > tmp_points; // loop over all faces and collect points on them @@ -3816,16 +3805,20 @@ add_quad_support_points(const Triangulation<3,3>::cell_iterator &cell, // On a quad, we have to check whether the manifold should determine the // point distribution or rather a weighted sum should be created. This // is the same logic as in the compute_mapping_support_points function - // below: if we have a StraightBoundary but no "real" boundary - // description doing something else, we should go for the smoothing of - // the support points through add_weighted_interior_points - if (smooth_support_points_on_flat_manifold == false || - std::string(typeid(face->get_manifold()).name()).find("StraightBoundary") == - std::string::npos) + // below + bool all_manifold_ids_are_equal = true; + for (unsigned int l=0; l::lines_per_face; ++l) + if (&face->line(l)->get_manifold() != &face->get_manifold()) + all_manifold_ids_are_equal = false; + if (smooth_support_points == false || all_manifold_ids_are_equal || + (dynamic_cast *>(&face->get_manifold()) && + std::string(typeid(face->get_manifold()).name()).find("StraightBoundary") == + std::string::npos)) { // ask the boundary/manifold object to return intermediate points on it - get_intermediate_points_on_object(face->get_manifold(), line_support_points, - face, quad_points); + quad_points.resize(0); + add_intermediate_points(face->get_manifold(), support_point_weights_interior[1], + 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 @@ -3852,9 +3845,9 @@ add_quad_support_points(const Triangulation<3,3>::cell_iterator &cell, a[GeometryInfo<3>::vertices_per_cell + (polynomial_degree-1)* GeometryInfo<3>::face_to_cell_lines(face_no,line) + i]; - add_weighted_interior_points (support_point_weights_on_quad, tmp_points); - a.insert(a.end(), tmp_points.begin()+4+4*(polynomial_degree-1), - tmp_points.end()); + face->get_manifold().add_new_points (tmp_points, + support_point_weights_perimeter_to_interior[1], + a); } } } @@ -3865,13 +3858,10 @@ template <> void MappingQGeneric<2,3>:: add_quad_support_points(const Triangulation<2,3>::cell_iterator &cell, - std::vector > &a) const + 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; iget_manifold(), support_point_weights_interior[1], + cell, a); } @@ -3893,72 +3883,83 @@ MappingQGeneric:: compute_mapping_support_points(const typename Triangulation::cell_iterator &cell) const { // get the vertices first - std::vector > a(GeometryInfo::vertices_per_cell); + std::vector > a; + a.reserve(Utilities::fixed_power(polynomial_degree+1)); for (unsigned int i=0; i::vertices_per_cell; ++i) - a[i] = cell->vertex(i); + a.push_back(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 have a - // StraightBoundary but no "real" boundary description doing something - // else, we should go for the smoothing of the support points through - // add_weighted_interior_points. we need to go through typeid because - // curved boundaries are derived from StraightBoundary and thus a - // plain dynamic_cast will not suffice. - if (dim != spacedim) - add_quad_support_points(cell, a); - else if (smooth_support_points_on_flat_manifold == false || - std::string(typeid(cell->get_manifold()).name()).find("StraightBoundary") == - std::string::npos) - { - std::vector > quad_points (Utilities::fixed_power(polynomial_degree-1)); - get_intermediate_points_on_object(cell->get_manifold(), line_support_points, cell, quad_points); - for (unsigned int i=0; ipolynomial_degree > 1) + { + // check if all entities have the same manifold id which is when we can + // simply ask the manifold for all points + Assert(dim<=3, ExcImpossibleInDim(dim)); + bool all_manifold_ids_are_equal = (dim == spacedim); + if (dim > 1) + for (unsigned int l=0; l::lines_per_cell; ++l) + if (&cell->line(l)->get_manifold() != &cell->get_manifold()) + all_manifold_ids_are_equal = false; + if (dim == 3) + if (check_identical_manifolds_of_quads(cell) == false) + all_manifold_ids_are_equal = false; + if (all_manifold_ids_are_equal) + { + std::vector > vertices(a); + cell->get_manifold().add_new_points(vertices, support_point_weights_cell, a); + return a; + } - 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. if we have a StraightBoundary but - // no "real" boundary description doing something else, we should go - // for the smoothing of the support points through - // add_weighted_interior_points. we need to go through typeid because - // curved boundaries are derived from StraightBoundary and thus a - // plain dynamic_cast will not suffice. - if (smooth_support_points_on_flat_manifold == false || - std::string(typeid(cell->get_manifold()).name()).find("StraightBoundary") == - std::string::npos) - { - std::vector > hex_points (Utilities::fixed_power(polynomial_degree-1)); - get_intermediate_points_on_object(cell->get_manifold(), line_support_points, cell, hex_points); - for (unsigned int i=0; iget_manifold(), support_point_weights_interior[1], + cell, a); + else + { + std::vector > tmp_points(a); + cell->get_manifold().add_new_points(tmp_points, + support_point_weights_perimeter_to_interior[1], + 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. decide whether to choose + // the smoothed variant that include all points and goes directly + // through the manifold or only the vertex points. + if (smooth_support_points == false) + add_intermediate_points(cell->get_manifold(), support_point_weights_interior[2], + cell, a); + else + { + std::vector > tmp_points(a); + cell->get_manifold().add_new_points(tmp_points, + support_point_weights_perimeter_to_interior[2], + a); + } + break; - default: - Assert(false, ExcNotImplemented()); - break; - } + default: + Assert(false, ExcNotImplemented()); + break; + } + } return a; } -- 2.39.5