]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Unify interpolation from surrounding points.
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Fri, 17 Mar 2017 14:21:31 +0000 (15:21 +0100)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Sun, 19 Mar 2017 09:17:25 +0000 (10:17 +0100)
Query multiple points from mapping.

include/deal.II/fe/mapping_q.h
include/deal.II/fe/mapping_q_generic.h
source/fe/mapping_q.cc
source/fe/mapping_q_generic.cc

index 2694e8521d040c4fa60ec7749bdc4f5feb5247ea..dd6946c3b69035f6992394696dd3b864b5c9194e 100644 (file)
@@ -72,13 +72,14 @@ DEAL_II_NAMESPACE_OPEN
  * <h4>Behavior along curved boundaries and with different manifolds</h4>
  *
  * 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.
index 1832e377865186f47378bbe442565bb79f061460..4c604c5c45a036ba4e8e150b8b3fe880da3fa4f8 100644 (file)
@@ -84,13 +84,14 @@ template <int,int> class MappingQ;
  * <h4>Behavior along curved boundaries and with different manifolds</h4>
  *
  * 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 <int,int> 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 <int dim, int spacedim=dim>
 class MappingQGeneric : public Mapping<dim,spacedim>
@@ -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<dim> > 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, <code>4 + 4*(degree-1)</code> for
+   * a quad, <code>8 + 12*(degree-1) + 6*(degree-1)*(degree-1)</code> 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<Table<2,double> > 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<Table<2,double> > 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., <code>(degree+1)^dim - 2^dim</code>.
    */
