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
+ 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,
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> &);
}
-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(
projects / dealii.git / commitdiff