]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Update documentation.
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Mon, 4 Dec 2017 14:28:14 +0000 (15:28 +0100)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Tue, 12 Dec 2017 08:58:05 +0000 (09:58 +0100)
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
include/deal.II/grid/manifold.h
include/deal.II/grid/tria_accessor.h
source/fe/mapping_q_generic.cc
source/grid/tria.cc
source/grid/tria_accessor.cc

index 7a5e14129ebfbf50cf48261997063fc334023dc4..0913ab053f0e9e968319feaedec77da87d438599 100644 (file)
@@ -99,9 +99,10 @@ template <int,int> 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 <int,int> 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<spacedim, it asks the underlying
-   * manifold for the locations of interior points.
+   * interpolation from the surrounding points using weights for transfinite
+   * interpolation, or if dim<spacedim, it asks the underlying manifold for
+   * the locations of interior points.
    */
   virtual
   std::vector<Point<spacedim> >
index 06c7c7795364d74f2972567824a915f2b7f75fa7..e34149a72dea8089f388a37e54a3b207b0ebc76e 100644 (file)
@@ -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<MeshIteratorType::AccessorType::space_dimension>
   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<MeshIteratorType>()>,
       std::array<double, n_default_points_per_cell<MeshIteratorType>()> >
       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 <typename MeshIteratorType>
   Quadrature<MeshIteratorType::AccessorType::space_dimension>
   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<spacedim>
            (std::vector<Point<spacedim>>(points_and_weights.first.begin(),
@@ -1157,7 +1159,7 @@ namespace Manifolds
       n_default_points_per_cell<MeshIteratorType>()>,
       std::array<double, n_default_points_per_cell<MeshIteratorType>()> >
       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<GeometryInfo<dim>::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;
index 22e9fd1aba559883567cede7b368bdbc3abb400a..2b14a773e41bfd9b6746abff723b3c5782b3d08b 100644 (file)
@@ -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 <code>1/2</code> on each of
+   * the four surrounding lines and a weight <code>-1/4</code> 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 <code>1/2</code>, the
+   * weights for lines are <code>-1/4</code>, and the weights for vertices are
+   * <code>1/8</code>. 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<spacedim> 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
    * <code>TriaAccessor<structdim,dim,spacedim></code>.
    */
   Point<spacedim> 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
index 29e35101aa4fa00c72d8cbae4ffe5520de557ce7..9c2647fd0db67ca0085cd117caef66c420c6c3e7 100644 (file)
@@ -925,8 +925,8 @@ namespace internal
        * This function is needed by the constructor of
        * <tt>MappingQ<dim,spacedim></tt> for <tt>dim=</tt> 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 <tt>MappingQ<3></tt>.
        *
-       * 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)
index 7488644353ab70540fd5ff75526741bb3da5f277..8ff8460535464013e2ffa2f8e46ccbd5190d82d6 100644 (file)
@@ -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<dim-1>::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;
index 4ea9917566a7e3f46745d75b6f3e461147fecfb8..9b655f2d442d07aaba47d0a8bb635d860faa4f6a 100644 (file)
@@ -1045,17 +1045,18 @@ namespace
 
   template <int structdim, int dim, int spacedim>
   Point<spacedim> get_new_point_on_object(const TriaAccessor<structdim, dim, spacedim> &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<const Boundary<dim,spacedim> *>(&obj.get_manifold()) != nullptr)
       return get_new_point_on_object(obj);
     else
       {
         TriaRawIterator<TriaAccessor<structdim, dim, spacedim> > 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<structdim, dim, spacedim>
 template <int structdim, int dim, int spacedim>
 Point<spacedim>
 TriaAccessor<structdim, dim, spacedim>::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<spacedim> p;
       for (unsigned int v=0; v<GeometryInfo<structdim>::vertices_per_cell; ++v)
         p += vertex(v);
       return p/GeometryInfo<structdim>::vertices_per_cell;
     }
   else
-    return get_new_point_on_object(*this, use_laplace);
+    return get_new_point_on_object(*this, use_interpolation);
 }
 
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.