From 021fd73d2b6dddc5c9ba9ca20c776099be34ef4f Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Mon, 4 Dec 2017 15:28:14 +0100 Subject: [PATCH] Update documentation. Replace all names mentioning the Laplace vector by the new approach of interpolating. Also update the description of how to ensure good convergence with the mapping, given that it was the Laplace smoothing in the old implementation that destroyed convergence rates. --- include/deal.II/fe/mapping_q_generic.h | 47 +++++++++----------- include/deal.II/grid/manifold.h | 59 ++++++++++++++------------ include/deal.II/grid/tria_accessor.h | 38 +++++++++++------ source/fe/mapping_q_generic.cc | 8 ++-- source/grid/tria.cc | 20 ++++----- source/grid/tria_accessor.cc | 17 ++++---- 6 files changed, 95 insertions(+), 94 deletions(-) diff --git a/include/deal.II/fe/mapping_q_generic.h b/include/deal.II/fe/mapping_q_generic.h index 7a5e14129e..0913ab053f 100644 --- a/include/deal.II/fe/mapping_q_generic.h +++ b/include/deal.II/fe/mapping_q_generic.h @@ -99,9 +99,10 @@ template class MappingQ; * transformation from the reference cell to such a cell would in general * contain inverted regions close to the boundary. * - * In order to avoid this situation, this class applies a smoothing on cells - * adjacent to the boundary by using so-called Laplace smoothing by - * default. In the algorithm that computes additional points, the + * In order to avoid this situation, this class applies an algorithm for + * making this transition smooth using a so-called transfinite interpolation + * that is essentially a linear blend between the descriptions along the + * surrounding entities. In the algorithm that computes 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. Points on objects higher up in the hierarchy are obtained @@ -113,24 +114,13 @@ template class MappingQ; * account when interpolating the position of the additional points inside the * quad and thus always result in a well-defined transformation. * - * 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 when switching from a - * curved manifold to a flat manifold over a face (and the associated - * smoothing). The finite element theory (see e.g. Strang and Fix, 1973, - * Sections 2.2 and 3.3 and in particular Theorem 3.6) requires the - * transformation to be globally C^0 continuous also over several elements and - * to be uniform as the mesh is refined. Even though the Laplace smoothing - * fixes the discontinuity within one layer of cells, it cannot provide - * uniformity as the change is always within one layer of elements only. 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 eventually covers a region of - * cells instead of only those that are immediately adjacent to the circular - * boundary. + * The interpolation scheme used in this class makes sure that curved + * descriptions can go over to flat descriptions within a single layer of + * elements, maintaining the overall optimal convergence rates of the finite + * element interpolation. However, one does often get better solution + * qualities if the transition between curved boundaries and flat interior + * domains is spread over a larger range as the mesh is refined. This is + * provided by the special manifold TransfiniteInterpolationManifold. * * @author Wolfgang Bangerth, 2015, Martin Kronbichler, 2017 */ @@ -664,19 +654,20 @@ protected: * 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. + * points in quads for 2d, in hexes for 3d) are then computed using an + * interpolation from the lower-dimensional entities (lines, quads) in order + * to make the transformation as smooth as possible without introducing + * additional boundary layers within the cells due to the placement of + * support points. * * 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 > diff --git a/include/deal.II/grid/manifold.h b/include/deal.II/grid/manifold.h index 06c7c77953..e34149a72d 100644 --- a/include/deal.II/grid/manifold.h +++ b/include/deal.II/grid/manifold.h @@ -71,12 +71,13 @@ namespace Manifolds * are the vertices, the line mid-points, and the face mid-points. * This results in a quadrature object with 26 (8+12+6) points. * - * The quadrature weights for these points are either chosen identically - * and equal to one over the number of quadrature points (if @p with_laplace - * is @p false), or in a way that gives points closer to the cell center - * (measured on the reference cell) a higher weight. These weights correspond - * to solving a Laplace equation and evaluating the solution at the quadrature - * points (if @p with_laplace is @p true). + * The quadrature weights for these points are either chosen identically and + * equal to one over the number of quadrature points (if @p + * with_interpolation is @p false), or in a way that gives points closer to + * the cell center (measured on the reference cell) a higher weight. These + * weights correspond to the weights applied to lines and vertices in + * transfinite interpolation (see TransfiniteInterpolationManifold for a + * more thorough description) if @p with_interpolation is @p true. * * The function is primarily used to construct the input argument * for the Manifold::get_new_point() function, which computes a new @@ -91,8 +92,8 @@ namespace Manifolds * * @param[in] iterator A mesh iterator that points to either a line, quad, * or hex. - * @param[in] with_laplace Whether or not to compute the quadrature weights - * by solving a Laplace equation, as discussed above. + * @param[in] with_interpolation Whether or not to compute the quadrature weights + * from transfinite interpolation, as discussed above. * @tparam MeshIteratorType An iterator type that corresponds to either * Triangulation::cell_iterator (or variants such as * Triangulation::active_cell_iterator or DoFHandler::cell_iterator) or @@ -103,7 +104,7 @@ namespace Manifolds DEAL_II_DEPRECATED Quadrature get_default_quadrature(const MeshIteratorType &iterator, - const bool with_laplace = false); + const bool with_interpolation = false); /** * Given a general mesh iterator, construct arrays of quadrature points and @@ -118,12 +119,13 @@ namespace Manifolds * are the vertices, the line mid-points, and the face mid-points. * This results in a points vector with 26 (8+12+6) points. * - * The quadrature weights for these points are either chosen identically - * and equal to one over the number of quadrature points (if @p with_laplace - * is @p false), or in a way that gives points closer to the cell center - * (measured on the reference cell) a higher weight. These weights correspond - * to solving a Laplace equation and evaluating the solution at the quadrature - * points (if @p with_laplace is @p true). + * The quadrature weights for these points are either chosen identically and + * equal to one over the number of quadrature points (if @p + * with_interpolation is @p false), or in a way that gives points closer to + * the cell center (measured on the reference cell) a higher weight. These + * weights correspond to the weights applied to lines and vertices in + * transfinite interpolation (see TransfiniteInterpolationManifold for a + * more thorough description) if @p with_interpolation is @p true. * * The function is primarily used to construct the input argument * for the Manifold::get_new_point() function, which computes a new @@ -138,8 +140,8 @@ namespace Manifolds * * @param[in] iterator A mesh iterator that points to either a line, quad, * or hex. - * @param[in] with_laplace Whether or not to compute the quadrature weights - * by solving a Laplace equation, as discussed above. + * @param[in] with_interpolation Whether or not to compute the quadrature weights + * from transfinite interpolation, as discussed above. * @tparam MeshIteratorType An iterator type that corresponds to either * Triangulation::cell_iterator (or variants such as * Triangulation::active_cell_iterator or DoFHandler::cell_iterator) or @@ -151,7 +153,7 @@ namespace Manifolds n_default_points_per_cell()>, std::array()> > get_default_points_and_weights(const MeshIteratorType &iterator, - const bool with_laplace = false); + const bool with_interpolation = false); } /** @@ -1139,9 +1141,9 @@ namespace Manifolds template Quadrature get_default_quadrature(const MeshIteratorType &iterator, - const bool with_laplace) + const bool with_interpolation) { - const auto points_and_weights = get_default_points_and_weights(iterator, with_laplace); + const auto points_and_weights = get_default_points_and_weights(iterator, with_interpolation); static const int spacedim = MeshIteratorType::AccessorType::space_dimension; return Quadrature (std::vector>(points_and_weights.first.begin(), @@ -1157,7 +1159,7 @@ namespace Manifolds n_default_points_per_cell()>, std::array()> > get_default_points_and_weights(const MeshIteratorType &iterator, - const bool with_laplace) + const bool with_interpolation) { const int dim = MeshIteratorType::AccessorType::structure_dimension; const int spacedim = MeshIteratorType::AccessorType::space_dimension; @@ -1170,7 +1172,7 @@ namespace Manifolds // note that the exact weights are chosen such as to minimize the // distortion of the four new quads from the optimal shape; their // derivation and values is copied over from the - // @p{MappingQ::set_laplace_on_vector} function + // interpolation function in the mapping switch (dim) { case 1: @@ -1193,7 +1195,7 @@ namespace Manifolds iterator->line(i)->get_manifold().get_new_point_on_line(iterator->line(i)) ); } - if (with_laplace) + if (with_interpolation) { std::fill(points_weights.second.begin(), points_weights.second.begin()+4, -0.25); std::fill(points_weights.second.begin()+4, points_weights.second.end(), 0.5); @@ -1217,9 +1219,10 @@ namespace Manifolds unsigned int j=0; // note that the exact weights are chosen such as to minimize the - // distortion of the eight new hexes from the optimal shape; their - // derivation and values is copied over from the - // @p{MappingQ::set_laplace_on_vector} function + // distortion of the eight new hexes from the optimal shape through + // transfinite interpolation from the faces and vertices, see + // TransfiniteInterpolationManifold for a deeper explanation of the + // mechanisms for (unsigned int i=0; i::vertices_per_cell; ++i, ++j) { (*sp3)[j] = hex->vertex(i); @@ -1240,8 +1243,8 @@ namespace Manifolds points_weights.second[j] = 1.0/2.0; } // Overwrite the weights with 1/np if we don't want to use - // laplace vectors. - if (with_laplace == false) + // interpolation. + if (with_interpolation == false) std::fill(points_weights.second.begin(), points_weights.second.end(), 1.0/np); } break; diff --git a/include/deal.II/grid/tria_accessor.h b/include/deal.II/grid/tria_accessor.h index 22e9fd1aba..2b14a773e4 100644 --- a/include/deal.II/grid/tria_accessor.h +++ b/include/deal.II/grid/tria_accessor.h @@ -1394,19 +1394,29 @@ public: * point to the physical space again; the resulting point is guaranteed to * lie within the manifold, even if the manifold is curved. * - * When the geometry of a TriaAccessor is not flat, or when part of the - * bounding objects of this TriaAccessor are not flat, the result given by - * the TriaAccessor::center() function may not be accurate enough, even when - * parameter @p respect_manifold is set to true. If you find this to be - * case, than you can further refine the computation of the center by - * setting to true the second additional parameter - * @p use_laplace_transformation, which will force this function to compute - * the location of the center by solving a linear elasticity problem with - * Dirichlet boundary conditions set to the location of the bounding - * vertices and the centers of the bounding lines and quads. + * When the object uses a different manifold description as its surrounding, + * like when part of the bounding objects of this TriaAccessor use a + * non-flat manifold description but the object itself is flat, the result + * given by the TriaAccessor::center() function may not be accurate enough, + * even when parameter @p respect_manifold is set to true. If you find this + * to be case, than you can further refine the computation of the center by + * setting to true the second additional parameter @p + * interpolate_from_surrounding. This computes the location of the center by + * a so-called transfinite interpolation from the center of all the bounding + * objects. For a 2D object, it puts a weight of 1/2 on each of + * the four surrounding lines and a weight -1/4 on the four + * vertices. This corresponds to a linear interpolation between the + * descriptions of the four faces, subtracting the contribution of the + * vertices that is added twice when coming through both lines adjacent to + * the vertex. In 3D, the weights for faces are 1/2, the + * weights for lines are -1/4, and the weights for vertices are + * 1/8. For further information, also confer to the + * TransfiniteInterpolationManifold class that is able to not only apply + * this beneficial description to a single cell but all children of a coarse + * cell. */ Point center (const bool respect_manifold=false, - const bool use_laplace_transformation=false) const; + const bool interpolate_from_surrounding=false) const; /** * Return the barycenter (also called centroid) @@ -1793,12 +1803,12 @@ public: /** * Return the center of this object, which of course coincides with the * location of the vertex this object refers to. The parameters @p - * respect_manifold and @p use_laplace_transformation are not used. They are - * there to provide the same interface as + * respect_manifold and @p interpolate_from_surrounding are not used. They + * are there to provide the same interface as * TriaAccessor. */ Point center (const bool respect_manifold=false, - const bool use_laplace_transformation=false) const; + const bool interpolate_from_surrounding=false) const; /** * Compute the dim-dimensional measure of the object. For a dim-dimensional diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc index 29e35101aa..9c2647fd0d 100644 --- a/source/fe/mapping_q_generic.cc +++ b/source/fe/mapping_q_generic.cc @@ -925,8 +925,8 @@ namespace internal * This function is needed by the constructor of * MappingQ for dim= 2 and 3. * - * For the definition of the @p support_point_weights_on_quad please refer to - * equation (8) of the `mapping' report. + * For the definition of the @p support_point_weights_on_quad please + * refer to the description of TransfiniteInterpolationManifold. */ dealii::Table<2,double> compute_support_point_weights_on_quad(const unsigned int polynomial_degree) @@ -974,8 +974,8 @@ namespace internal /** * This function is needed by the constructor of MappingQ<3>. * - * For the definition of the @p support_point_weights_on_hex please refer to - * equation (8) of the `mapping' report. + * For the definition of the @p support_point_weights_on_quad please + * refer to the description of TransfiniteInterpolationManifold. */ dealii::Table<2,double> compute_support_point_weights_on_hex(const unsigned int polynomial_degree) diff --git a/source/grid/tria.cc b/source/grid/tria.cc index 7488644353..8ff8460535 100644 --- a/source/grid/tria.cc +++ b/source/grid/tria.cc @@ -4002,17 +4002,13 @@ namespace internal // if the user_flag is set, i.e. if the cell is at the // boundary, use a different calculation of the middle vertex - // here. this is of advantage, if the boundary is strongly - // curved and the cell has a high aspect ratio. this can - // happen for example, if it was refined anisotropically - // before. + // here. this is of advantage if the boundary is strongly + // curved (whereas the cell is not) and the cell has a high + // aspect ratio. if (cell->user_flag_set()) { - // first reset the user_flag + // first reset the user_flag and then refine cell->clear_user_flag(); - - // the user flag indicates: at least one face is at the - // boundary. triangulation.vertices[next_unused_vertex] =cell->center(true, true); } @@ -5793,7 +5789,7 @@ namespace internal // isotropic refinement Assert(quad_ref_case==RefinementCase::no_refinement, ExcInternalError()); - // set the middle vertex appropriately it might be that + // set the middle vertex appropriately: it might be that // the quad itself is not at the boundary, but that one of // its lines actually is. in this case, the newly created // vertices at the centers of the lines are not @@ -5808,9 +5804,9 @@ namespace internal // // note that the exact weights are chosen such as to // minimize the distortion of the four new quads from the - // optimal shape; their derivation and values is copied - // over from the @p{MappingQ::set_laplace_on_vector} - // function + // optimal shape. their description uses the formulas + // underlying the TransfiniteInterpolationManifold + // implementation triangulation.vertices[next_unused_vertex] = quad->center(true, true); triangulation.vertices_used[next_unused_vertex] = true; diff --git a/source/grid/tria_accessor.cc b/source/grid/tria_accessor.cc index 4ea9917566..9b655f2d44 100644 --- a/source/grid/tria_accessor.cc +++ b/source/grid/tria_accessor.cc @@ -1045,17 +1045,18 @@ namespace template Point get_new_point_on_object(const TriaAccessor &obj, - const bool use_laplace) + const bool use_interpolation) { - // if we should not do mesh smoothing or if the object is of the old - // Boundary type, do only use the old points - if (use_laplace == false || + // if we should not interpolate from the surroundings of the current + // object or if the object is of the old Boundary type, only use the + // vertices points + if (use_interpolation == false || dynamic_cast *>(&obj.get_manifold()) != nullptr) return get_new_point_on_object(obj); else { TriaRawIterator > it(obj); - const auto points_and_weights = Manifolds::get_default_points_and_weights(it, use_laplace); + const auto points_and_weights = Manifolds::get_default_points_and_weights(it, use_interpolation); return obj.get_manifold().get_new_point(make_array_view(points_and_weights.first.begin(), points_and_weights.first.end()), make_array_view(points_and_weights.second.begin(), @@ -1393,18 +1394,18 @@ TriaAccessor template Point TriaAccessor::center (const bool respect_manifold, - const bool use_laplace) const + const bool use_interpolation) const { if (respect_manifold == false) { - Assert(use_laplace == false, ExcNotImplemented()); + Assert(use_interpolation == false, ExcNotImplemented()); Point p; for (unsigned int v=0; v::vertices_per_cell; ++v) p += vertex(v); return p/GeometryInfo::vertices_per_cell; } else - return get_new_point_on_object(*this, use_laplace); + return get_new_point_on_object(*this, use_interpolation); } -- 2.39.5