From: David Wells Date: Sun, 15 Oct 2023 22:30:58 +0000 (-0400) Subject: Move geometric GridTools functions to a new file. X-Git-Tag: relicensing~389^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F16149%2Fhead;p=dealii.git Move geometric GridTools functions to a new file. grid_tools.cc is presently our most expensive file to compile. This is the first commit in a sequence which will chop it up in the same way we chopped up VectorTools. --- diff --git a/include/deal.II/fe/mapping_q_internal.h b/include/deal.II/fe/mapping_q_internal.h index 80250bbdd7..dd56d422f0 100644 --- a/include/deal.II/fe/mapping_q_internal.h +++ b/include/deal.II/fe/mapping_q_internal.h @@ -31,7 +31,7 @@ #include #include -#include +#include #include #include diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index 7faf8c8afc..b42ec6043f 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -32,6 +32,7 @@ #include #include +#include #include #include #include @@ -132,286 +133,6 @@ namespace internal */ namespace GridTools { - /** - * @name Information about meshes and cells - */ - /** @{ */ - - /** - * Return the diameter of a triangulation. The diameter is computed using - * only the vertices, i.e. if the diameter should be larger than the maximal - * distance between boundary vertices due to a higher order mapping, then - * this function will not catch this. - */ - template - double - diameter(const Triangulation &tria); - - /** - * Compute the volume (i.e. the dim-dimensional measure) of the - * triangulation. We compute the measure using the integral $\sum_K \int_K 1 - * \; dx$ where $K$ are the cells of the given triangulation. The integral - * is approximated via quadrature. This version of the function uses a - * linear mapping to compute the JxW values on each cell. - * - * If the triangulation is a dim-dimensional one embedded in a higher - * dimensional space of dimension spacedim, then the value returned is the - * dim-dimensional measure. For example, for a two-dimensional triangulation - * in three-dimensional space, the value returned is the area of the surface - * so described. (This obviously makes sense since the spacedim-dimensional - * measure of a dim-dimensional triangulation would always be zero if dim @< - * spacedim). - * - * This function also works for objects of type - * parallel::distributed::Triangulation, in which case the function is a - * collective operation. - * - * @param tria The triangulation. - * @return The dim-dimensional measure of the domain described by the - * triangulation, as discussed above. - */ - template - double - volume(const Triangulation &tria); - - /** - * Compute the volume (i.e. the dim-dimensional measure) of the - * triangulation. We compute the measure using the integral $\sum_K \int_K 1 - * \; dx$ where $K$ are the cells of the given triangulation. The integral - * is approximated via quadrature for which we use the mapping argument. - * - * If the triangulation is a dim-dimensional one embedded in a higher - * dimensional space of dimension spacedim, then the value returned is the - * dim-dimensional measure. For example, for a two-dimensional triangulation - * in three-dimensional space, the value returned is the area of the surface - * so described. (This obviously makes sense since the spacedim-dimensional - * measure of a dim-dimensional triangulation would always be zero if dim @< - * spacedim. - * - * This function also works for objects of type - * parallel::distributed::Triangulation, in which case the function is a - * collective operation. - * - * @param tria The triangulation. - * @param mapping The Mapping which computes the Jacobians used to - * approximate the volume via quadrature. Explicitly using a higher-order - * Mapping (i.e., instead of using the other version of this function) will - * result in a more accurate approximation of the volume on Triangulations - * with curvature described by Manifold objects. - * @return The dim-dimensional measure of the domain described by the - * triangulation, as discussed above. - */ - template - double - volume(const Triangulation &tria, - const Mapping &mapping); - - /** - * Return an approximation of the diameter of the smallest active cell of a - * triangulation. See step-24 for an example of use of this function. - * - * Notice that, even if you pass a non-trivial mapping, the returned value is - * computed only using information on the vertices of the triangulation, - * possibly transformed by the mapping. While this is accurate most of the - * times, it may fail to give the correct result when the triangulation - * contains very distorted cells. - */ - template - double - minimal_cell_diameter( - const Triangulation &triangulation, - const Mapping &mapping = - (ReferenceCells::get_hypercube() -#ifndef _MSC_VER - .template get_default_linear_mapping() -#else - .ReferenceCell::get_default_linear_mapping() -#endif - )); - - /** - * Return an approximation of the diameter of the largest active cell of a - * triangulation. - * - * Notice that, even if you pass a non-trivial mapping to this function, the - * returned value is computed only using information on the vertices of the - * triangulation, possibly transformed by the mapping. While this is accurate - * most of the times, it may fail to give the correct result when the - * triangulation contains very distorted cells. - */ - template - double - maximal_cell_diameter( - const Triangulation &triangulation, - const Mapping &mapping = - (ReferenceCells::get_hypercube() -#ifndef _MSC_VER - .template get_default_linear_mapping() -#else - .ReferenceCell::get_default_linear_mapping() -#endif - )); - - /** - * Given a list of vertices (typically obtained using - * Triangulation::get_vertices()) as the first, and a list of vertex indices - * that characterize a single cell as the second argument, return the - * measure (area, volume) of this cell. If this is a real cell, then you can - * get the same result using cell-@>measure(), but this - * function also works for cells that do not exist except that you make it - * up by naming its vertices from the list. - * - * The size of @p vertex_indices, combined with `dim`, implicitly encodes - * the ReferenceCell type of the provided cell. For example, if `dim == 2` and - * `vertex_indices.size() == 3` then the cell is a triangle, but if - * `dim == 2` and `vertex_indices.size() == 4` then the cell is a - * quadrilateral. A std::vector is implicitly convertible to an ArrayView, so - * it can be passed directly to this function. See the ArrayView class for - * more information. - * - * @note This function is only implemented for codimension zero objects. - */ - template - double - cell_measure(const std::vector> &all_vertices, - const ArrayView &vertex_indices); - - /** - * 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 - * @p mapping. The size of the vector that is returned equals the number of - * active cells. The vector contains zero for non locally-owned cells. The - * aspect ratio of a cell is defined as the ratio of the maximum to minimum - * singular value of the Jacobian, taking the maximum over all quadrature - * points of a quadrature rule specified via @p quadrature. For example, for - * the special case of rectangular elements in 2d with dimensions $a$ and $b$ - * ($a \geq b$), this function returns the usual aspect ratio definition - * $a/b$. The above definition using singular values is a generalization to - * arbitrarily deformed elements. This function is intended to be used for - * $d=2,3$ space dimensions, but it can also be used for $d=1$ returning a - * value of 1. - * - * @note Inverted elements do not throw an exception. Instead, a value of inf - * is written into the vector in case of inverted elements. - * - * @note Make sure to use enough quadrature points for a precise calculation - * of the aspect ratio in case of deformed elements. - * - * @note In parallel computations the return value will have the length - * n_active_cells but the aspect ratio is only computed for the cells that - * are locally owned and placed at index CellAccessor::active_cell_index(), - * respectively. All other values are set to 0. - * - * @note This function can only be used if deal.II was configured with - * support for LAPACK. - */ - template - Vector - compute_aspect_ratio_of_cells(const Mapping &mapping, - const Triangulation &triangulation, - const Quadrature &quadrature); - - /** - * Computes the maximum aspect ratio by taking the maximum over all cells. - * - * @note When running in parallel with a Triangulation that supports MPI, - * this is a collective call and the return value is the maximum over all - * processors. - */ - template - double - compute_maximum_aspect_ratio(const Mapping &mapping, - const Triangulation &triangulation, - const Quadrature &quadrature); - - /** - * Compute the smallest box containing the entire triangulation. - * - * If the input triangulation is a `parallel::distributed::Triangulation`, - * then each processor will compute a bounding box enclosing all locally - * owned, ghost, and artificial cells. In the case of a domain without curved - * boundaries, these bounding boxes will all agree between processors because - * the union of the areas occupied by artificial and ghost cells equals the - * union of the areas occupied by the cells that other processors own. - * However, if the domain has curved boundaries, this is no longer the case. - * The bounding box returned may be appropriate for the current processor, - * but different from the bounding boxes computed on other processors. - */ - template - BoundingBox - compute_bounding_box(const Triangulation &triangulation); - - /** - * Return the point on the geometrical object @p object closest to the given - * point @p trial_point. For example, if @p object is a one-dimensional line - * or edge, then the returned point will be a point on the geodesic that - * connects the vertices as the manifold associated with the object sees it - * (i.e., the geometric line may be curved if it lives in a higher - * dimensional space). If the iterator points to a quadrilateral in a higher - * dimensional space, then the returned point lies within the convex hull of - * the vertices of the quad as seen by the associated manifold. - * - * @note This projection is usually not well-posed since there may be - * multiple points on the object that minimize the distance. The algorithm - * used in this function is robust (and the output is guaranteed to be on - * the given @p object) but may only provide a few correct digits if the - * object has high curvature. If your manifold supports it then the - * specialized function Manifold::project_to_manifold() may perform better. - */ - template - Point - project_to_object( - const Iterator &object, - const Point &trial_point); - - /** - * Return the arrays that define the coarse mesh of a Triangulation. This - * function is the inverse of Triangulation::create_triangulation() in the - * sense that if one called this function on a triangulation, then that - * triangulation could be recreated by some kind of refinement from the - * results of this function. - * - * The return value is a tuple with the vector of vertices, the vector of - * cells, and a SubCellData structure. The latter contains additional - * information about faces and lines. These three objects are exactly - * the arguments to Triangulation::create_triangulation(). - * - * This function is useful in cases where one needs to deconstruct a - * Triangulation or manipulate the numbering of the vertices in some way: an - * example is GridGenerator::merge_triangulations(). - */ - template - std:: - tuple>, std::vector>, SubCellData> - get_coarse_mesh_description(const Triangulation &tria); - - /** @} */ /** * @name Functions supporting the creation of meshes */ @@ -536,6 +257,26 @@ namespace GridTools * @name Rotating, stretching and otherwise transforming meshes */ /** @{ */ + /** + * Return the arrays that define the coarse mesh of a Triangulation. This + * function is the inverse of Triangulation::create_triangulation() in the + * sense that if one called this function on a triangulation, then that + * triangulation could be recreated by some kind of refinement from the + * results of this function. + * + * The return value is a tuple with the vector of vertices, the vector of + * cells, and a SubCellData structure. The latter contains additional + * information about faces and lines. These three objects are exactly + * the arguments to Triangulation::create_triangulation(). + * + * This function is useful in cases where one needs to deconstruct a + * Triangulation or manipulate the numbering of the vertices in some way: an + * example is GridGenerator::merge_triangulations(). + */ + template + std:: + tuple>, std::vector>, SubCellData> + get_coarse_mesh_description(const Triangulation &tria); /** * Transform the vertices of the given triangulation by applying the @@ -2130,34 +1871,6 @@ namespace GridTools const MeshType &mesh, const double layer_thickness); - /** - * Compute and return a bounding box, defined through a pair of points - * bottom left and top right, that surrounds a subdomain of the @p mesh. - * Here, the "subdomain" consists of exactly all of those - * active cells for which the @p predicate returns @p true. - * - * For a description of how @p predicate works, - * see compute_active_cell_halo_layer(). - * - * @note This function was written before the BoundingBox class was invented. - * Consequently, it returns a pair of points, rather than a BoundingBox - * object as one may expect. However, BoundingBox has a conversion constructor - * from pairs of points, so the result of this function can still be assigned - * to a BoundingBox object. - * - * @dealiiConceptRequires{concepts::is_triangulation_or_dof_handler} - */ - template - DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler) - std::pair< - Point, - Point> compute_bounding_box(const MeshType &mesh, - const std::function - &predicate); - /** * Compute a collection of bounding boxes so that all active cells for which * the given predicate is true, are completely enclosed in at least one of the @@ -4054,549 +3767,6 @@ namespace GridTools - namespace internal - { - namespace ProjectToObject - { - /** - * The method GridTools::project_to_object requires taking derivatives - * along the surface of a simplex. In general these cannot be - * approximated with finite differences but special differences of the - * form - * - * df/dx_i - df/dx_j - * - * can be approximated. This struct just stores - * the two derivatives approximated by the stencil (in the case of the - * example above i and j). - */ - struct CrossDerivative - { - const unsigned int direction_0; - const unsigned int direction_1; - - CrossDerivative(const unsigned int d0, const unsigned int d1); - }; - - inline CrossDerivative::CrossDerivative(const unsigned int d0, - const unsigned int d1) - : direction_0(d0) - , direction_1(d1) - {} - - - - /** - * Standard second-order approximation to the first derivative with a - * two-point centered scheme. This is used below in a 1d Newton method. - */ - template - inline auto - centered_first_difference(const double center, - const double step, - const F &f) -> decltype(f(center) - f(center)) - { - return (f(center + step) - f(center - step)) / (2.0 * step); - } - - - - /** - * Standard second-order approximation to the second derivative with a - * three-point centered scheme. This is used below in a 1d Newton method. - */ - template - inline auto - centered_second_difference(const double center, - const double step, - const F &f) -> decltype(f(center) - f(center)) - { - return (f(center + step) - 2.0 * f(center) + f(center - step)) / - (step * step); - } - - - - /** - * Fourth order approximation of the derivative - * - * df/dx_i - df/dx_j - * - * where i and j are specified by @p - * cross_derivative. The derivative approximation is at @p center with a - * step size of @p step and function @p f. - */ - template - inline auto - cross_stencil( - const CrossDerivative cross_derivative, - const Tensor<1, GeometryInfo::vertices_per_cell> ¢er, - const double step, - const F &f) -> decltype(f(center) - f(center)) - { - Tensor<1, GeometryInfo::vertices_per_cell> simplex_vector; - simplex_vector[cross_derivative.direction_0] = 0.5 * step; - simplex_vector[cross_derivative.direction_1] = -0.5 * step; - return (-4.0 * f(center) - 1.0 * f(center + simplex_vector) - - 1.0 / 3.0 * f(center - simplex_vector) + - 16.0 / 3.0 * f(center + 0.5 * simplex_vector)) / - step; - } - - - - /** - * The optimization algorithm used in GridTools::project_to_object is - * essentially a gradient descent method. This function computes entries - * in the gradient of the objective function; see the description in the - * comments inside GridTools::project_to_object for more information. - */ - template - inline double - gradient_entry( - const unsigned int row_n, - const unsigned int dependent_direction, - const Point &p0, - const Tensor<1, GeometryInfo::vertices_per_cell> ¢er, - const double step, - const F &f) - { - Assert(row_n < GeometryInfo::vertices_per_cell && - dependent_direction < - GeometryInfo::vertices_per_cell, - ExcMessage("This function assumes that the last weight is a " - "dependent variable (and hence we cannot take its " - "derivative directly).")); - Assert(row_n != dependent_direction, - ExcMessage( - "We cannot differentiate with respect to the variable " - "that is assumed to be dependent.")); - - const Point manifold_point = f(center); - const Tensor<1, spacedim> stencil_value = cross_stencil( - {row_n, dependent_direction}, center, step, f); - double entry = 0.0; - for (unsigned int dim_n = 0; dim_n < spacedim; ++dim_n) - entry += - -2.0 * (p0[dim_n] - manifold_point[dim_n]) * stencil_value[dim_n]; - return entry; - } - - /** - * Project onto a d-linear object. This is more accurate than the - * general algorithm in project_to_object but only works for geometries - * described by linear, bilinear, or trilinear mappings. - */ - template - Point - project_to_d_linear_object(const Iterator &object, - const Point &trial_point) - { - // let's look at this for simplicity for a quadrilateral - // (structdim==2) in a space with spacedim>2 (notate trial_point by - // y): all points on the surface are given by - // x(\xi) = sum_i v_i phi_x(\xi) - // where v_i are the vertices of the quadrilateral, and - // \xi=(\xi_1,\xi_2) are the reference coordinates of the - // quadrilateral. so what we are trying to do is find a point x on the - // surface that is closest to the point y. there are different ways to - // solve this problem, but in the end it's a nonlinear problem and we - // have to find reference coordinates \xi so that J(\xi) = 1/2 || - // x(\xi)-y ||^2 is minimal. x(\xi) is a function that is - // structdim-linear in \xi, so J(\xi) is a polynomial of degree - // 2*structdim that we'd like to minimize. unless structdim==1, we'll - // have to use a Newton method to find the answer. This leads to the - // following formulation of Newton steps: - // - // Given \xi_k, find \delta\xi_k so that - // H_k \delta\xi_k = - F_k - // where H_k is an approximation to the second derivatives of J at - // \xi_k, and F_k is the first derivative of J. We'll iterate this a - // number of times until the right hand side is small enough. As a - // stopping criterion, we terminate if ||\delta\xi|| xi; - for (unsigned int d = 0; d < structdim; ++d) - xi[d] = 0.5; - - Point x_k; - for (const unsigned int i : GeometryInfo::vertex_indices()) - x_k += object->vertex(i) * - GeometryInfo::d_linear_shape_function(xi, i); - - do - { - Tensor<1, structdim> F_k; - for (const unsigned int i : - GeometryInfo::vertex_indices()) - F_k += - (x_k - trial_point) * object->vertex(i) * - GeometryInfo::d_linear_shape_function_gradient(xi, - i); - - Tensor<2, structdim> H_k; - for (const unsigned int i : - GeometryInfo::vertex_indices()) - for (const unsigned int j : - GeometryInfo::vertex_indices()) - { - Tensor<2, structdim> tmp = outer_product( - GeometryInfo::d_linear_shape_function_gradient( - xi, i), - GeometryInfo::d_linear_shape_function_gradient( - xi, j)); - H_k += (object->vertex(i) * object->vertex(j)) * tmp; - } - - const Tensor<1, structdim> delta_xi = -invert(H_k) * F_k; - xi += delta_xi; - - x_k = Point(); - for (const unsigned int i : - GeometryInfo::vertex_indices()) - x_k += object->vertex(i) * - GeometryInfo::d_linear_shape_function(xi, i); - - if (delta_xi.norm() < 1e-7) - break; - } - while (true); - - return x_k; - } - } // namespace ProjectToObject - } // namespace internal - - - - namespace internal - { - // We hit an internal compiler error in ICC 15 if we define this as a lambda - // inside the project_to_object function below. - template - inline bool - weights_are_ok( - const Tensor<1, GeometryInfo::vertices_per_cell> &v) - { - // clang has trouble figuring out structdim here, so define it - // again: - static const std::size_t n_vertices_per_cell = - Tensor<1, GeometryInfo::vertices_per_cell>:: - n_independent_components; - std::array copied_weights; - for (unsigned int i = 0; i < n_vertices_per_cell; ++i) - { - copied_weights[i] = v[i]; - if (v[i] < 0.0 || v[i] > 1.0) - return false; - } - - // check the sum: try to avoid some roundoff errors by summing in order - std::sort(copied_weights.begin(), copied_weights.end()); - const double sum = - std::accumulate(copied_weights.begin(), copied_weights.end(), 0.0); - return std::abs(sum - 1.0) < 1e-10; // same tolerance used in manifold.cc - } - } // namespace internal - - template - Point - project_to_object( - const Iterator &object, - const Point &trial_point) - { - const int spacedim = Iterator::AccessorType::space_dimension; - const int structdim = Iterator::AccessorType::structure_dimension; - - Point projected_point = trial_point; - - if (structdim >= spacedim) - return projected_point; - else if (structdim == 1 || structdim == 2) - { - using namespace internal::ProjectToObject; - // Try to use the special flat algorithm for quads (this is better - // than the general algorithm in 3d). This does not take into account - // whether projected_point is outside the quad, but we optimize along - // lines below anyway: - const int dim = Iterator::AccessorType::dimension; - const Manifold &manifold = object->get_manifold(); - if (structdim == 2 && dynamic_cast *>( - &manifold) != nullptr) - { - projected_point = - project_to_d_linear_object( - object, trial_point); - } - else - { - // We want to find a point on the convex hull (defined by the - // vertices of the object and the manifold description) that is - // relatively close to the trial point. This has a few issues: - // - // 1. For a general convex hull we are not guaranteed that a unique - // minimum exists. - // 2. The independent variables in the optimization process are the - // weights given to Manifold::get_new_point, which must sum to 1, - // so we cannot use standard finite differences to approximate a - // gradient. - // - // There is not much we can do about 1., but for 2. we can derive - // finite difference stencils that work on a structdim-dimensional - // simplex and rewrite the optimization problem to use those - // instead. Consider the structdim 2 case and let - // - // F(c0, c1, c2, c3) = Manifold::get_new_point(vertices, {c0, c1, - // c2, c3}) - // - // where {c0, c1, c2, c3} are the weights for the four vertices on - // the quadrilateral. We seek to minimize the Euclidean distance - // between F(...) and trial_point. We can solve for c3 in terms of - // the other weights and get, for one coordinate direction - // - // d/dc0 ((x0 - F(c0, c1, c2, 1 - c0 - c1 - c2))^2) - // = -2(x0 - F(...)) (d/dc0 F(...) - d/dc3 F(...)) - // - // where we substitute back in for c3 after taking the - // derivative. We can compute a stencil for the cross derivative - // d/dc0 - d/dc3: this is exactly what cross_stencil approximates - // (and gradient_entry computes the sum over the independent - // variables). Below, we somewhat arbitrarily pick the last - // component as the dependent one. - // - // Since we can now calculate derivatives of the objective - // function we can use gradient descent to minimize it. - // - // Of course, this is much simpler in the structdim = 1 case (we - // could rewrite the projection as a 1d optimization problem), but - // to reduce the potential for bugs we use the same code in both - // cases. - const double step_size = object->diameter() / 64.0; - - constexpr unsigned int n_vertices_per_cell = - GeometryInfo::vertices_per_cell; - - std::array, n_vertices_per_cell> vertices; - for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell; - ++vertex_n) - vertices[vertex_n] = object->vertex(vertex_n); - - auto get_point_from_weights = - [&](const Tensor<1, n_vertices_per_cell> &weights) - -> Point { - return object->get_manifold().get_new_point( - make_array_view(vertices.begin(), vertices.end()), - make_array_view(weights.begin_raw(), weights.end_raw())); - }; - - // pick the initial weights as (normalized) inverse distances from - // the trial point: - Tensor<1, n_vertices_per_cell> guess_weights; - double guess_weights_sum = 0.0; - for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell; - ++vertex_n) - { - const double distance = - vertices[vertex_n].distance(trial_point); - if (distance == 0.0) - { - guess_weights = 0.0; - guess_weights[vertex_n] = 1.0; - guess_weights_sum = 1.0; - break; - } - else - { - guess_weights[vertex_n] = 1.0 / distance; - guess_weights_sum += guess_weights[vertex_n]; - } - } - guess_weights /= guess_weights_sum; - Assert(internal::weights_are_ok(guess_weights), - ExcInternalError()); - - // The optimization algorithm consists of two parts: - // - // 1. An outer loop where we apply the gradient descent algorithm. - // 2. An inner loop where we do a line search to find the optimal - // length of the step one should take in the gradient direction. - // - for (unsigned int outer_n = 0; outer_n < 40; ++outer_n) - { - const unsigned int dependent_direction = - n_vertices_per_cell - 1; - Tensor<1, n_vertices_per_cell> current_gradient; - for (unsigned int row_n = 0; row_n < n_vertices_per_cell; - ++row_n) - { - if (row_n != dependent_direction) - { - current_gradient[row_n] = - gradient_entry( - row_n, - dependent_direction, - trial_point, - guess_weights, - step_size, - get_point_from_weights); - - current_gradient[dependent_direction] -= - current_gradient[row_n]; - } - } - - // We need to travel in the -gradient direction, as noted - // above, but we may not want to take a full step in that - // direction; instead, guess that we will go -0.5*gradient and - // do quasi-Newton iteration to pick the best multiplier. The - // goal is to find a scalar alpha such that - // - // F(x - alpha g) - // - // is minimized, where g is the gradient and F is the - // objective function. To find the optimal value we find roots - // of the derivative of the objective function with respect to - // alpha by Newton iteration, where we approximate the first - // and second derivatives of F(x - alpha g) with centered - // finite differences. - double gradient_weight = -0.5; - auto gradient_weight_objective_function = - [&](const double gradient_weight_guess) -> double { - return (trial_point - - get_point_from_weights(guess_weights + - gradient_weight_guess * - current_gradient)) - .norm_square(); - }; - - for (unsigned int inner_n = 0; inner_n < 10; ++inner_n) - { - const double update_numerator = centered_first_difference( - gradient_weight, - step_size, - gradient_weight_objective_function); - const double update_denominator = - centered_second_difference( - gradient_weight, - step_size, - gradient_weight_objective_function); - - // avoid division by zero. Note that we limit the gradient - // weight below - if (std::abs(update_denominator) == 0.0) - break; - gradient_weight = - gradient_weight - update_numerator / update_denominator; - - // Put a fairly lenient bound on the largest possible - // gradient (things tend to be locally flat, so the gradient - // itself is usually small) - if (std::abs(gradient_weight) > 10) - { - gradient_weight = -10.0; - break; - } - } - - // It only makes sense to take convex combinations with weights - // between zero and one. If the update takes us outside of this - // region then rescale the update to stay within the region and - // try again - Tensor<1, n_vertices_per_cell> tentative_weights = - guess_weights + gradient_weight * current_gradient; - - double new_gradient_weight = gradient_weight; - for (unsigned int iteration_count = 0; iteration_count < 40; - ++iteration_count) - { - if (internal::weights_are_ok(tentative_weights)) - break; - - for (unsigned int i = 0; i < n_vertices_per_cell; ++i) - { - if (tentative_weights[i] < 0.0) - { - tentative_weights -= - (tentative_weights[i] / current_gradient[i]) * - current_gradient; - } - if (tentative_weights[i] < 0.0 || - 1.0 < tentative_weights[i]) - { - new_gradient_weight /= 2.0; - tentative_weights = - guess_weights + - new_gradient_weight * current_gradient; - } - } - } - - // the update might still send us outside the valid region, so - // check again and quit if the update is still not valid - if (!internal::weights_are_ok(tentative_weights)) - break; - - // if we cannot get closer by traveling in the gradient - // direction then quit - if (get_point_from_weights(tentative_weights) - .distance(trial_point) < - get_point_from_weights(guess_weights).distance(trial_point)) - guess_weights = tentative_weights; - else - break; - Assert(internal::weights_are_ok(guess_weights), - ExcInternalError()); - } - Assert(internal::weights_are_ok(guess_weights), - ExcInternalError()); - projected_point = get_point_from_weights(guess_weights); - } - - // if structdim == 2 and the optimal point is not on the interior then - // we may be able to get a more accurate result by projecting onto the - // lines. - if (structdim == 2) - { - std::array, GeometryInfo::lines_per_cell> - line_projections; - for (unsigned int line_n = 0; - line_n < GeometryInfo::lines_per_cell; - ++line_n) - { - line_projections[line_n] = - project_to_object(object->line(line_n), trial_point); - } - std::sort(line_projections.begin(), - line_projections.end(), - [&](const Point &a, const Point &b) { - return a.distance(trial_point) < - b.distance(trial_point); - }); - if (line_projections[0].distance(trial_point) < - projected_point.distance(trial_point)) - projected_point = line_projections[0]; - } - } - else - { - Assert(false, ExcNotImplemented()); - return projected_point; - } - - return projected_point; - } - - - namespace internal { template + +#include +#include +#include +#include + +#include + +#include +#include + +#include +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +/** + * This namespace is a collection of algorithms working on triangulations, + * such as shifting or rotating triangulations, but also finding a cell that + * contains a given point. See the descriptions of the individual functions + * for more information. + * + * @ingroup grid + */ +namespace GridTools +{ + /** + * @name Information about meshes and cells + */ + /** @{ */ + + /** + * Return the diameter of a triangulation. The diameter is computed using + * only the vertices, i.e. if the diameter should be larger than the maximal + * distance between boundary vertices due to a higher order mapping, then + * this function will not catch this. + */ + template + double + diameter(const Triangulation &tria); + + /** + * Compute the volume (i.e. the dim-dimensional measure) of the + * triangulation. We compute the measure using the integral $\sum_K \int_K 1 + * \; dx$ where $K$ are the cells of the given triangulation. The integral + * is approximated via quadrature. This version of the function uses a + * linear mapping to compute the JxW values on each cell. + * + * If the triangulation is a dim-dimensional one embedded in a higher + * dimensional space of dimension spacedim, then the value returned is the + * dim-dimensional measure. For example, for a two-dimensional triangulation + * in three-dimensional space, the value returned is the area of the surface + * so described. (This obviously makes sense since the spacedim-dimensional + * measure of a dim-dimensional triangulation would always be zero if dim @< + * spacedim). + * + * This function also works for objects of type + * parallel::distributed::Triangulation, in which case the function is a + * collective operation. + * + * @param tria The triangulation. + * @return The dim-dimensional measure of the domain described by the + * triangulation, as discussed above. + */ + template + double + volume(const Triangulation &tria); + + /** + * Compute the volume (i.e. the dim-dimensional measure) of the + * triangulation. We compute the measure using the integral $\sum_K \int_K 1 + * \; dx$ where $K$ are the cells of the given triangulation. The integral + * is approximated via quadrature for which we use the mapping argument. + * + * If the triangulation is a dim-dimensional one embedded in a higher + * dimensional space of dimension spacedim, then the value returned is the + * dim-dimensional measure. For example, for a two-dimensional triangulation + * in three-dimensional space, the value returned is the area of the surface + * so described. (This obviously makes sense since the spacedim-dimensional + * measure of a dim-dimensional triangulation would always be zero if dim @< + * spacedim. + * + * This function also works for objects of type + * parallel::distributed::Triangulation, in which case the function is a + * collective operation. + * + * @param tria The triangulation. + * @param mapping The Mapping which computes the Jacobians used to + * approximate the volume via quadrature. Explicitly using a higher-order + * Mapping (i.e., instead of using the other version of this function) will + * result in a more accurate approximation of the volume on Triangulations + * with curvature described by Manifold objects. + * @return The dim-dimensional measure of the domain described by the + * triangulation, as discussed above. + */ + template + double + volume(const Triangulation &tria, + const Mapping &mapping); + + /** + * Return an approximation of the diameter of the smallest active cell of a + * triangulation. See step-24 for an example of use of this function. + * + * Notice that, even if you pass a non-trivial mapping, the returned value is + * computed only using information on the vertices of the triangulation, + * possibly transformed by the mapping. While this is accurate most of the + * times, it may fail to give the correct result when the triangulation + * contains very distorted cells. + */ + template + double + minimal_cell_diameter( + const Triangulation &triangulation, + const Mapping &mapping = + (ReferenceCells::get_hypercube() +#ifndef _MSC_VER + .template get_default_linear_mapping() +#else + .ReferenceCell::get_default_linear_mapping() +#endif + )); + + /** + * Return an approximation of the diameter of the largest active cell of a + * triangulation. + * + * Notice that, even if you pass a non-trivial mapping to this function, the + * returned value is computed only using information on the vertices of the + * triangulation, possibly transformed by the mapping. While this is accurate + * most of the times, it may fail to give the correct result when the + * triangulation contains very distorted cells. + */ + template + double + maximal_cell_diameter( + const Triangulation &triangulation, + const Mapping &mapping = + (ReferenceCells::get_hypercube() +#ifndef _MSC_VER + .template get_default_linear_mapping() +#else + .ReferenceCell::get_default_linear_mapping() +#endif + )); + + /** + * Given a list of vertices (typically obtained using + * Triangulation::get_vertices()) as the first, and a list of vertex indices + * that characterize a single cell as the second argument, return the + * measure (area, volume) of this cell. If this is a real cell, then you can + * get the same result using cell-@>measure(), but this + * function also works for cells that do not exist except that you make it + * up by naming its vertices from the list. + * + * The size of @p vertex_indices, combined with `dim`, implicitly encodes + * the ReferenceCell type of the provided cell. For example, if `dim == 2` and + * `vertex_indices.size() == 3` then the cell is a triangle, but if + * `dim == 2` and `vertex_indices.size() == 4` then the cell is a + * quadrilateral. A std::vector is implicitly convertible to an ArrayView, so + * it can be passed directly to this function. See the ArrayView class for + * more information. + * + * @note This function is only implemented for codimension zero objects. + */ + template + double + cell_measure(const std::vector> &all_vertices, + const ArrayView &vertex_indices); + + /** + * 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 + * @p mapping. The size of the vector that is returned equals the number of + * active cells. The vector contains zero for non locally-owned cells. The + * aspect ratio of a cell is defined as the ratio of the maximum to minimum + * singular value of the Jacobian, taking the maximum over all quadrature + * points of a quadrature rule specified via @p quadrature. For example, for + * the special case of rectangular elements in 2d with dimensions $a$ and $b$ + * ($a \geq b$), this function returns the usual aspect ratio definition + * $a/b$. The above definition using singular values is a generalization to + * arbitrarily deformed elements. This function is intended to be used for + * $d=2,3$ space dimensions, but it can also be used for $d=1$ returning a + * value of 1. + * + * @note Inverted elements do not throw an exception. Instead, a value of inf + * is written into the vector in case of inverted elements. + * + * @note Make sure to use enough quadrature points for a precise calculation + * of the aspect ratio in case of deformed elements. + * + * @note In parallel computations the return value will have the length + * n_active_cells but the aspect ratio is only computed for the cells that + * are locally owned and placed at index CellAccessor::active_cell_index(), + * respectively. All other values are set to 0. + * + * @note This function can only be used if deal.II was configured with + * support for LAPACK. + */ + template + Vector + compute_aspect_ratio_of_cells(const Mapping &mapping, + const Triangulation &triangulation, + const Quadrature &quadrature); + + /** + * Computes the maximum aspect ratio by taking the maximum over all cells. + * + * @note When running in parallel with a Triangulation that supports MPI, + * this is a collective call and the return value is the maximum over all + * processors. + */ + template + double + compute_maximum_aspect_ratio(const Mapping &mapping, + const Triangulation &triangulation, + const Quadrature &quadrature); + + /** + * Compute the smallest box containing the entire triangulation. + * + * If the input triangulation is a `parallel::distributed::Triangulation`, + * then each processor will compute a bounding box enclosing all locally + * owned, ghost, and artificial cells. In the case of a domain without curved + * boundaries, these bounding boxes will all agree between processors because + * the union of the areas occupied by artificial and ghost cells equals the + * union of the areas occupied by the cells that other processors own. + * However, if the domain has curved boundaries, this is no longer the case. + * The bounding box returned may be appropriate for the current processor, + * but different from the bounding boxes computed on other processors. + */ + template + BoundingBox + compute_bounding_box(const Triangulation &triangulation); + + /** + * Compute and return a bounding box, defined through a pair of points + * bottom left and top right, that surrounds a subdomain of the @p mesh. + * Here, the "subdomain" consists of exactly all of those + * active cells for which the @p predicate returns @p true. + * + * For a description of how @p predicate works, + * see compute_active_cell_halo_layer(). + * + * @note This function was written before the BoundingBox class was invented. + * Consequently, it returns a pair of points, rather than a BoundingBox + * object as one may expect. However, BoundingBox has a conversion constructor + * from pairs of points, so the result of this function can still be assigned + * to a BoundingBox object. + * + * @dealiiConceptRequires{concepts::is_triangulation_or_dof_handler} + */ + template + DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler) + std::pair< + Point, + Point> compute_bounding_box(const MeshType &mesh, + const std::function + &predicate); + + /** + * Return the point on the geometrical object @p object closest to the given + * point @p trial_point. For example, if @p object is a one-dimensional line + * or edge, then the returned point will be a point on the geodesic that + * connects the vertices as the manifold associated with the object sees it + * (i.e., the geometric line may be curved if it lives in a higher + * dimensional space). If the iterator points to a quadrilateral in a higher + * dimensional space, then the returned point lies within the convex hull of + * the vertices of the quad as seen by the associated manifold. + * + * @note This projection is usually not well-posed since there may be + * multiple points on the object that minimize the distance. The algorithm + * used in this function is robust (and the output is guaranteed to be on + * the given @p object) but may only provide a few correct digits if the + * object has high curvature. If your manifold supports it then the + * specialized function Manifold::project_to_manifold() may perform better. + */ + template + Point + project_to_object( + const Iterator &object, + const Point &trial_point); + /** @} */ +} // namespace GridTools + +#ifndef DOXYGEN +namespace GridTools +{ + namespace internal + { + namespace ProjectToObject + { + /** + * The method GridTools::project_to_object requires taking derivatives + * along the surface of a simplex. In general these cannot be + * approximated with finite differences but special differences of the + * form + * + * df/dx_i - df/dx_j + * + * can be approximated. This struct just stores + * the two derivatives approximated by the stencil (in the case of the + * example above i and j). + */ + struct CrossDerivative + { + const unsigned int direction_0; + const unsigned int direction_1; + + CrossDerivative(const unsigned int d0, const unsigned int d1); + }; + + inline CrossDerivative::CrossDerivative(const unsigned int d0, + const unsigned int d1) + : direction_0(d0) + , direction_1(d1) + {} + + + + /** + * Standard second-order approximation to the first derivative with a + * two-point centered scheme. This is used below in a 1d Newton method. + */ + template + inline auto + centered_first_difference(const double center, + const double step, + const F &f) -> decltype(f(center) - f(center)) + { + return (f(center + step) - f(center - step)) / (2.0 * step); + } + + + + /** + * Standard second-order approximation to the second derivative with a + * three-point centered scheme. This is used below in a 1d Newton method. + */ + template + inline auto + centered_second_difference(const double center, + const double step, + const F &f) -> decltype(f(center) - f(center)) + { + return (f(center + step) - 2.0 * f(center) + f(center - step)) / + (step * step); + } + + + + /** + * Fourth order approximation of the derivative + * + * df/dx_i - df/dx_j + * + * where i and j are specified by @p + * cross_derivative. The derivative approximation is at @p center with a + * step size of @p step and function @p f. + */ + template + inline auto + cross_stencil( + const CrossDerivative cross_derivative, + const Tensor<1, GeometryInfo::vertices_per_cell> ¢er, + const double step, + const F &f) -> decltype(f(center) - f(center)) + { + Tensor<1, GeometryInfo::vertices_per_cell> simplex_vector; + simplex_vector[cross_derivative.direction_0] = 0.5 * step; + simplex_vector[cross_derivative.direction_1] = -0.5 * step; + return (-4.0 * f(center) - 1.0 * f(center + simplex_vector) - + 1.0 / 3.0 * f(center - simplex_vector) + + 16.0 / 3.0 * f(center + 0.5 * simplex_vector)) / + step; + } + + + + /** + * The optimization algorithm used in GridTools::project_to_object is + * essentially a gradient descent method. This function computes entries + * in the gradient of the objective function; see the description in the + * comments inside GridTools::project_to_object for more information. + */ + template + inline double + gradient_entry( + const unsigned int row_n, + const unsigned int dependent_direction, + const Point &p0, + const Tensor<1, GeometryInfo::vertices_per_cell> ¢er, + const double step, + const F &f) + { + Assert(row_n < GeometryInfo::vertices_per_cell && + dependent_direction < + GeometryInfo::vertices_per_cell, + ExcMessage("This function assumes that the last weight is a " + "dependent variable (and hence we cannot take its " + "derivative directly).")); + Assert(row_n != dependent_direction, + ExcMessage( + "We cannot differentiate with respect to the variable " + "that is assumed to be dependent.")); + + const Point manifold_point = f(center); + const Tensor<1, spacedim> stencil_value = cross_stencil( + {row_n, dependent_direction}, center, step, f); + double entry = 0.0; + for (unsigned int dim_n = 0; dim_n < spacedim; ++dim_n) + entry += + -2.0 * (p0[dim_n] - manifold_point[dim_n]) * stencil_value[dim_n]; + return entry; + } + + /** + * Project onto a d-linear object. This is more accurate than the + * general algorithm in project_to_object but only works for geometries + * described by linear, bilinear, or trilinear mappings. + */ + template + Point + project_to_d_linear_object(const Iterator &object, + const Point &trial_point) + { + // let's look at this for simplicity for a quadrilateral + // (structdim==2) in a space with spacedim>2 (notate trial_point by + // y): all points on the surface are given by + // x(\xi) = sum_i v_i phi_x(\xi) + // where v_i are the vertices of the quadrilateral, and + // \xi=(\xi_1,\xi_2) are the reference coordinates of the + // quadrilateral. so what we are trying to do is find a point x on the + // surface that is closest to the point y. there are different ways to + // solve this problem, but in the end it's a nonlinear problem and we + // have to find reference coordinates \xi so that J(\xi) = 1/2 || + // x(\xi)-y ||^2 is minimal. x(\xi) is a function that is + // structdim-linear in \xi, so J(\xi) is a polynomial of degree + // 2*structdim that we'd like to minimize. unless structdim==1, we'll + // have to use a Newton method to find the answer. This leads to the + // following formulation of Newton steps: + // + // Given \xi_k, find \delta\xi_k so that + // H_k \delta\xi_k = - F_k + // where H_k is an approximation to the second derivatives of J at + // \xi_k, and F_k is the first derivative of J. We'll iterate this a + // number of times until the right hand side is small enough. As a + // stopping criterion, we terminate if ||\delta\xi|| xi; + for (unsigned int d = 0; d < structdim; ++d) + xi[d] = 0.5; + + Point x_k; + for (const unsigned int i : GeometryInfo::vertex_indices()) + x_k += object->vertex(i) * + GeometryInfo::d_linear_shape_function(xi, i); + + do + { + Tensor<1, structdim> F_k; + for (const unsigned int i : + GeometryInfo::vertex_indices()) + F_k += + (x_k - trial_point) * object->vertex(i) * + GeometryInfo::d_linear_shape_function_gradient(xi, + i); + + Tensor<2, structdim> H_k; + for (const unsigned int i : + GeometryInfo::vertex_indices()) + for (const unsigned int j : + GeometryInfo::vertex_indices()) + { + Tensor<2, structdim> tmp = outer_product( + GeometryInfo::d_linear_shape_function_gradient( + xi, i), + GeometryInfo::d_linear_shape_function_gradient( + xi, j)); + H_k += (object->vertex(i) * object->vertex(j)) * tmp; + } + + const Tensor<1, structdim> delta_xi = -invert(H_k) * F_k; + xi += delta_xi; + + x_k = Point(); + for (const unsigned int i : + GeometryInfo::vertex_indices()) + x_k += object->vertex(i) * + GeometryInfo::d_linear_shape_function(xi, i); + + if (delta_xi.norm() < 1e-7) + break; + } + while (true); + + return x_k; + } + } // namespace ProjectToObject + + // We hit an internal compiler error in ICC 15 if we define this as a lambda + // inside the project_to_object function below. + template + inline bool + weights_are_ok( + const Tensor<1, GeometryInfo::vertices_per_cell> &v) + { + // clang has trouble figuring out structdim here, so define it + // again: + static const std::size_t n_vertices_per_cell = + Tensor<1, GeometryInfo::vertices_per_cell>:: + n_independent_components; + std::array copied_weights; + for (unsigned int i = 0; i < n_vertices_per_cell; ++i) + { + copied_weights[i] = v[i]; + if (v[i] < 0.0 || v[i] > 1.0) + return false; + } + + // check the sum: try to avoid some roundoff errors by summing in order + std::sort(copied_weights.begin(), copied_weights.end()); + const double sum = + std::accumulate(copied_weights.begin(), copied_weights.end(), 0.0); + return std::abs(sum - 1.0) < 1e-10; // same tolerance used in manifold.cc + } + } // namespace internal + + template + Point + project_to_object( + const Iterator &object, + const Point &trial_point) + { + const int spacedim = Iterator::AccessorType::space_dimension; + const int structdim = Iterator::AccessorType::structure_dimension; + + Point projected_point = trial_point; + + if (structdim >= spacedim) + return projected_point; + else if (structdim == 1 || structdim == 2) + { + using namespace internal::ProjectToObject; + // Try to use the special flat algorithm for quads (this is better + // than the general algorithm in 3d). This does not take into account + // whether projected_point is outside the quad, but we optimize along + // lines below anyway: + const int dim = Iterator::AccessorType::dimension; + const Manifold &manifold = object->get_manifold(); + if (structdim == 2 && dynamic_cast *>( + &manifold) != nullptr) + { + projected_point = + project_to_d_linear_object( + object, trial_point); + } + else + { + // We want to find a point on the convex hull (defined by the + // vertices of the object and the manifold description) that is + // relatively close to the trial point. This has a few issues: + // + // 1. For a general convex hull we are not guaranteed that a unique + // minimum exists. + // 2. The independent variables in the optimization process are the + // weights given to Manifold::get_new_point, which must sum to 1, + // so we cannot use standard finite differences to approximate a + // gradient. + // + // There is not much we can do about 1., but for 2. we can derive + // finite difference stencils that work on a structdim-dimensional + // simplex and rewrite the optimization problem to use those + // instead. Consider the structdim 2 case and let + // + // F(c0, c1, c2, c3) = Manifold::get_new_point(vertices, {c0, c1, + // c2, c3}) + // + // where {c0, c1, c2, c3} are the weights for the four vertices on + // the quadrilateral. We seek to minimize the Euclidean distance + // between F(...) and trial_point. We can solve for c3 in terms of + // the other weights and get, for one coordinate direction + // + // d/dc0 ((x0 - F(c0, c1, c2, 1 - c0 - c1 - c2))^2) + // = -2(x0 - F(...)) (d/dc0 F(...) - d/dc3 F(...)) + // + // where we substitute back in for c3 after taking the + // derivative. We can compute a stencil for the cross derivative + // d/dc0 - d/dc3: this is exactly what cross_stencil approximates + // (and gradient_entry computes the sum over the independent + // variables). Below, we somewhat arbitrarily pick the last + // component as the dependent one. + // + // Since we can now calculate derivatives of the objective + // function we can use gradient descent to minimize it. + // + // Of course, this is much simpler in the structdim = 1 case (we + // could rewrite the projection as a 1d optimization problem), but + // to reduce the potential for bugs we use the same code in both + // cases. + const double step_size = object->diameter() / 64.0; + + constexpr unsigned int n_vertices_per_cell = + GeometryInfo::vertices_per_cell; + + std::array, n_vertices_per_cell> vertices; + for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell; + ++vertex_n) + vertices[vertex_n] = object->vertex(vertex_n); + + auto get_point_from_weights = + [&](const Tensor<1, n_vertices_per_cell> &weights) + -> Point { + return object->get_manifold().get_new_point( + make_array_view(vertices.begin(), vertices.end()), + make_array_view(weights.begin_raw(), weights.end_raw())); + }; + + // pick the initial weights as (normalized) inverse distances from + // the trial point: + Tensor<1, n_vertices_per_cell> guess_weights; + double guess_weights_sum = 0.0; + for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell; + ++vertex_n) + { + const double distance = + vertices[vertex_n].distance(trial_point); + if (distance == 0.0) + { + guess_weights = 0.0; + guess_weights[vertex_n] = 1.0; + guess_weights_sum = 1.0; + break; + } + else + { + guess_weights[vertex_n] = 1.0 / distance; + guess_weights_sum += guess_weights[vertex_n]; + } + } + guess_weights /= guess_weights_sum; + Assert(internal::weights_are_ok(guess_weights), + ExcInternalError()); + + // The optimization algorithm consists of two parts: + // + // 1. An outer loop where we apply the gradient descent algorithm. + // 2. An inner loop where we do a line search to find the optimal + // length of the step one should take in the gradient direction. + // + for (unsigned int outer_n = 0; outer_n < 40; ++outer_n) + { + const unsigned int dependent_direction = + n_vertices_per_cell - 1; + Tensor<1, n_vertices_per_cell> current_gradient; + for (unsigned int row_n = 0; row_n < n_vertices_per_cell; + ++row_n) + { + if (row_n != dependent_direction) + { + current_gradient[row_n] = + gradient_entry( + row_n, + dependent_direction, + trial_point, + guess_weights, + step_size, + get_point_from_weights); + + current_gradient[dependent_direction] -= + current_gradient[row_n]; + } + } + + // We need to travel in the -gradient direction, as noted + // above, but we may not want to take a full step in that + // direction; instead, guess that we will go -0.5*gradient and + // do quasi-Newton iteration to pick the best multiplier. The + // goal is to find a scalar alpha such that + // + // F(x - alpha g) + // + // is minimized, where g is the gradient and F is the + // objective function. To find the optimal value we find roots + // of the derivative of the objective function with respect to + // alpha by Newton iteration, where we approximate the first + // and second derivatives of F(x - alpha g) with centered + // finite differences. + double gradient_weight = -0.5; + auto gradient_weight_objective_function = + [&](const double gradient_weight_guess) -> double { + return (trial_point - + get_point_from_weights(guess_weights + + gradient_weight_guess * + current_gradient)) + .norm_square(); + }; + + for (unsigned int inner_n = 0; inner_n < 10; ++inner_n) + { + const double update_numerator = centered_first_difference( + gradient_weight, + step_size, + gradient_weight_objective_function); + const double update_denominator = + centered_second_difference( + gradient_weight, + step_size, + gradient_weight_objective_function); + + // avoid division by zero. Note that we limit the gradient + // weight below + if (std::abs(update_denominator) == 0.0) + break; + gradient_weight = + gradient_weight - update_numerator / update_denominator; + + // Put a fairly lenient bound on the largest possible + // gradient (things tend to be locally flat, so the gradient + // itself is usually small) + if (std::abs(gradient_weight) > 10) + { + gradient_weight = -10.0; + break; + } + } + + // It only makes sense to take convex combinations with weights + // between zero and one. If the update takes us outside of this + // region then rescale the update to stay within the region and + // try again + Tensor<1, n_vertices_per_cell> tentative_weights = + guess_weights + gradient_weight * current_gradient; + + double new_gradient_weight = gradient_weight; + for (unsigned int iteration_count = 0; iteration_count < 40; + ++iteration_count) + { + if (internal::weights_are_ok(tentative_weights)) + break; + + for (unsigned int i = 0; i < n_vertices_per_cell; ++i) + { + if (tentative_weights[i] < 0.0) + { + tentative_weights -= + (tentative_weights[i] / current_gradient[i]) * + current_gradient; + } + if (tentative_weights[i] < 0.0 || + 1.0 < tentative_weights[i]) + { + new_gradient_weight /= 2.0; + tentative_weights = + guess_weights + + new_gradient_weight * current_gradient; + } + } + } + + // the update might still send us outside the valid region, so + // check again and quit if the update is still not valid + if (!internal::weights_are_ok(tentative_weights)) + break; + + // if we cannot get closer by traveling in the gradient + // direction then quit + if (get_point_from_weights(tentative_weights) + .distance(trial_point) < + get_point_from_weights(guess_weights).distance(trial_point)) + guess_weights = tentative_weights; + else + break; + Assert(internal::weights_are_ok(guess_weights), + ExcInternalError()); + } + Assert(internal::weights_are_ok(guess_weights), + ExcInternalError()); + projected_point = get_point_from_weights(guess_weights); + } + + // if structdim == 2 and the optimal point is not on the interior then + // we may be able to get a more accurate result by projecting onto the + // lines. + if (structdim == 2) + { + std::array, GeometryInfo::lines_per_cell> + line_projections; + for (unsigned int line_n = 0; + line_n < GeometryInfo::lines_per_cell; + ++line_n) + { + line_projections[line_n] = + project_to_object(object->line(line_n), trial_point); + } + std::sort(line_projections.begin(), + line_projections.end(), + [&](const Point &a, const Point &b) { + return a.distance(trial_point) < + b.distance(trial_point); + }); + if (line_projections[0].distance(trial_point) < + projected_point.distance(trial_point)) + projected_point = line_projections[0]; + } + } + else + { + Assert(false, ExcNotImplemented()); + return projected_point; + } + + return projected_point; + } +} // namespace GridTools +#endif // DOXYGEN + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/grid/CMakeLists.txt b/source/grid/CMakeLists.txt index 3d70ba8095..4c31e658aa 100644 --- a/source/grid/CMakeLists.txt +++ b/source/grid/CMakeLists.txt @@ -52,9 +52,10 @@ set(_separate_src grid_generator_pipe_junction.cc grid_in.cc grid_out.cc - grid_tools_cache.cc grid_tools.cc + grid_tools_cache.cc grid_tools_dof_handlers.cc + grid_tools_geometry.cc grid_tools_nontemplates.cc tria.cc ) @@ -78,8 +79,9 @@ set(_inst grid_out.inst.in grid_refinement.inst.in grid_tools.inst.in - grid_tools_dof_handlers.inst.in grid_tools_cache.inst.in + grid_tools_dof_handlers.inst.in + grid_tools_geometry.inst.in intergrid_map.inst.in manifold.inst.in manifold_lib.inst.in diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 6abda211d9..ba458f8fa2 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -106,338 +106,6 @@ DEAL_II_NAMESPACE_OPEN namespace GridTools { - template - double - diameter(const Triangulation &tria) - { - // we can't deal with distributed meshes since we don't have all - // vertices locally. there is one exception, however: if the mesh has - // never been refined. the way to test this is not to ask - // tria.n_levels()==1, since this is something that can happen on one - // processor without being true on all. however, we can ask for the - // global number of active cells and use that -#if defined(DEAL_II_WITH_P4EST) && defined(DEBUG) - if (const parallel::distributed::Triangulation *p_tria = - dynamic_cast< - const parallel::distributed::Triangulation *>(&tria)) - Assert(p_tria->n_global_active_cells() == tria.n_cells(0), - ExcNotImplemented()); -#endif - - // the algorithm used simply traverses all cells and picks out the - // boundary vertices. it may or may not be faster to simply get all - // vectors, don't mark boundary vertices, and compute the distances - // thereof, but at least as the mesh is refined, it seems better to - // first mark boundary nodes, as marking is O(N) in the number of - // cells/vertices, while computing the maximal distance is O(N*N) - const std::vector> &vertices = tria.get_vertices(); - std::vector boundary_vertices(vertices.size(), false); - - typename Triangulation::active_cell_iterator cell = - tria.begin_active(); - const typename Triangulation::active_cell_iterator endc = - tria.end(); - for (; cell != endc; ++cell) - for (const unsigned int face : cell->face_indices()) - if (cell->face(face)->at_boundary()) - for (unsigned int i = 0; i < cell->face(face)->n_vertices(); ++i) - boundary_vertices[cell->face(face)->vertex_index(i)] = true; - - // now traverse the list of boundary vertices and check distances. - // since distances are symmetric, we only have to check one half - double max_distance_sqr = 0; - std::vector::const_iterator pi = boundary_vertices.begin(); - const unsigned int N = boundary_vertices.size(); - for (unsigned int i = 0; i < N; ++i, ++pi) - { - std::vector::const_iterator pj = pi + 1; - for (unsigned int j = i + 1; j < N; ++j, ++pj) - if ((*pi == true) && (*pj == true) && - ((vertices[i] - vertices[j]).norm_square() > max_distance_sqr)) - max_distance_sqr = (vertices[i] - vertices[j]).norm_square(); - } - - return std::sqrt(max_distance_sqr); - } - - - - template - double - volume(const Triangulation &triangulation) - { - Assert(triangulation.get_reference_cells().size() == 1, - ExcNotImplemented()); - const ReferenceCell reference_cell = triangulation.get_reference_cells()[0]; - return volume( - triangulation, - reference_cell.template get_default_linear_mapping()); - } - - - - template - double - volume(const Triangulation &triangulation, - const Mapping &mapping) - { - // get the degree of the mapping if possible. if not, just assume 1 - unsigned int mapping_degree = 1; - if (const auto *p = dynamic_cast *>(&mapping)) - mapping_degree = p->get_degree(); - else if (const auto *p = - dynamic_cast *>(&mapping)) - mapping_degree = p->get_degree(); - - // then initialize an appropriate quadrature formula - Assert(triangulation.get_reference_cells().size() == 1, - ExcNotImplemented()); - const ReferenceCell reference_cell = triangulation.get_reference_cells()[0]; - const Quadrature quadrature_formula = - reference_cell.template get_gauss_type_quadrature(mapping_degree + - 1); - const unsigned int n_q_points = quadrature_formula.size(); - - // we really want the JxW values from the FEValues object, but it - // wants a finite element. create a cheap element as a dummy - // element - FE_Nothing dummy_fe(reference_cell); - FEValues fe_values(mapping, - dummy_fe, - quadrature_formula, - update_JxW_values); - - double local_volume = 0; - - // compute the integral quantities by quadrature - for (const auto &cell : triangulation.active_cell_iterators()) - if (cell->is_locally_owned()) - { - fe_values.reinit(cell); - for (unsigned int q = 0; q < n_q_points; ++q) - local_volume += fe_values.JxW(q); - } - - const double global_volume = - Utilities::MPI::sum(local_volume, triangulation.get_communicator()); - - return global_volume; - } - - - - 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, - const Triangulation &triangulation, - const Quadrature &quadrature) - { - FE_Nothing fe; - FEValues fe_values(mapping, fe, quadrature, update_jacobians); - - Vector aspect_ratio_vector(triangulation.n_active_cells()); - - // loop over cells of processor - for (const auto &cell : triangulation.active_cell_iterators()) - { - if (cell->is_locally_owned()) - { - double aspect_ratio_cell = 0.0; - - fe_values.reinit(cell); - - // loop over quadrature points - for (unsigned int q = 0; q < quadrature.size(); ++q) - { - const Tensor<2, dim, double> jacobian = - Tensor<2, dim, double>(fe_values.jacobian(q)); - - // We intentionally do not want to throw an exception in case of - // inverted elements since this is not the task of this - // function. Instead, inf is written into the vector in case of - // inverted elements. - if (determinant(jacobian) <= 0) - { - aspect_ratio_cell = std::numeric_limits::infinity(); - } - else - { - LAPACKFullMatrix J = LAPACKFullMatrix(dim); - for (unsigned int i = 0; i < dim; ++i) - for (unsigned int j = 0; j < dim; ++j) - J(i, j) = jacobian[i][j]; - - J.compute_svd(); - - const double max_sv = J.singular_value(0); - const double min_sv = J.singular_value(dim - 1); - const double ar = max_sv / min_sv; - - // Take the max between the previous and the current - // aspect ratio value; if we had previously encountered - // an inverted cell, we will have placed an infinity - // in the aspect_ratio_cell variable, and that value - // will survive this max operation. - aspect_ratio_cell = std::max(aspect_ratio_cell, ar); - } - } - - // fill vector - aspect_ratio_vector(cell->active_cell_index()) = aspect_ratio_cell; - } - } - - return aspect_ratio_vector; - } - - - - template - double - compute_maximum_aspect_ratio(const Mapping &mapping, - const Triangulation &triangulation, - const Quadrature &quadrature) - { - Vector aspect_ratio_vector = - compute_aspect_ratio_of_cells(mapping, triangulation, quadrature); - - return VectorTools::compute_global_error(triangulation, - aspect_ratio_vector, - VectorTools::Linfty_norm); - } - - - - template - BoundingBox - compute_bounding_box(const Triangulation &tria) - { - using iterator = - typename Triangulation::active_cell_iterator; - const auto predicate = [](const iterator &) { return true; }; - - return compute_bounding_box( - tria, std::function(predicate)); - } - - - // Generic functions for appending face data in 2d or 3d. TODO: we can // remove these once we have 'if constexpr'. namespace internal @@ -4442,40 +4110,6 @@ namespace GridTools - template - double - minimal_cell_diameter(const Triangulation &triangulation, - const Mapping &mapping) - { - double min_diameter = std::numeric_limits::max(); - for (const auto &cell : triangulation.active_cell_iterators()) - if (!cell->is_artificial()) - min_diameter = std::min(min_diameter, cell->diameter(mapping)); - - const double global_min_diameter = - Utilities::MPI::min(min_diameter, triangulation.get_communicator()); - return global_min_diameter; - } - - - - template - double - maximal_cell_diameter(const Triangulation &triangulation, - const Mapping &mapping) - { - double max_diameter = 0.; - for (const auto &cell : triangulation.active_cell_iterators()) - if (!cell->is_artificial()) - max_diameter = std::max(max_diameter, cell->diameter(mapping)); - - const double global_max_diameter = - Utilities::MPI::max(max_diameter, triangulation.get_communicator()); - return global_max_diameter; - } - - - namespace internal { namespace FixUpDistortedChildCells diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in index 446fd33c27..4824f73888 100644 --- a/source/grid/grid_tools.inst.in +++ b/source/grid/grid_tools.inst.in @@ -255,16 +255,6 @@ for (deal_II_space_dimension : SPACE_DIMENSIONS) GridTools::build_global_description_tree( const std::vector> &, const MPI_Comm); - - template Vector GridTools::compute_aspect_ratio_of_cells( - const Mapping &, - const Triangulation &, - const Quadrature &); - - template double GridTools::compute_maximum_aspect_ratio( - const Mapping &, - const Triangulation &, - const Quadrature &); } @@ -273,27 +263,6 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) #if deal_II_dimension <= deal_II_space_dimension namespace GridTools \{ - template double - diameter( - const Triangulation &); - - template double - volume(const Triangulation &); - - template double - 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 &); - template std::tuple>, std::vector>, SubCellData> @@ -441,18 +410,6 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) get_locally_owned_vertices( const Triangulation &); - template double - minimal_cell_diameter( - const Triangulation - &triangulation, - const Mapping &); - - template double - maximal_cell_diameter( - const Triangulation - &triangulation, - const Mapping &); - template std::map> get_all_vertices_at_boundary( const Triangulation &tria); diff --git a/source/grid/grid_tools_geometry.cc b/source/grid/grid_tools_geometry.cc new file mode 100644 index 0000000000..fd20d88573 --- /dev/null +++ b/source/grid/grid_tools_geometry.cc @@ -0,0 +1,410 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2001 - 2023 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#include +#include +#include +#include + +#include + +#include +#include +#include +#include + +#include +#include + +#include + +#include + +#include + +DEAL_II_NAMESPACE_OPEN + + +namespace GridTools +{ + template + double + diameter(const Triangulation &tria) + { + // we can't deal with distributed meshes since we don't have all + // vertices locally. there is one exception, however: if the mesh has + // never been refined. the way to test this is not to ask + // tria.n_levels()==1, since this is something that can happen on one + // processor without being true on all. however, we can ask for the + // global number of active cells and use that +#ifdef DEBUG + if (const auto *p_tria = dynamic_cast< + const parallel::DistributedTriangulationBase *>(&tria)) + Assert(p_tria->n_global_active_cells() == tria.n_cells(0), + ExcNotImplemented()); +#endif + + // the algorithm used simply traverses all cells and picks out the + // boundary vertices. it may or may not be faster to simply get all + // vectors, don't mark boundary vertices, and compute the distances + // thereof, but at least as the mesh is refined, it seems better to + // first mark boundary nodes, as marking is O(N) in the number of + // cells/vertices, while computing the maximal distance is O(N*N) + const std::vector> &vertices = tria.get_vertices(); + std::vector boundary_vertices(vertices.size(), false); + + typename Triangulation::active_cell_iterator cell = + tria.begin_active(); + const typename Triangulation::active_cell_iterator endc = + tria.end(); + for (; cell != endc; ++cell) + for (const unsigned int face : cell->face_indices()) + if (cell->face(face)->at_boundary()) + for (unsigned int i = 0; i < cell->face(face)->n_vertices(); ++i) + boundary_vertices[cell->face(face)->vertex_index(i)] = true; + + // now traverse the list of boundary vertices and check distances. + // since distances are symmetric, we only have to check one half + double max_distance_sqr = 0; + std::vector::const_iterator pi = boundary_vertices.begin(); + const unsigned int N = boundary_vertices.size(); + for (unsigned int i = 0; i < N; ++i, ++pi) + { + std::vector::const_iterator pj = pi + 1; + for (unsigned int j = i + 1; j < N; ++j, ++pj) + if ((*pi == true) && (*pj == true) && + ((vertices[i] - vertices[j]).norm_square() > max_distance_sqr)) + max_distance_sqr = (vertices[i] - vertices[j]).norm_square(); + } + + return std::sqrt(max_distance_sqr); + } + + + + template + double + volume(const Triangulation &triangulation) + { + Assert(triangulation.get_reference_cells().size() == 1, + ExcNotImplemented()); + const ReferenceCell reference_cell = triangulation.get_reference_cells()[0]; + return volume( + triangulation, + reference_cell.template get_default_linear_mapping()); + } + + + + template + double + volume(const Triangulation &triangulation, + const Mapping &mapping) + { + // get the degree of the mapping if possible. if not, just assume 1 + unsigned int mapping_degree = 1; + if (const auto *p = dynamic_cast *>(&mapping)) + mapping_degree = p->get_degree(); + else if (const auto *p = + dynamic_cast *>(&mapping)) + mapping_degree = p->get_degree(); + + // then initialize an appropriate quadrature formula + Assert(triangulation.get_reference_cells().size() == 1, + ExcNotImplemented()); + const ReferenceCell reference_cell = triangulation.get_reference_cells()[0]; + const Quadrature quadrature_formula = + reference_cell.template get_gauss_type_quadrature(mapping_degree + + 1); + const unsigned int n_q_points = quadrature_formula.size(); + + // we really want the JxW values from the FEValues object, but it + // wants a finite element. create a cheap element as a dummy + // element + FE_Nothing dummy_fe(reference_cell); + FEValues fe_values(mapping, + dummy_fe, + quadrature_formula, + update_JxW_values); + + double local_volume = 0; + + // compute the integral quantities by quadrature + for (const auto &cell : triangulation.active_cell_iterators()) + if (cell->is_locally_owned()) + { + fe_values.reinit(cell); + for (unsigned int q = 0; q < n_q_points; ++q) + local_volume += fe_values.JxW(q); + } + + const double global_volume = + Utilities::MPI::sum(local_volume, triangulation.get_communicator()); + + return global_volume; + } + + + + 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, + const Triangulation &triangulation, + const Quadrature &quadrature) + { + FE_Nothing fe; + FEValues fe_values(mapping, fe, quadrature, update_jacobians); + + Vector aspect_ratio_vector(triangulation.n_active_cells()); + + // loop over cells of processor + for (const auto &cell : triangulation.active_cell_iterators()) + { + if (cell->is_locally_owned()) + { + double aspect_ratio_cell = 0.0; + + fe_values.reinit(cell); + + // loop over quadrature points + for (unsigned int q = 0; q < quadrature.size(); ++q) + { + const Tensor<2, dim, double> jacobian = + Tensor<2, dim, double>(fe_values.jacobian(q)); + + // We intentionally do not want to throw an exception in case of + // inverted elements since this is not the task of this + // function. Instead, inf is written into the vector in case of + // inverted elements. + if (determinant(jacobian) <= 0) + { + aspect_ratio_cell = std::numeric_limits::infinity(); + } + else + { + LAPACKFullMatrix J = LAPACKFullMatrix(dim); + for (unsigned int i = 0; i < dim; ++i) + for (unsigned int j = 0; j < dim; ++j) + J(i, j) = jacobian[i][j]; + + J.compute_svd(); + + const double max_sv = J.singular_value(0); + const double min_sv = J.singular_value(dim - 1); + const double ar = max_sv / min_sv; + + // Take the max between the previous and the current + // aspect ratio value; if we had previously encountered + // an inverted cell, we will have placed an infinity + // in the aspect_ratio_cell variable, and that value + // will survive this max operation. + aspect_ratio_cell = std::max(aspect_ratio_cell, ar); + } + } + + // fill vector + aspect_ratio_vector(cell->active_cell_index()) = aspect_ratio_cell; + } + } + + return aspect_ratio_vector; + } + + + + template + double + compute_maximum_aspect_ratio(const Mapping &mapping, + const Triangulation &triangulation, + const Quadrature &quadrature) + { + Vector aspect_ratio_vector = + compute_aspect_ratio_of_cells(mapping, triangulation, quadrature); + + return VectorTools::compute_global_error(triangulation, + aspect_ratio_vector, + VectorTools::Linfty_norm); + } + + + + template + BoundingBox + compute_bounding_box(const Triangulation &tria) + { + using iterator = + typename Triangulation::active_cell_iterator; + const auto predicate = [](const iterator &) { return true; }; + + return compute_bounding_box( + tria, std::function(predicate)); + } + + + + template + double + minimal_cell_diameter(const Triangulation &triangulation, + const Mapping &mapping) + { + double min_diameter = std::numeric_limits::max(); + for (const auto &cell : triangulation.active_cell_iterators()) + if (!cell->is_artificial()) + min_diameter = std::min(min_diameter, cell->diameter(mapping)); + + const double global_min_diameter = + Utilities::MPI::min(min_diameter, triangulation.get_communicator()); + return global_min_diameter; + } + + + + template + double + maximal_cell_diameter(const Triangulation &triangulation, + const Mapping &mapping) + { + double max_diameter = 0.; + for (const auto &cell : triangulation.active_cell_iterators()) + if (!cell->is_artificial()) + max_diameter = std::max(max_diameter, cell->diameter(mapping)); + + const double global_max_diameter = + Utilities::MPI::max(max_diameter, triangulation.get_communicator()); + return global_max_diameter; + } +} /* namespace GridTools */ + + +// explicit instantiations +#include "grid_tools_geometry.inst" + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/grid/grid_tools_geometry.inst.in b/source/grid/grid_tools_geometry.inst.in new file mode 100644 index 0000000000..5695304d95 --- /dev/null +++ b/source/grid/grid_tools_geometry.inst.in @@ -0,0 +1,73 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2001 - 2023 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) + { +#if deal_II_dimension <= deal_II_space_dimension + namespace GridTools + \{ + template double + diameter( + const Triangulation &); + + template double + volume(const Triangulation &); + + template double + volume(const Triangulation &, + const Mapping &); + + template double + minimal_cell_diameter( + const Triangulation + &triangulation, + const Mapping &); + + template double + maximal_cell_diameter( + const Triangulation + &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 &); + \} +#endif + } + +for (deal_II_space_dimension : SPACE_DIMENSIONS) + { + namespace GridTools + \{ + template Vector + GridTools::compute_aspect_ratio_of_cells( + const Mapping &, + const Triangulation &, + const Quadrature &); + + template double + GridTools::compute_maximum_aspect_ratio( + const Mapping &, + const Triangulation &, + const Quadrature &); + \} + } diff --git a/source/grid/tria_accessor.cc b/source/grid/tria_accessor.cc index 8723a402af..0836ac904c 100644 --- a/source/grid/tria_accessor.cc +++ b/source/grid/tria_accessor.cc @@ -21,7 +21,7 @@ #include #include -#include +#include #include #include #include