From 9733420f7aaf20104730ee44492e57e720a761fb Mon Sep 17 00:00:00 2001 From: Niklas Fehn Date: Tue, 8 Oct 2019 13:10:56 +0200 Subject: [PATCH] new functionality to compute mesh aspect ratio --- doc/news/changes/minor/20191009NiklasFehn | 4 + include/deal.II/grid/grid_tools.h | 49 +++++ source/grid/grid_tools.cc | 74 +++++++ source/grid/grid_tools.inst.in | 10 + tests/grid/grid_tools_aspect_ratio.cc | 234 ++++++++++++++++++++++ tests/grid/grid_tools_aspect_ratio.output | 109 ++++++++++ 6 files changed, 480 insertions(+) create mode 100644 doc/news/changes/minor/20191009NiklasFehn create mode 100644 tests/grid/grid_tools_aspect_ratio.cc create mode 100644 tests/grid/grid_tools_aspect_ratio.output diff --git a/doc/news/changes/minor/20191009NiklasFehn b/doc/news/changes/minor/20191009NiklasFehn new file mode 100644 index 0000000000..fe96e57c03 --- /dev/null +++ b/doc/news/changes/minor/20191009NiklasFehn @@ -0,0 +1,4 @@ +New: The function GridTools::compute_maximum_aspect_ratio() computes the +maximum aspect ratio of a mesh with arbitrarily deformed elements. +
+(Niklas Fehn, 2019/10/08) diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index 307c305dbe..8cb2bcae5b 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -230,6 +230,55 @@ namespace GridTools double cell_measure(const T &, ...); + /** + * 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. + * + * @author Niklas Fehn, Martin Kronbichler, 2019 + */ + template + Vector + compute_aspect_ratio_of_cells(const Triangulation &triangulation, + const Mapping & mapping, + 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. + * + * @author Niklas Fehn, Martin Kronbichler, 2019 + */ + template + double + compute_maximum_aspect_ratio(const Triangulation &triangulation, + const Mapping & mapping, + const Quadrature & quadrature); + /** * Compute the smallest box containing the entire triangulation. * diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index bbbaeecc20..b824ba3192 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -52,6 +52,7 @@ #include #include +#include #include #include @@ -414,6 +415,79 @@ namespace GridTools } + template + Vector + compute_aspect_ratio_of_cells(const Triangulation &triangulation, + const Mapping & mapping, + 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) + { + Tensor<2, dim, double> jacobian = + Tensor<2, dim, double>(fe_values.jacobian(q)); + + 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]; + + // 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 + { + J.compute_svd(); + + double const max_sv = J.singular_value(0); + double const min_sv = J.singular_value(dim - 1); + double const ar = max_sv / min_sv; + + 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 Triangulation &triangulation, + const Mapping & mapping, + const Quadrature & quadrature) + { + Vector aspect_ratio_vector = + compute_aspect_ratio_of_cells(triangulation, mapping, quadrature); + + return VectorTools::compute_global_error(triangulation, + aspect_ratio_vector, + VectorTools::Linfty_norm); + } + template BoundingBox diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in index 900a1caf99..edf8ab5280 100644 --- a/source/grid/grid_tools.inst.in +++ b/source/grid/grid_tools.inst.in @@ -175,6 +175,16 @@ for (deal_II_space_dimension : SPACE_DIMENSIONS) std::pair, unsigned int>> GridTools::build_global_description_tree( const std::vector> &, MPI_Comm); + + template Vector GridTools::compute_aspect_ratio_of_cells( + const Triangulation &, + const Mapping &, + const Quadrature &); + + template double GridTools::compute_maximum_aspect_ratio( + const Triangulation &, + const Mapping &, + const Quadrature &); } diff --git a/tests/grid/grid_tools_aspect_ratio.cc b/tests/grid/grid_tools_aspect_ratio.cc new file mode 100644 index 0000000000..54507c9b35 --- /dev/null +++ b/tests/grid/grid_tools_aspect_ratio.cc @@ -0,0 +1,234 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 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. +// +// --------------------------------------------------------------------- + +// Short test to validate GridTools::compute_maximum_aspect_ratio() + +#include +#include +#include +#include +#include + +#include + +#include + +#include +#include + +#include +#include + +#include "../tests.h" + +using namespace dealii; + +/* + * The parameters left, right, refinements follow the hyper-rectangle + * notation. + * + * degree denotes the polynomial degree of the mapping. + * + * n_q_points denotes the number of 1D quadrature points used during + * the computation of the aspect ratio. + * + * If deform is true, the first vertex of the first element of the + * hyper_rectangle is shifted, with factor being a scaling factor to + * be able to generate inverted elements. + */ +template +double +compute_aspect_ratio_hyper_rectangle( + Point const & left, + Point const & right, + std::vector const &refinements, + unsigned int degree = 1, + unsigned int n_q_points = 2, + bool deform = false, + double factor = 1.0) +{ + Triangulation tria; + GridGenerator::subdivided_hyper_rectangle( + tria, refinements, left, right, false); + + if (deform) + { + Tensor<1, dim> diag = right - left; + double l = diag.norm(); + Point shift; + for (unsigned int d = 0; d < dim; ++d) + shift[d] = l * factor * (0.05 + d * 0.01); + tria.begin_active()->vertex(0) += shift; + } + + MappingQGeneric const mapping(degree); + QGauss const gauss(n_q_points); + + Vector ratios = + GridTools::compute_aspect_ratio_of_cells(tria, mapping, gauss); + + deallog << std::endl + << "Parameters:" + << " d = " << dim << ", degree = " << degree + << ", n_q_points = " << n_q_points << ", deform = " << deform + << ", factor = " << factor << std::endl; + + deallog << "Aspect ratio vector = " << std::endl; + ratios.print(deallog.get_file_stream()); + + return GridTools::compute_maximum_aspect_ratio(tria, mapping, gauss); +} + +int +main(int argc, char **argv) +{ + initlog(); + deallog << std::setprecision(6); + + try + { + Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1); + + // 1D + { + Point<1> left = Point<1>(0.0); + Point<1> right = Point<1>(1.0); + + double ar = 0.0; + + std::vector refine(1, 2); + + ar = compute_aspect_ratio_hyper_rectangle(left, right, refine, 1, 2); + deallog << "aspect ratio max = " << ar << std::endl << std::endl; + + refine[0] = 5; + + ar = compute_aspect_ratio_hyper_rectangle(left, right, refine, 1, 2); + deallog << "aspect ratio max = " << ar << std::endl << std::endl; + } + + // 2D + { + Point<2> left = Point<2>(0.0, 0.0); + Point<2> right = Point<2>(1.0, 1.0); + + double ar = 0.0; + + std::vector refine(2, 2); + + ar = compute_aspect_ratio_hyper_rectangle(left, right, refine, 1, 2); + deallog << "aspect ratio max = " << ar << std::endl << std::endl; + + ar = + compute_aspect_ratio_hyper_rectangle(left, right, refine, 1, 2, true); + deallog << "aspect ratio max = " << ar << std::endl << std::endl; + + ar = compute_aspect_ratio_hyper_rectangle( + left, right, refine, 1, 2, true, 10); + deallog << "aspect ratio max = " << ar << std::endl << std::endl; + + refine[0] = 5; + + ar = compute_aspect_ratio_hyper_rectangle(left, right, refine, 1, 2); + deallog << "aspect ratio max = " << ar << std::endl << std::endl; + + ar = + compute_aspect_ratio_hyper_rectangle(left, right, refine, 1, 2, true); + deallog << "aspect ratio max = " << ar << std::endl << std::endl; + + ar = + compute_aspect_ratio_hyper_rectangle(left, right, refine, 1, 5, true); + deallog << "aspect ratio max = " << ar << std::endl << std::endl; + + ar = compute_aspect_ratio_hyper_rectangle( + left, right, refine, 1, 15, true); + deallog << "aspect ratio max = " << ar << std::endl << std::endl; + + ar = compute_aspect_ratio_hyper_rectangle( + left, right, refine, 1, 2, true, 10); + deallog << "aspect ratio max = " << ar << std::endl << std::endl; + } + + // 3D + { + Point<3> left = Point<3>(0.0, 0.0, 0.0); + Point<3> right = Point<3>(1.0, 1.0, 1.0); + + double ar = 0.0; + + std::vector refine(3, 2); + + ar = compute_aspect_ratio_hyper_rectangle(left, right, refine, 1, 2); + deallog << "aspect ratio max = " << ar << std::endl << std::endl; + + ar = + compute_aspect_ratio_hyper_rectangle(left, right, refine, 1, 2, true); + deallog << "aspect ratio max = " << ar << std::endl << std::endl; + + ar = compute_aspect_ratio_hyper_rectangle( + left, right, refine, 1, 2, true, 10); + deallog << "aspect ratio max = " << ar << std::endl << std::endl; + + refine[0] = 5; + refine[1] = 3; + + ar = compute_aspect_ratio_hyper_rectangle(left, right, refine, 1, 2); + deallog << "aspect ratio max = " << ar << std::endl << std::endl; + + ar = + compute_aspect_ratio_hyper_rectangle(left, right, refine, 1, 2, true); + deallog << "aspect ratio max = " << ar << std::endl << std::endl; + + ar = + compute_aspect_ratio_hyper_rectangle(left, right, refine, 1, 5, true); + deallog << "aspect ratio max = " << ar << std::endl << std::endl; + + ar = compute_aspect_ratio_hyper_rectangle( + left, right, refine, 1, 15, true); + deallog << "aspect ratio max = " << ar << std::endl << std::endl; + + ar = compute_aspect_ratio_hyper_rectangle( + left, right, refine, 1, 2, true, 10); + deallog << "aspect ratio max = " << ar << std::endl << std::endl; + } + } + catch (std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +} diff --git a/tests/grid/grid_tools_aspect_ratio.output b/tests/grid/grid_tools_aspect_ratio.output new file mode 100644 index 0000000000..395806aec4 --- /dev/null +++ b/tests/grid/grid_tools_aspect_ratio.output @@ -0,0 +1,109 @@ + +DEAL:: +DEAL::Parameters: d = 1, degree = 1, n_q_points = 2, deform = 0, factor = 1.00000 +DEAL::Aspect ratio vector = +1.000e+00 1.000e+00 +DEAL::aspect ratio max = 1.00000 +DEAL:: +DEAL:: +DEAL::Parameters: d = 1, degree = 1, n_q_points = 2, deform = 0, factor = 1.00000 +DEAL::Aspect ratio vector = +1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 +DEAL::aspect ratio max = 1.00000 +DEAL:: +DEAL:: +DEAL::Parameters: d = 2, degree = 1, n_q_points = 2, deform = 0, factor = 1.00000 +DEAL::Aspect ratio vector = +1.000e+00 1.000e+00 1.000e+00 1.000e+00 +DEAL::aspect ratio max = 1.00000 +DEAL:: +DEAL:: +DEAL::Parameters: d = 2, degree = 1, n_q_points = 2, deform = 1, factor = 1.00000 +DEAL::Aspect ratio vector = +1.327e+00 1.000e+00 1.000e+00 1.000e+00 +DEAL::aspect ratio max = 1.32670 +DEAL:: +DEAL:: +DEAL::Parameters: d = 2, degree = 1, n_q_points = 2, deform = 1, factor = 10.0000 +DEAL::Aspect ratio vector = +inf 1.000e+00 1.000e+00 1.000e+00 +DEAL::aspect ratio max = inf +DEAL:: +DEAL:: +DEAL::Parameters: d = 2, degree = 1, n_q_points = 2, deform = 0, factor = 1.00000 +DEAL::Aspect ratio vector = +2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 +DEAL::aspect ratio max = 2.50000 +DEAL:: +DEAL:: +DEAL::Parameters: d = 2, degree = 1, n_q_points = 2, deform = 1, factor = 1.00000 +DEAL::Aspect ratio vector = +3.476e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 +DEAL::aspect ratio max = 3.47552 +DEAL:: +DEAL:: +DEAL::Parameters: d = 2, degree = 1, n_q_points = 5, deform = 1, factor = 1.00000 +DEAL::Aspect ratio vector = +3.866e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 +DEAL::aspect ratio max = 3.86561 +DEAL:: +DEAL:: +DEAL::Parameters: d = 2, degree = 1, n_q_points = 15, deform = 1, factor = 1.00000 +DEAL::Aspect ratio vector = +3.971e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 +DEAL::aspect ratio max = 3.97069 +DEAL:: +DEAL:: +DEAL::Parameters: d = 2, degree = 1, n_q_points = 2, deform = 1, factor = 10.0000 +DEAL::Aspect ratio vector = +inf 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 +DEAL::aspect ratio max = inf +DEAL:: +DEAL:: +DEAL::Parameters: d = 3, degree = 1, n_q_points = 2, deform = 0, factor = 1.00000 +DEAL::Aspect ratio vector = +1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 +DEAL::aspect ratio max = 1.00000 +DEAL:: +DEAL:: +DEAL::Parameters: d = 3, degree = 1, n_q_points = 2, deform = 1, factor = 1.00000 +DEAL::Aspect ratio vector = +1.641e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 +DEAL::aspect ratio max = 1.64083 +DEAL:: +DEAL:: +DEAL::Parameters: d = 3, degree = 1, n_q_points = 2, deform = 1, factor = 10.0000 +DEAL::Aspect ratio vector = +inf 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 +DEAL::aspect ratio max = inf +DEAL:: +DEAL:: +DEAL::Parameters: d = 3, degree = 1, n_q_points = 2, deform = 0, factor = 1.00000 +DEAL::Aspect ratio vector = +2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 +DEAL::aspect ratio max = 2.50000 +DEAL:: +DEAL:: +DEAL::Parameters: d = 3, degree = 1, n_q_points = 2, deform = 1, factor = 1.00000 +DEAL::Aspect ratio vector = +4.611e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 +DEAL::aspect ratio max = 4.61134 +DEAL:: +DEAL:: +DEAL::Parameters: d = 3, degree = 1, n_q_points = 5, deform = 1, factor = 1.00000 +DEAL::Aspect ratio vector = +1.613e+01 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 +DEAL::aspect ratio max = 16.1349 +DEAL:: +DEAL:: +DEAL::Parameters: d = 3, degree = 1, n_q_points = 15, deform = 1, factor = 1.00000 +DEAL::Aspect ratio vector = +6.696e+01 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 +DEAL::aspect ratio max = 66.9573 +DEAL:: +DEAL:: +DEAL::Parameters: d = 3, degree = 1, n_q_points = 2, deform = 1, factor = 10.0000 +DEAL::Aspect ratio vector = +inf 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 2.500e+00 +DEAL::aspect ratio max = inf +DEAL:: \ No newline at end of file -- 2.39.5