From 53f0593feef89951e1b05eac1f7c6cdd37f1f14d Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Sat, 3 Jun 2017 16:39:14 +0200 Subject: [PATCH] Introduce function real_to_unit_cell_affine_approximation. --- include/deal.II/grid/tria_accessor.h | 25 +++++ source/grid/tria_accessor.cc | 133 +++++++++++++++++++++++++++ 2 files changed, 158 insertions(+) diff --git a/include/deal.II/grid/tria_accessor.h b/include/deal.II/grid/tria_accessor.h index d945e0a005..6bf44eacbc 100644 --- a/include/deal.II/grid/tria_accessor.h +++ b/include/deal.II/grid/tria_accessor.h @@ -1295,6 +1295,31 @@ public: */ Point intermediate_point(const Point &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 + real_to_unit_cell_affine_approximation (const Point &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$ diff --git a/source/grid/tria_accessor.cc b/source/grid/tria_accessor.cc index f569bbb85e..207fc909e1 100644 --- a/source/grid/tria_accessor.cc +++ b/source/grid/tria_accessor.cc @@ -1223,6 +1223,139 @@ TriaAccessor::intermediate_point (const Point + *
  • 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. + *
  • Solve the equation A x_hat + b = p for x_hat + * + * + * 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 + struct TransformR2UAffine + { + static const double KA[GeometryInfo::vertices_per_cell][dim]; + static const double Kb[GeometryInfo::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 +Point +TriaAccessor +::real_to_unit_cell_affine_approximation (const Point &point) const +{ + // A = vertex * KA + DerivativeForm<1,structdim,spacedim> A; + + // copy vertices to avoid expensive resolution of vertex index inside loop + std::array, GeometryInfo::vertices_per_cell> vertices; + for (unsigned int v=0; v::vertices_per_cell; ++v) + vertices[v] = this->vertex(v); + for (unsigned int d=0; d::vertices_per_cell; ++v) + for (unsigned int e=0; e::KA[v][e]; + + // b = vertex * Kb + Tensor<1,spacedim> b = point; + for (unsigned int v=0; v::vertices_per_cell; ++v) + b -= this->vertex(v) * TransformR2UAffine::Kb[v]; + + DerivativeForm<1,spacedim,structdim> A_inv = A.covariant_form().transpose(); + return Point(apply_transformation(A_inv, b)); +} + + template Point TriaAccessor::center (const bool respect_manifold, -- 2.39.5