From: Yann Jobic Date: Tue, 10 Dec 2024 08:06:01 +0000 (+0100) Subject: add a tolerance parameter with a default value to the match of periodic nodes. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F17921%2Fhead;p=dealii.git add a tolerance parameter with a default value to the match of periodic nodes. modif variable name and more explanatory comments add the test to the new collect_periodic_faces modification of the test source file to include a distort mesh, and to possibly make the test fail, in debug mode add a sequential test program, and add the comments of what the tests are doing add initlog and the output file remove parallel version of the test files modified by the indent script --- diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index df298799ca..79941f36eb 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -2358,11 +2358,13 @@ namespace GridTools * If no such relation exists then the returned std::optional object is empty * (i.e., has_value() will return `false`). * - * Here, two vertices v_1 and v_2 are considered equal, if - * $M\cdot v_1 + offset - v_2$ is parallel to the unit vector in unit - * direction @p direction. If the parameter @p matrix is a reference to a - * spacedim x spacedim matrix, $M$ is set to @p matrix, otherwise $M$ is the - * identity matrix. + * Here, two vertices v_1 and v_2 are considered equal, to + * the nearest tolerance, e.g. distance < abs_tol, if $M\cdot v_1 + offset - + * v_2$ is parallel to the + * unit vector in unit direction @p direction. Note that this tolerance is not + * scaled with the mesh, it's an absolute tolerance. If the parameter @p matrix + * is a reference to a spacedim x spacedim matrix, $M$ is set to @p matrix, + * otherwise $M$ is the identity matrix. * * If the matching was successful, the _relative_ orientation of @p face1 with * respect to @p face2 is returned a @@ -2379,7 +2381,8 @@ namespace GridTools const unsigned int direction, const Tensor<1, FaceIterator::AccessorType::space_dimension> &offset = Tensor<1, FaceIterator::AccessorType::space_dimension>(), - const FullMatrix &matrix = FullMatrix()); + const FullMatrix &matrix = FullMatrix(), + const double abs_tol = 1e-10); /** * This function will collect periodic face pairs on the coarsest mesh level @@ -2403,8 +2406,9 @@ namespace GridTools * * The @p offset is a vector tangential to the faces that is added to the * location of vertices of the 'first' boundary when attempting to match - * them to the corresponding vertices of the 'second' boundary. This can be - * used to implement conditions such as $u(0,y)=u(1,y+1)$. + * them (to the nearest absolute tolerance, distance < abs_tol) to the + * corresponding vertices of the 'second' boundary. This can be used to + * implement conditions such as $u(0,y)=u(1,y+1)$. * * Optionally, a $dim\times dim$ rotation @p matrix can be specified that * describes how vector valued DoFs of the first face should be modified @@ -2450,7 +2454,8 @@ namespace GridTools &matched_pairs, const Tensor<1, MeshType::space_dimension> &offset = dealii::Tensor<1, MeshType::space_dimension>(), - const FullMatrix &matrix = FullMatrix()); + const FullMatrix &matrix = FullMatrix(), + const double abs_tol = 1e-10); /** @@ -2483,7 +2488,8 @@ namespace GridTools &matched_pairs, const dealii::Tensor<1, MeshType::space_dimension> &offset = dealii::Tensor<1, MeshType::space_dimension>(), - const FullMatrix &matrix = FullMatrix()); + const FullMatrix &matrix = FullMatrix(), + const double abs_tol = 1e-10); /** @} */ /** diff --git a/source/grid/grid_tools_dof_handlers.cc b/source/grid/grid_tools_dof_handlers.cc index 3f7da7e914..1a27a0e588 100644 --- a/source/grid/grid_tools_dof_handlers.cc +++ b/source/grid/grid_tools_dof_handlers.cc @@ -2127,7 +2127,8 @@ namespace GridTools std::vector> &matched_pairs, const dealii::Tensor<1, CellIterator::AccessorType::space_dimension> &offset, - const FullMatrix &matrix) + const FullMatrix &matrix, + const double abs_tol = 1e-10) { static const int space_dim = CellIterator::AccessorType::space_dimension; (void)space_dim; @@ -2173,7 +2174,8 @@ namespace GridTools cell2->face(face_idx2), direction, offset, - matrix)) + matrix, + abs_tol)) { // We have a match, so insert the matching pairs and // remove the matched cell in pairs2 to speed up the @@ -2223,7 +2225,8 @@ namespace GridTools std::vector> &matched_pairs, const Tensor<1, MeshType::space_dimension> &offset, - const FullMatrix &matrix) + const FullMatrix &matrix, + const double abs_tol) { static const int dim = MeshType::dimension; static const int space_dim = MeshType::space_dimension; @@ -2275,7 +2278,7 @@ namespace GridTools // and call match_periodic_face_pairs that does the actual matching: match_periodic_face_pairs( - pairs1, pairs2, direction, matched_pairs, offset, matrix); + pairs1, pairs2, direction, matched_pairs, offset, matrix, abs_tol); if constexpr (running_in_debug_mode()) { @@ -2305,7 +2308,8 @@ namespace GridTools std::vector> &matched_pairs, const Tensor<1, MeshType::space_dimension> &offset, - const FullMatrix &matrix) + const FullMatrix &matrix, + const double abs_tol) { static const int dim = MeshType::dimension; static const int space_dim = MeshType::space_dimension; @@ -2370,7 +2374,7 @@ namespace GridTools // and call match_periodic_face_pairs that does the actual matching: match_periodic_face_pairs( - pairs1, pairs2, direction, matched_pairs, offset, matrix); + pairs1, pairs2, direction, matched_pairs, offset, matrix, abs_tol); } @@ -2390,7 +2394,8 @@ namespace GridTools const Point &point2, const unsigned int direction, const Tensor<1, spacedim> &offset, - const FullMatrix &matrix) + const FullMatrix &matrix, + const double abs_tol = 1e-10) { AssertIndexRange(direction, spacedim); @@ -2413,7 +2418,7 @@ namespace GridTools if (i == direction) continue; - if (std::abs(distance[i]) > 1.e-10) + if (std::abs(distance[i]) > abs_tol) return false; } @@ -2429,7 +2434,8 @@ namespace GridTools const FaceIterator &face2, const unsigned int direction, const Tensor<1, FaceIterator::AccessorType::space_dimension> &offset, - const FullMatrix &matrix) + const FullMatrix &matrix, + const double abs_tol) { Assert(matrix.m() == matrix.n(), ExcMessage("The supplied matrix must be a square matrix")); @@ -2458,7 +2464,8 @@ namespace GridTools face2->vertex(*it), direction, offset, - matrix)) + matrix, + abs_tol)) { face1_vertices[i] = *it; face2_vertices[i] = i; diff --git a/source/grid/grid_tools_dof_handlers.inst.in b/source/grid/grid_tools_dof_handlers.inst.in index f03a2ac367..194530d489 100644 --- a/source/grid/grid_tools_dof_handlers.inst.in +++ b/source/grid/grid_tools_dof_handlers.inst.in @@ -255,7 +255,8 @@ for (X : SEQUENTIAL_TRIANGULATION_AND_DOFHANDLER; const X::active_face_iterator &, const unsigned int, const Tensor<1, deal_II_space_dimension> &, - const FullMatrix &); + const FullMatrix &, + const double); template std::optional orthogonal_equality( @@ -263,7 +264,8 @@ for (X : SEQUENTIAL_TRIANGULATION_AND_DOFHANDLER; const X::face_iterator &, const unsigned int, const Tensor<1, deal_II_space_dimension> &, - const FullMatrix &); + const FullMatrix &, + const double); \} #endif } @@ -283,7 +285,8 @@ for (X : SEQUENTIAL_TRIANGULATION_AND_DOFHANDLERS; const unsigned int, std::vector> &, const Tensor<1, X::space_dimension> &, - const FullMatrix &); + const FullMatrix &, + const double); template void collect_periodic_faces( @@ -292,7 +295,8 @@ for (X : SEQUENTIAL_TRIANGULATION_AND_DOFHANDLERS; const unsigned int, std::vector> &, const Tensor<1, X::space_dimension> &, - const FullMatrix &); + const FullMatrix &, + const double); \} #endif @@ -323,7 +327,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) parallel::distributed::Triangulation< deal_II_dimension, deal_II_space_dimension>::space_dimension> &, - const FullMatrix &); + const FullMatrix &, + const double); template void collect_periodic_faces< @@ -340,7 +345,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) parallel::distributed::Triangulation< deal_II_dimension, deal_II_space_dimension>::space_dimension> &, - const FullMatrix &); + const FullMatrix &, + const double); \} # endif #endif @@ -369,7 +375,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) parallel::shared::Triangulation< deal_II_dimension, deal_II_space_dimension>::space_dimension> &, - const FullMatrix &); + const FullMatrix &, + const double); template void collect_periodic_faces< parallel::shared::Triangulation::space_dimension> &, - const FullMatrix &); + const FullMatrix &, + const double); \} #endif } diff --git a/tests/grid/grid_tools_collect_periodic_faces_tolerance.cc b/tests/grid/grid_tools_collect_periodic_faces_tolerance.cc new file mode 100644 index 0000000000..4d3a73dd51 --- /dev/null +++ b/tests/grid/grid_tools_collect_periodic_faces_tolerance.cc @@ -0,0 +1,90 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2023 by the deal.II authors +// +// This file is part of the deal.II library. +// +// Part of the source code is dual licensed under Apache-2.0 WITH +// LLVM-exception OR LGPL-2.1-or-later. Detailed license information +// governing the source code and code contributions can be found in +// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. +// +// ------------------------------------------------------------------------ +// +// Test a tolerance parameter to find the periodic nodes of a triangulation +// using the collect_periodic_faces function. +// Useful in case of CAD geometry that is not precise enough. The default +// value is at 1e-10 (absolute value). Sequential version. + +#include +#include +#include + +#include +#include +#include + +#include + +#include "../tests.h" + +template +void +test() +{ + deallog << "dim = " << dim << std::endl; + Triangulation tria; + + std::vector::cell_iterator>> + periodicity_vector; + + unsigned int num_refinements = 1 << 4; + Point p1; + Point p2; + std::vector repetitions(dim); + + for (unsigned int i = 0; i < dim; ++i) + { + p1[i] = 0.; + p2[i] = 1.; + repetitions[i] = num_refinements; + } + GridGenerator::subdivided_hyper_rectangle(tria, repetitions, p1, p2, true); + + // Here we are only interested in using the tolerance parameter of the + // collect_periodic_faces function, the other default parameter stay as + // default one. + const Tensor<1, dim> defaultOffset = dealii::Tensor<1, dim>(); + const FullMatrix matrix = FullMatrix(); + const double tolerance = 1e-7; + const unsigned int direction = 0; // x-direction + + // We distort the mesh to test the tolerance parameter + GridTools::distort_random(tolerance / 10, tria, /* keep boundary */ false); + + // Collect periodic faces in the x-direction + // activate debug to see it fail + GridTools::collect_periodic_faces(tria, + 0, // boundary Id 1 + 1, // boundary Id 2 + direction, + periodicity_vector, + defaultOffset, + matrix, + tolerance); + // Check the size of the matched_pairs vector + deallog << periodicity_vector.size() << std::endl; +} + +int +main(int argc, char *argv[]) +{ + initlog(); + + test<2>(); + test<3>(); + + return 0; +} diff --git a/tests/grid/grid_tools_collect_periodic_faces_tolerance.output b/tests/grid/grid_tools_collect_periodic_faces_tolerance.output new file mode 100644 index 0000000000..46c6716483 --- /dev/null +++ b/tests/grid/grid_tools_collect_periodic_faces_tolerance.output @@ -0,0 +1,5 @@ + +DEAL::dim = 2 +DEAL::16 +DEAL::dim = 3 +DEAL::256