]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Allow to compute affine approximation on arbitrary set of points
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Fri, 2 Oct 2020 10:14:26 +0000 (12:14 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Sat, 3 Oct 2020 19:13:58 +0000 (21:13 +0200)
include/deal.II/grid/grid_tools.h
source/grid/grid_tools.cc
source/grid/grid_tools.inst.in
source/grid/tria_accessor.cc

index d781ceba28ec6a671ccb3d7c87236da8582440d8..f7a07d15acb3ba982c0dc812061fc549a035a0e7 100644 (file)
@@ -247,6 +247,32 @@ namespace GridTools
   double
   cell_measure(const T &, ...);
 
+  /**
+   * This function computes an affine approximation of the map from the unit
+   * coordinates to the real coordinates of the form $p_\text{real} = A
+   * p_\text{unit} + b $ by a least squares fit of this affine function to the
+   * $2^\text{dim}$ vertices representing a quadrilateral or hexahedral cell
+   * in `spacedim` dimensions. The result is returned as a pair with the
+   * matrix <i>A</i> as the first argument and the vector <i>b</i> describing
+   * distance of the plane to the origin.
+   *
+   * For any valid mesh cell whose geometry is not degenerate, this operation
+   * results in a unique affine mapping, even in cases where the actual
+   * transformation by a bi-/trilinear or higher order mapping might be
+   * singular. The result is exact in case the transformation from the unit to
+   * the real cell is indeed affine, such as in one dimension or for Cartesian
+   * and affine (parallelogram) meshes in 2D/3D.
+   *
+   * This approximation is underlying the function
+   * TriaAccessor::real_to_unit_cell_affine_approximation() function.
+   *
+   * For exact transformations to the unit cell, use
+   * Mapping::transform_real_to_unit_cell().
+   */
+  template <int dim, int spacedim>
+  std::pair<DerivativeForm<1, dim, spacedim>, Tensor<1, spacedim>>
+  affine_cell_approximation(const ArrayView<const Point<spacedim>> &vertices);
+
   /**
    * Computes an aspect ratio measure for all locally-owned active cells and
    * fills a vector with one entry per cell, given a @p triangulation and
index 62212f88816fe1f5ab1c0bf3486bf8c434aff0fa..2499c4f2721bf0ead3f80a151e37bde3f6c0d3cb 100644 (file)
@@ -188,6 +188,122 @@ namespace GridTools
 
 
 
+  namespace
+  {
+    /**
+     * The algorithm to compute the affine approximation to a cell finds an
+     * affine map A x_hat + b from the reference cell to the real space.
+     *
+     * Some details about how we compute the least square plane. We look
+     * for a spacedim x (dim + 1) matrix X such that X * M = Y where M is
+     * a (dim+1) x n_vertices matrix and Y a spacedim x n_vertices.  And:
+     * The i-th column of M is unit_vertex[i] and the last row all
+     * 1's. The i-th column of Y is real_vertex[i].  If we split X=[A|b],
+     * the least square approx is A x_hat+b Classically X = Y * (M^t (M
+     * M^t)^{-1}) Let K = M^t * (M M^t)^{-1} = [KA Kb] this can be
+     * precomputed, and that is exactly what we do.  Finally A = Y*KA and
+     * b = Y*Kb.
+     */
+    template <int dim>
+    struct TransformR2UAffine
+    {
+      static const double KA[GeometryInfo<dim>::vertices_per_cell][dim];
+      static const double Kb[GeometryInfo<dim>::vertices_per_cell];
+    };
+
+
+    /*
+      Octave code:
+      M=[0 1; 1 1];
+      K1 = transpose(M) * inverse (M*transpose(M));
+      printf ("{%f, %f},\n", K1' );
+    */
+    template <>
+    const double TransformR2UAffine<1>::KA[GeometryInfo<1>::vertices_per_cell]
+                                          [1] = {{-1.000000}, {1.000000}};
+
+    template <>
+    const double TransformR2UAffine<1>::Kb[GeometryInfo<1>::vertices_per_cell] =
+      {1.000000, 0.000000};
+
+
+    /*
+      Octave code:
+      M=[0 1 0 1;0 0 1 1;1 1 1 1];
+      K2 = transpose(M) * inverse (M*transpose(M));
+      printf ("{%f, %f, %f},\n", K2' );
+    */
+    template <>
+    const double TransformR2UAffine<2>::KA[GeometryInfo<2>::vertices_per_cell]
+                                          [2] = {{-0.500000, -0.500000},
+                                                 {0.500000, -0.500000},
+                                                 {-0.500000, 0.500000},
+                                                 {0.500000, 0.500000}};
+
+    /*
+      Octave code:
+      M=[0 1 0 1 0 1 0 1;0 0 1 1 0 0 1 1; 0 0 0 0 1 1 1 1; 1 1 1 1 1 1 1 1];
+      K3 = transpose(M) * inverse (M*transpose(M))
+      printf ("{%f, %f, %f, %f},\n", K3' );
+    */
+    template <>
+    const double TransformR2UAffine<2>::Kb[GeometryInfo<2>::vertices_per_cell] =
+      {0.750000, 0.250000, 0.250000, -0.250000};
+
+
+    template <>
+    const double TransformR2UAffine<3>::KA[GeometryInfo<3>::vertices_per_cell]
+                                          [3] = {
+                                            {-0.250000, -0.250000, -0.250000},
+                                            {0.250000, -0.250000, -0.250000},
+                                            {-0.250000, 0.250000, -0.250000},
+                                            {0.250000, 0.250000, -0.250000},
+                                            {-0.250000, -0.250000, 0.250000},
+                                            {0.250000, -0.250000, 0.250000},
+                                            {-0.250000, 0.250000, 0.250000},
+                                            {0.250000, 0.250000, 0.250000}
+
+    };
+
+
+    template <>
+    const double TransformR2UAffine<3>::Kb[GeometryInfo<3>::vertices_per_cell] =
+      {0.500000,
+       0.250000,
+       0.250000,
+       0.000000,
+       0.250000,
+       0.000000,
+       0.000000,
+       -0.250000};
+  } // namespace
+
+
+
+  template <int dim, int spacedim>
+  std::pair<DerivativeForm<1, dim, spacedim>, Tensor<1, spacedim>>
+  affine_cell_approximation(const ArrayView<const Point<spacedim>> &vertices)
+  {
+    AssertDimension(vertices.size(), GeometryInfo<dim>::vertices_per_cell);
+
+    // A = vertex * KA
+    DerivativeForm<1, dim, spacedim> A;
+
+    for (unsigned int d = 0; d < spacedim; ++d)
+      for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
+        for (unsigned int e = 0; e < dim; ++e)
+          A[d][e] += vertices[v][d] * TransformR2UAffine<dim>::KA[v][e];
+
+    // b = vertex * Kb
+    Tensor<1, spacedim> b;
+    for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
+      b += vertices[v] * TransformR2UAffine<dim>::Kb[v];
+
+    return std::make_pair(A, b);
+  }
+
+
+
   template <int dim>
   Vector<double>
   compute_aspect_ratio_of_cells(const Mapping<dim> &      mapping,
index fb09c045f64c46df78bd981e411657aca33883c7..713d8c77bda1ee8c305c1321ab460ba45b64374f 100644 (file)
@@ -204,6 +204,12 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
       volume(const Triangulation<deal_II_dimension, deal_II_space_dimension> &,
              const Mapping<deal_II_dimension, deal_II_space_dimension> &);
 
+      template std::pair<
+        DerivativeForm<1, deal_II_dimension, deal_II_space_dimension>,
+        Tensor<1, deal_II_space_dimension>>
+      affine_cell_approximation<deal_II_dimension, deal_II_space_dimension>(
+        const ArrayView<const Point<deal_II_space_dimension>> &);
+
       template BoundingBox<deal_II_space_dimension>
       compute_bounding_box(
         const Triangulation<deal_II_dimension, deal_II_space_dimension> &);
index 80ac4c3c49aca749e18c4b17a2c20dd96e7b6810..002031c0b3148627270fc8aa7fcf58434f087b78 100644 (file)
@@ -1700,131 +1700,26 @@ TriaAccessor<structdim, dim, spacedim>::intermediate_point(
 }
 
 
-namespace
-{
-  /**
-   * The algorithm to compute the affine approximation to the point on the
-   * unit cell does the following steps:
-   * <ul>
-   * <li> find the least square dim-dimensional plane approximating the cell
-   * vertices, i.e. we find an affine map A x_hat + b from the reference cell
-   * to the real space.
-   * <li> Solve the equation A x_hat + b = p for x_hat
-   * </ul>
-   *
-   * Some details about how we compute the least square plane. We look
-   * for a spacedim x (dim + 1) matrix X such that X * M = Y where M is
-   * a (dim+1) x n_vertices matrix and Y a spacedim x n_vertices.  And:
-   * The i-th column of M is unit_vertex[i] and the last row all
-   * 1's. The i-th column of Y is real_vertex[i].  If we split X=[A|b],
-   * the least square approx is A x_hat+b Classically X = Y * (M^t (M
-   * M^t)^{-1}) Let K = M^t * (M M^t)^{-1} = [KA Kb] this can be
-   * precomputed, and that is exactly what we do.  Finally A = Y*KA and
-   * b = Y*Kb.
-   */
-  template <int dim>
-  struct TransformR2UAffine
-  {
-    static const double KA[GeometryInfo<dim>::vertices_per_cell][dim];
-    static const double Kb[GeometryInfo<dim>::vertices_per_cell];
-  };
-
-
-  /*
-    Octave code:
-    M=[0 1; 1 1];
-    K1 = transpose(M) * inverse (M*transpose(M));
-    printf ("{%f, %f},\n", K1' );
-  */
-  template <>
-  const double TransformR2UAffine<1>::KA[GeometryInfo<1>::vertices_per_cell]
-                                        [1] = {{-1.000000}, {1.000000}};
-
-  template <>
-  const double TransformR2UAffine<1>::Kb[GeometryInfo<1>::vertices_per_cell] =
-    {1.000000, 0.000000};
-
-
-  /*
-    Octave code:
-    M=[0 1 0 1;0 0 1 1;1 1 1 1];
-    K2 = transpose(M) * inverse (M*transpose(M));
-    printf ("{%f, %f, %f},\n", K2' );
-  */
-  template <>
-  const double TransformR2UAffine<2>::KA[GeometryInfo<2>::vertices_per_cell]
-                                        [2] = {{-0.500000, -0.500000},
-                                               {0.500000, -0.500000},
-                                               {-0.500000, 0.500000},
-                                               {0.500000, 0.500000}};
-
-  /*
-    Octave code:
-    M=[0 1 0 1 0 1 0 1;0 0 1 1 0 0 1 1; 0 0 0 0 1 1 1 1; 1 1 1 1 1 1 1 1];
-    K3 = transpose(M) * inverse (M*transpose(M))
-    printf ("{%f, %f, %f, %f},\n", K3' );
-  */
-  template <>
-  const double TransformR2UAffine<2>::Kb[GeometryInfo<2>::vertices_per_cell] =
-    {0.750000, 0.250000, 0.250000, -0.250000};
-
-
-  template <>
-  const double TransformR2UAffine<3>::KA[GeometryInfo<3>::vertices_per_cell]
-                                        [3] = {
-                                          {-0.250000, -0.250000, -0.250000},
-                                          {0.250000, -0.250000, -0.250000},
-                                          {-0.250000, 0.250000, -0.250000},
-                                          {0.250000, 0.250000, -0.250000},
-                                          {-0.250000, -0.250000, 0.250000},
-                                          {0.250000, -0.250000, 0.250000},
-                                          {-0.250000, 0.250000, 0.250000},
-                                          {0.250000, 0.250000, 0.250000}
-
-  };
-
-
-  template <>
-  const double TransformR2UAffine<3>::Kb[GeometryInfo<3>::vertices_per_cell] = {
-    0.500000,
-    0.250000,
-    0.250000,
-    0.000000,
-    0.250000,
-    0.000000,
-    0.000000,
-    -0.250000};
-} // namespace
-
 
 template <int structdim, int dim, int spacedim>
 Point<structdim>
 TriaAccessor<structdim, dim, spacedim>::real_to_unit_cell_affine_approximation(
   const Point<spacedim> &point) const
 {
-  // A = vertex * KA
-  DerivativeForm<1, structdim, spacedim> A;
-
-  // copy vertices to avoid expensive resolution of vertex index inside loop
   std::array<Point<spacedim>, GeometryInfo<structdim>::vertices_per_cell>
     vertices;
   for (const unsigned int v : this->vertex_indices())
     vertices[v] = this->vertex(v);
-  for (unsigned int d = 0; d < spacedim; ++d)
-    for (const unsigned int v : this->vertex_indices())
-      for (unsigned int e = 0; e < structdim; ++e)
-        A[d][e] += vertices[v][d] * TransformR2UAffine<structdim>::KA[v][e];
 
-  // b = vertex * Kb
-  Tensor<1, spacedim> b = point;
-  for (const unsigned int v : this->vertex_indices())
-    b -= vertices[v] * TransformR2UAffine<structdim>::Kb[v];
-
-  DerivativeForm<1, spacedim, structdim> A_inv = A.covariant_form().transpose();
-  return Point<structdim>(apply_transformation(A_inv, b));
+  const auto A_b =
+    GridTools::affine_cell_approximation<structdim, spacedim>(vertices);
+  DerivativeForm<1, spacedim, structdim> A_inv =
+    A_b.first.covariant_form().transpose();
+  return Point<structdim>(apply_transformation(A_inv, point - A_b.second));
 }
 
 
+
 template <int structdim, int dim, int spacedim>
 Point<spacedim>
 TriaAccessor<structdim, dim, spacedim>::center(

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.