-  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
index 87ed327bb8d5251a7fac24b7e0fd9798924f80bb..4b0b849ac7b8a32bb66c42adbeed7b49dc3f40fd 100644 (file)
@@ -55,7 +55,7 @@ MappingQ<dim,spacedim>::InternalData::memory_consumption () const
 template<int dim, int spacedim>
 MappingQ<dim,spacedim>::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<dim,spacedim>::MappingQ (const unsigned int degree,
   // created via the shared_ptr objects
   qp_mapping (this->polynomial_degree>1
               ?
-              std_cxx11::shared_ptr<const MappingQGeneric<dim,spacedim> >(new MappingQGeneric<dim,spacedim>(degree, smooth_support_points_on_flat_manifold))
+              std_cxx11::shared_ptr<const MappingQGeneric<dim,spacedim> >(new MappingQGeneric<dim,spacedim>(degree, smooth_support_points))
               :
               q1_mapping)
 {}
index 27a1fbb99f0bf42d043f028e89be7301519015b4..4392dd8551fc32beb5219d3413ea4956a6b784e5 100644 (file)
@@ -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<dim>
-        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<int dim, int spacedim>
 MappingQGeneric<dim,spacedim>::InternalData::InternalData (const unsigned int polynomial_degree)
   :
   polynomial_degree (polynomial_degree),
-  n_shape_functions (Utilities::fixed_power<dim>(polynomial_degree+1)),
-  line_support_points (polynomial_degree + 1)
+  n_shape_functions (Utilities::fixed_power<dim>(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<int dim>
   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 <int dim>
   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<dim>(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<dim>(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<Table<2,double> >
+  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<Table<2,double> > 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<polynomial_degree-1; ++q)
+      for (unsigned int i=0; i<GeometryInfo<1>::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<Table<2,double> >
+  compute_support_point_weights_interior(const unsigned int polynomial_degree,
+                                         const unsigned int dim)
+  {
+    Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim));
+    std::vector<Table<2,double> > 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<polynomial_degree-1; ++q)
+      for (unsigned int i=0; i<GeometryInfo<1>::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<polynomial_degree-1; ++q2)
+          for (unsigned int q1=0; q1<polynomial_degree-1; ++q1, ++q)
+            {
+              Point<2> point(quadrature.point(q1+1)[0],
+                             quadrature.point(q2+1)[0]);
+              for (unsigned int i=0; i<GeometryInfo<2>::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<polynomial_degree-1; ++q3)
+          for (unsigned int q2=0; q2<polynomial_degree-1; ++q2)
+            for (unsigned int q1=0; q1<polynomial_degree-1; ++q1, ++q)
+              {
+                Point<3> point(quadrature.point(q1+1)[0],
+                               quadrature.point(q2+1)[0],
+                               quadrature.point(q3+1)[0]);
+                for (unsigned int i=0; i<GeometryInfo<3>::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 <int dim>
+  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<dim> fe(polynomial_degree);
+    std::vector<Point<dim> > unit_support_points = fe.get_unit_support_points();
+
+    Table<2,double> output(fe.dofs_per_cell - GeometryInfo<dim>::vertices_per_cell,
+                           GeometryInfo<dim>::vertices_per_cell);
+    for (unsigned int q=0; q<output.size(0); ++q)
+      for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+        output(q,i) = GeometryInfo<dim>::d_linear_shape_function(unit_support_points[q+GeometryInfo<dim>::vertices_per_cell],
+                                                                 i);
+
+    return output;
+  }
 }
 
 
@@ -1207,14 +1306,15 @@ namespace
 
 template<int dim, int spacedim>
 MappingQGeneric<dim,spacedim>::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<dim>(this->polynomial_degree) : 0),
-  support_point_weights_on_quad (compute_support_point_weights_on_quad<dim>(this->polynomial_degree)),
-  support_point_weights_on_hex (compute_support_point_weights_on_hex<dim>(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<dim>(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<int dim, int spacedim>
 MappingQGeneric<dim,spacedim>::MappingQGeneric (const MappingQGeneric<dim,spacedim> &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<dim>(*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<const  Tensor<3,dim> >                  &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<int dim, int spacedim>
-  void
-  get_intermediate_points (const Manifold<dim, spacedim> &manifold,
-                           const QGaussLobatto<1>        &line_support_points,
-                           const std::vector<Point<spacedim> > &surrounding_points,
-                           std::vector<Point<spacedim> > &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<double> 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<n; ++i)
-          {
-            const double x = line_support_points.point(i+1)[0];
-            w[1] = x;
-            w[0] = (1-x);
-            points[i] = manifold.get_new_point(surrounding_points, w);
-          }
-        break;
-      }
-
-      case 4:
-      {
-        Assert(spacedim >= 2, ExcImpossibleInDim(spacedim));
-        const unsigned m=
-          static_cast<unsigned int>(std::sqrt(static_cast<double>(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<m; ++i)
-          {
-            p[1] = line_support_points.point(1+i)[0];
-            for (unsigned int j=0; j<m; ++j, ++c)
-              {
-                p[0] = line_support_points.point(1+j)[0];
-
-                for (unsigned int l=0; l<4; ++l)
-                  w[l] = GeometryInfo<2>::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<m; ++k)
-          {
-            p[2] = line_support_points.point(1+k)[0];
-            for (unsigned int i=0; i<m; ++i)
-              {
-                p[1] = line_support_points.point(1+i)[0];
-                for (unsigned int j=0; j<m; ++j, ++c)
-                  {
-                    p[0] = line_support_points.point(1+j)[0];
-
-                    for (unsigned int l=0; l<8; ++l)
-                      w[l] = GeometryInfo<3>::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<dim,spacedim>.
    */
   template <int dim, int spacedim, class TriaIterator>
-  void get_intermediate_points_on_object(const Manifold<dim, spacedim> &manifold,
-                                         const QGaussLobatto<1>        &line_support_points,
-                                         const TriaIterator &iter,
-                                         std::vector<Point<spacedim> > &points)
+  void add_intermediate_points(const Manifold<dim, spacedim> &manifold,
+                               const Table<2,double>         &weight_table,
+                               const TriaIterator            &iter,
+                               std::vector<Point<spacedim> > &points)
   {
     const unsigned int structdim = TriaIterator::AccessorType::structure_dimension;
 
     // Try backward compatibility option.
-    if (const Boundary<dim,spacedim> *boundary
-        = dynamic_cast<const Boundary<dim,spacedim> *>(&manifold))
+    const Boundary<dim,spacedim> *boundary
+      = dynamic_cast<const Boundary<dim,spacedim> *>(&manifold);
+    if (structdim < spacedim && boundary != NULL)
       // This is actually a boundary. Call old methods.
       {
+        std::vector<Point<spacedim> > new_points(weight_table.size(0));
         switch (structdim)
           {
           case 1:
           {
             const typename Triangulation<dim,spacedim>::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<dim,spacedim>::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<new_points.size(); ++i)
+          points.push_back(new_points[i]);
       }
     else
       {
         std::vector<Point<spacedim> > sp(GeometryInfo<structdim>::vertices_per_cell);
         for (unsigned int i=0; i<sp.size(); ++i)
           sp[i] = iter->vertex(i);
-        get_intermediate_points(manifold, line_support_points, sp, points);
+        manifold.add_new_points(sp, weight_table, points);
       }
   }
 
 
-  /**
-   * Take a <tt>support_point_weights_on_hex(quad)</tt> 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 <int spacedim>
-  void add_weighted_interior_points(const Table<2,double>   &lvs,
-                                    std::vector<Point<spacedim> > &a)
+
+  template <typename Iterator>
+  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<n_inner_apply; ++unit_point)
-      {
-        Assert(lvs.n_cols()==n_outer_apply, ExcInternalError());
-        Point<spacedim> p;
-        for (unsigned int k=0; k<n_outer_apply; ++k)
-          p+=lvs[unit_point][k]*a[k];
+    Assert(false, ExcNotImplemented());
+    return true;
+  }
 
-        a.push_back(p);
-      }
+  bool check_identical_manifolds_of_quads(const Triangulation<3,3>::cell_iterator &cell)
+  {
+    for (unsigned int f=0; f<GeometryInfo<3>::faces_per_cell; ++f)
+      if (&cell->face(f)->get_manifold() != &cell->get_manifold())
+        return false;
+    return true;
   }
 }
 
 
+
 template <int dim, int spacedim>
 void
 MappingQGeneric<dim,spacedim>::
@@ -3749,7 +3736,9 @@ add_line_support_points (const typename Triangulation<dim,spacedim>::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<Point<3> >                &a) const
+                        std::vector<Point<3> >                  &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<Point<3> > quad_points ((polynomial_degree-1)*(polynomial_degree-1));
+  std::vector<Point<3> > quad_points;
   std::vector<Point<3> > 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<GeometryInfo<3>::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<const Boundary<3,3> *>(&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<Point<3> >                &a) const
+                        std::vector<Point<3> >                  &a) const
 {
-  std::vector<Point<3> > 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<quad_points.size(); ++i)
-    a.push_back(quad_points[i]);
+  add_intermediate_points (cell->get_manifold(), support_point_weights_interior[1],
+                           cell, a);
 }
 
 
@@ -3893,72 +3883,83 @@ MappingQGeneric<dim,spacedim>::
 compute_mapping_support_points(const typename Triangulation<dim,spacedim>::cell_iterator &cell) const
 {
   // get the vertices first
-  std::vector<Point<spacedim> > a(GeometryInfo<dim>::vertices_per_cell);
+  std::vector<Point<spacedim> > a;
+  a.reserve(Utilities::fixed_power<dim>(polynomial_degree+1));
   for (unsigned int i=0; i<GeometryInfo<dim>::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<Point<spacedim> > quad_points (Utilities::fixed_power<dim>(polynomial_degree-1));
-            get_intermediate_points_on_object(cell->get_manifold(), line_support_points, cell, quad_points);
-            for (unsigned int i=0; i<quad_points.size(); ++i)
-              a.push_back(quad_points[i]);
-          }
-        else
-          add_weighted_interior_points (support_point_weights_on_quad, a);
-        break;
+  if (this->polynomial_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<GeometryInfo<dim>::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<Point<spacedim> > 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<Point<spacedim> > hex_points (Utilities::fixed_power<dim>(polynomial_degree-1));
-            get_intermediate_points_on_object(cell->get_manifold(), line_support_points, cell, hex_points);
-            for (unsigned int i=0; i<hex_points.size(); ++i)
-              a.push_back(hex_points[i]);
-          }
-        else
-          add_weighted_interior_points (support_point_weights_on_hex, a);
-        break;
-      }
+      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. decide whether to choose
+          // the smoothed variant that include all points and goes directly
+          // through the manifold or only the vertex points.
+          if (dim != spacedim)
+            add_quad_support_points(cell, a);
+          else if (smooth_support_points == false)
+            add_intermediate_points(cell->get_manifold(), support_point_weights_interior[1],
+                                    cell, a);
+          else
+            {
+              std::vector<Point<spacedim> > 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<Point<spacedim> > 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;
 }

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.