From: Martin Kronbichler Date: Fri, 2 Oct 2020 10:14:26 +0000 (+0200) Subject: Allow to compute affine approximation on arbitrary set of points X-Git-Tag: v9.3.0-rc1~1038^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8d1b3a2d5a941501ce5c9a2cb498e45553972822;p=dealii.git Allow to compute affine approximation on arbitrary set of points --- diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index d781ceba28..f7a07d15ac 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -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 A as the first argument and the vector b 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 + std::pair, Tensor<1, spacedim>> + affine_cell_approximation(const ArrayView> &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 diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 62212f8881..2499c4f272 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -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 + 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}; + } // namespace + + + + template + std::pair, Tensor<1, spacedim>> + affine_cell_approximation(const ArrayView> &vertices) + { + AssertDimension(vertices.size(), GeometryInfo::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::vertices_per_cell; ++v) + for (unsigned int e = 0; e < dim; ++e) + A[d][e] += vertices[v][d] * TransformR2UAffine::KA[v][e]; + + // b = vertex * Kb + Tensor<1, spacedim> b; + for (unsigned int v = 0; v < GeometryInfo::vertices_per_cell; ++v) + b += vertices[v] * TransformR2UAffine::Kb[v]; + + return std::make_pair(A, b); + } + + + template Vector compute_aspect_ratio_of_cells(const Mapping & mapping, diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in index fb09c045f6..713d8c77bd 100644 --- a/source/grid/grid_tools.inst.in +++ b/source/grid/grid_tools.inst.in @@ -204,6 +204,12 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) volume(const Triangulation &, const Mapping &); + template std::pair< + DerivativeForm<1, deal_II_dimension, deal_II_space_dimension>, + Tensor<1, deal_II_space_dimension>> + affine_cell_approximation( + const ArrayView> &); + template BoundingBox compute_bounding_box( const Triangulation &); diff --git a/source/grid/tria_accessor.cc b/source/grid/tria_accessor.cc index 80ac4c3c49..002031c0b3 100644 --- a/source/grid/tria_accessor.cc +++ b/source/grid/tria_accessor.cc @@ -1700,131 +1700,26 @@ TriaAccessor::intermediate_point( } -namespace -{ - /** - * The algorithm to compute the affine approximation to the point on the - * unit cell does the following steps: - *
    - *
  • 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}; -} // namespace - 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 (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::KA[v][e]; - // b = vertex * Kb - Tensor<1, spacedim> b = point; - for (const unsigned int v : this->vertex_indices()) - b -= vertices[v] * TransformR2UAffine::Kb[v]; - - DerivativeForm<1, spacedim, structdim> A_inv = A.covariant_form().transpose(); - return Point(apply_transformation(A_inv, b)); + const auto A_b = + GridTools::affine_cell_approximation(vertices); + DerivativeForm<1, spacedim, structdim> A_inv = + A_b.first.covariant_form().transpose(); + return Point(apply_transformation(A_inv, point - A_b.second)); } + template Point TriaAccessor::center(