From decee3c427b841309e3e0f3fe8d74b72b902a8f3 Mon Sep 17 00:00:00 2001 From: Marco Feder Date: Sat, 18 Jun 2022 13:06:54 +0200 Subject: [PATCH] compute_affine_transformation in codim1 and codim2 case --- include/deal.II/base/quadrature_lib.h | 12 +++- source/base/quadrature_lib.cc | 43 ++++++++++---- tests/base/quadrature_simplex_01_codim1.cc | 59 +++++++++++++++++++ .../base/quadrature_simplex_01_codim1.output | 13 ++++ tests/base/quadrature_simplex_02_codim2.cc | 56 ++++++++++++++++++ .../base/quadrature_simplex_02_codim2.output | 7 +++ 6 files changed, 175 insertions(+), 15 deletions(-) create mode 100644 tests/base/quadrature_simplex_01_codim1.cc create mode 100644 tests/base/quadrature_simplex_01_codim1.output create mode 100644 tests/base/quadrature_simplex_02_codim2.cc create mode 100644 tests/base/quadrature_simplex_02_codim2.output diff --git a/include/deal.II/base/quadrature_lib.h b/include/deal.II/base/quadrature_lib.h index 87b2718a06..6f1ad61c53 100644 --- a/include/deal.II/base/quadrature_lib.h +++ b/include/deal.II/base/quadrature_lib.h @@ -609,6 +609,7 @@ public: QSimplex(const Quadrature &quad); /** + * * Return an affine transformation of this quadrature, that can be used to * integrate on the simplex identified by `vertices`. * @@ -626,14 +627,19 @@ public: * that is $J \dealcoloneq |\text{det}(B)|$. If $J$ is zero, an empty * quadrature is returned. This may happen, in two dimensions, if the three * vertices are aligned, or in three dimensions if the four vertices are on - * the same plane. + * the same plane. The present function works also in the codimension one and + * codimension two case. For instance, when `dim=2` and `spacedim=3`, we can + * map the quadrature points so that they live on the physical triangle + * embedded in the three dimensional space. In such a case, the matrix $B$ is + * not square anymore. * * @param[in] vertices The vertices of the simplex you wish to integrate on * @return A quadrature object that can be used to integrate on the simplex */ - Quadrature + template + Quadrature compute_affine_transformation( - const std::array, dim + 1> &vertices) const; + const std::array, dim + 1> &vertices) const; }; /** diff --git a/source/base/quadrature_lib.cc b/source/base/quadrature_lib.cc index 368cd980fe..bfe1d52b23 100644 --- a/source/base/quadrature_lib.cc +++ b/source/base/quadrature_lib.cc @@ -1227,31 +1227,35 @@ QSimplex::QSimplex(const Quadrature &quad) template -Quadrature +template +Quadrature QSimplex::compute_affine_transformation( - const std::array, dim + 1> &vertices) const + const std::array, dim + 1> &vertices) const { - Tensor<2, dim> B; + Assert(dim <= spacedim, + ExcMessage("Invalid combination of dim and spacedim .")); + DerivativeForm<1, spacedim, dim> Bt; for (unsigned int d = 0; d < dim; ++d) - B[d] = vertices[d + 1] - vertices[0]; + Bt[d] = vertices[d + 1] - vertices[0]; - B = transpose(B); - const double J = std::abs(determinant(B)); + const auto B = Bt.transpose(); + const double J = std::abs(B.determinant()); // if the determinant is zero, we return an empty quadrature if (J < 1e-12) - return Quadrature(); + return Quadrature(); - std::vector> qp(this->size()); - std::vector w(this->size()); + std::vector> qp(this->size()); + std::vector w(this->size()); for (unsigned int i = 0; i < this->size(); ++i) { - qp[i] = Point(vertices[0] + B * this->point(i)); - w[i] = J * this->weight(i); + qp[i] = + Point(vertices[0] + apply_transformation(B, this->point(i))); + w[i] = J * this->weight(i); } - return Quadrature(qp, w); + return Quadrature(qp, w); } @@ -2208,3 +2212,18 @@ template class QWitherdenVincentSimplex<2>; template class QWitherdenVincentSimplex<3>; DEAL_II_NAMESPACE_CLOSE + +namespace dealii +{ + template Quadrature<2> + QSimplex<1>::compute_affine_transformation( + const std::array, 1 + 1> &vertices) const; + + template Quadrature<3> + QSimplex<1>::compute_affine_transformation( + const std::array, 1 + 1> &vertices) const; + + template Quadrature<3> + QSimplex<2>::compute_affine_transformation( + const std::array, 2 + 1> &vertices) const; +} // namespace dealii diff --git a/tests/base/quadrature_simplex_01_codim1.cc b/tests/base/quadrature_simplex_01_codim1.cc new file mode 100644 index 0000000000..e50ff089ed --- /dev/null +++ b/tests/base/quadrature_simplex_01_codim1.cc @@ -0,0 +1,59 @@ +// --------------------------------------------------------------------- + +// Copyright (C) 2022 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. + +// --------------------------------------------------------------------- + +// construct a simplex quadrature, and check that the affine transformation +// gives the correct results in the codimension one case + +#include + +#include + +#include "../tests.h" + +// #include "simplex.h" + +template +void +test(const int degree, const std::array, spacedim> &vertices) +{ + QGaussSimplex quad(degree); + + deallog << "# dim = " << dim << std::endl; + deallog << "# spacedim = " << spacedim << std::endl; + + auto quad2 = quad.compute_affine_transformation(vertices); + + for (auto p : quad2.get_points()) + deallog << p << std::endl; + + deallog << std::endl + << "# Area: " << std::setprecision(15) + << std::accumulate(quad2.get_weights().begin(), + quad2.get_weights().end(), + 0.0) + << std::endl + << std::endl; +} + + +int +main() +{ + initlog(); + + test<1, 2>(1, {{Point<2>(0., 2.), Point<2>(1., 0.)}}); + test<2, 3>( + 1, {{Point<3>(0., 2., 0.), Point<3>(1., 2., 4.), Point<3>(1., 1., 1.)}}); +} diff --git a/tests/base/quadrature_simplex_01_codim1.output b/tests/base/quadrature_simplex_01_codim1.output new file mode 100644 index 0000000000..268045a994 --- /dev/null +++ b/tests/base/quadrature_simplex_01_codim1.output @@ -0,0 +1,13 @@ + +DEAL::# dim = 1 +DEAL::# spacedim = 2 +DEAL::0.500000 1.00000 +DEAL:: +DEAL::# Area: 2.23606797749979 +DEAL:: +DEAL::# dim = 2 +DEAL::# spacedim = 3 +DEAL::0.666666666666667 1.66666666666667 1.66666666666667 +DEAL:: +DEAL::# Area: 2.54950975679639 +DEAL:: diff --git a/tests/base/quadrature_simplex_02_codim2.cc b/tests/base/quadrature_simplex_02_codim2.cc new file mode 100644 index 0000000000..be07845df4 --- /dev/null +++ b/tests/base/quadrature_simplex_02_codim2.cc @@ -0,0 +1,56 @@ +// --------------------------------------------------------------------- + +// Copyright (C) 2022 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. + +// --------------------------------------------------------------------- + +// construct a simplex quadrature, and check that the affine transformation +// gives the correct results in the codimension one case + +#include + +#include + +#include "../tests.h" + +// #include "simplex.h" + +template +void +test(const int degree, const std::array, dim + 1> &vertices) +{ + QGaussSimplex quad(degree); + + deallog << "# dim = " << dim << std::endl; + deallog << "# spacedim = " << spacedim << std::endl; + + auto quad2 = quad.compute_affine_transformation(vertices); + + for (auto p : quad2.get_points()) + deallog << p << std::endl; + + deallog << std::endl + << "# Area: " << std::setprecision(15) + << std::accumulate(quad2.get_weights().begin(), + quad2.get_weights().end(), + 0.0) + << std::endl + << std::endl; +} + + +int +main() +{ + initlog(); + test<1, 3>(1, {{Point<3>(0., 1., 0.), Point<3>(0., 0., 1.)}}); +} diff --git a/tests/base/quadrature_simplex_02_codim2.output b/tests/base/quadrature_simplex_02_codim2.output new file mode 100644 index 0000000000..afcbea4275 --- /dev/null +++ b/tests/base/quadrature_simplex_02_codim2.output @@ -0,0 +1,7 @@ + +DEAL::# dim = 1 +DEAL::# spacedim = 3 +DEAL::0.00000 0.500000 0.500000 +DEAL:: +DEAL::# Area: 1.41421356237310 +DEAL:: -- 2.39.5