--- /dev/null
+New: The function GridTools::compute_maximum_aspect_ratio() computes the
+maximum aspect ratio of a mesh with arbitrarily deformed elements.
+<br>
+(Niklas Fehn, 2019/10/08)
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 <int dim>
+ Vector<double>
+ compute_aspect_ratio_of_cells(const Triangulation<dim> &triangulation,
+ const Mapping<dim> & mapping,
+ const Quadrature<dim> & 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 <int dim>
+ double
+ compute_maximum_aspect_ratio(const Triangulation<dim> &triangulation,
+ const Mapping<dim> & mapping,
+ const Quadrature<dim> & quadrature);
+
/**
* Compute the smallest box containing the entire triangulation.
*
#include <deal.II/lac/vector_memory.h>
#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/uniform_real_distribution.hpp>
}
+ template <int dim>
+ Vector<double>
+ compute_aspect_ratio_of_cells(const Triangulation<dim> &triangulation,
+ const Mapping<dim> & mapping,
+ const Quadrature<dim> & quadrature)
+ {
+ FE_Nothing<dim> fe;
+ FEValues<dim> fe_values(mapping, fe, quadrature, update_jacobians);
+
+ Vector<double> 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<double> J = LAPACKFullMatrix<double>(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<double>::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 <int dim>
+ double
+ compute_maximum_aspect_ratio(const Triangulation<dim> &triangulation,
+ const Mapping<dim> & mapping,
+ const Quadrature<dim> & quadrature)
+ {
+ Vector<double> aspect_ratio_vector =
+ compute_aspect_ratio_of_cells(triangulation, mapping, quadrature);
+
+ return VectorTools::compute_global_error(triangulation,
+ aspect_ratio_vector,
+ VectorTools::Linfty_norm);
+ }
+
template <int dim, int spacedim>
BoundingBox<spacedim>
std::pair<BoundingBox<deal_II_space_dimension>, unsigned int>>
GridTools::build_global_description_tree(
const std::vector<BoundingBox<deal_II_space_dimension>> &, MPI_Comm);
+
+ template Vector<double> GridTools::compute_aspect_ratio_of_cells(
+ const Triangulation<deal_II_space_dimension> &,
+ const Mapping<deal_II_space_dimension> &,
+ const Quadrature<deal_II_space_dimension> &);
+
+ template double GridTools::compute_maximum_aspect_ratio(
+ const Triangulation<deal_II_space_dimension> &,
+ const Mapping<deal_II_space_dimension> &,
+ const Quadrature<deal_II_space_dimension> &);
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/exceptions.h>
+#include <deal.II/base/mpi.h>
+#include <deal.II/base/partitioner.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/mapping_q_generic.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+
+#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 <int dim>
+double
+compute_aspect_ratio_hyper_rectangle(
+ Point<dim> const & left,
+ Point<dim> const & right,
+ std::vector<unsigned int> const &refinements,
+ unsigned int degree = 1,
+ unsigned int n_q_points = 2,
+ bool deform = false,
+ double factor = 1.0)
+{
+ Triangulation<dim> tria;
+ GridGenerator::subdivided_hyper_rectangle(
+ tria, refinements, left, right, false);
+
+ if (deform)
+ {
+ Tensor<1, dim> diag = right - left;
+ double l = diag.norm();
+ Point<dim> 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<dim> const mapping(degree);
+ QGauss<dim> const gauss(n_q_points);
+
+ Vector<double> 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<unsigned int> 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<unsigned int> 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<unsigned int> 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;
+}
--- /dev/null
+
+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
projects / dealii.git / commitdiff