]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Introduce function real_to_unit_cell_affine_approximation.
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Sat, 3 Jun 2017 14:39:14 +0000 (16:39 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Mon, 5 Jun 2017 14:08:48 +0000 (16:08 +0200)
include/deal.II/grid/tria_accessor.h
source/grid/tria_accessor.cc

index d945e0a00574d599b836a92d783f409057d29ab9..6bf44eacbccd2ad911ad2992f97d05fce11f19eb 100644 (file)
@@ -1295,6 +1295,31 @@ public:
    */
   Point<spacedim> intermediate_point(const Point<structdim> &coordinates) const;
 
+  /**
+   * This function computes a fast approximate transformation from the real to
+   * the unit cell by inversion of an affine approximation of the $d$-linear
+   * function from the reference $d$-dimensional cell.
+   *
+   * The affine approximation of the unit to real cell mapping is found by a
+   * least squares fit of an affine function to the $2^d$ vertices of the
+   * present object. For any valid mesh cell whose geometry is not degenerate,
+   * this operation results in a unique affine mapping. Thus, this function
+   * will return a finite result for all given input points, even in cases
+   * where the actual transformation by an actual bi-/trilinear or higher
+   * order mapping might be singular. Besides only approximating the mapping
+   * from the vertex points, this function also ignores the attached manifold
+   * descriptions. The result is only 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.
+   *
+   * For exact transformations to the unit cell, use
+   * Mapping::transform_real_to_unit_cell().
+   *
+   * @note If dim<spacedim we first project p onto the plane.
+   */
+  Point<structdim>
+  real_to_unit_cell_affine_approximation (const Point<spacedim> &point) const;
+
   /**
    * Center of the object. The center of an object is defined to be the
    * average of the locations of the vertices, which is also where a $Q_1$
index f569bbb85e566786bfcc419f42c6c28ee733a36f..207fc909e15fb668f9a248c5ece43977e3ec4dfb 100644 (file)
@@ -1223,6 +1223,139 @@ TriaAccessor<structdim, dim, spacedim>::intermediate_point (const Point<structdi
 }
 
 
+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};
+}
+
+
+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 (unsigned int v=0; v<GeometryInfo<structdim>::vertices_per_cell; ++v)
+    vertices[v] = this->vertex(v);
+  for (unsigned int d=0; d<spacedim; ++d)
+    for (unsigned int v=0; v<GeometryInfo<structdim>::vertices_per_cell; ++v)
+      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 (unsigned int v=0; v<GeometryInfo<structdim>::vertices_per_cell; ++v)
+    b -= this->vertex(v) * TransformR2UAffine<structdim>::Kb[v];
+
+  DerivativeForm<1,spacedim,structdim> A_inv = A.covariant_form().transpose();
+  return Point<structdim>(apply_transformation(A_inv, b));
+}
+
+
 template <int structdim, int dim, int spacedim>
 Point<spacedim>
 TriaAccessor<structdim, dim, spacedim>::center (const bool respect_manifold,

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.