]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add a tolerance parameter with a default value to the match of periodic nodes. 17921/head
authorYann Jobic <yjobic@gmail.com>
Tue, 10 Dec 2024 08:06:01 +0000 (09:06 +0100)
committerYann Jobic <yjobic@gmail.com>
Mon, 2 Jun 2025 12:04:17 +0000 (14:04 +0200)
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

include/deal.II/grid/grid_tools.h
source/grid/grid_tools_dof_handlers.cc
source/grid/grid_tools_dof_handlers.inst.in
tests/grid/grid_tools_collect_periodic_faces_tolerance.cc [new file with mode: 0644]
tests/grid/grid_tools_collect_periodic_faces_tolerance.output [new file with mode: 0644]

index df298799cad047619e61737f6a0c4d94c52f7c5e..79941f36eb1b6d76e3ee35404d90ad636e591617 100644 (file)
@@ -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 <tt>v_1</tt> and <tt>v_2</tt> 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 <tt>v_1</tt> and <tt>v_2</tt> 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<double> &matrix = FullMatrix<double>());
+    const FullMatrix<double> &matrix  = FullMatrix<double>(),
+    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<double> &matrix = FullMatrix<double>());
+    const FullMatrix<double> &matrix  = FullMatrix<double>(),
+    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<double> &matrix = FullMatrix<double>());
+    const FullMatrix<double> &matrix  = FullMatrix<double>(),
+    const double              abs_tol = 1e-10);
 
   /** @} */
   /**
index 3f7da7e914eb96831044688a1c4cea56318f7bd4..1a27a0e58842d08330948327a7612a6a7eed0f2e 100644 (file)
@@ -2127,7 +2127,8 @@ namespace GridTools
     std::vector<PeriodicFacePair<CellIterator>> &matched_pairs,
     const dealii::Tensor<1, CellIterator::AccessorType::space_dimension>
                              &offset,
-    const FullMatrix<double> &matrix)
+    const FullMatrix<double> &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<PeriodicFacePair<typename MeshType::cell_iterator>>
                                                &matched_pairs,
     const Tensor<1, MeshType::space_dimension> &offset,
-    const FullMatrix<double>                   &matrix)
+    const FullMatrix<double>                   &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<PeriodicFacePair<typename MeshType::cell_iterator>>
                                                &matched_pairs,
     const Tensor<1, MeshType::space_dimension> &offset,
-    const FullMatrix<double>                   &matrix)
+    const FullMatrix<double>                   &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<spacedim>     &point2,
                       const unsigned int         direction,
                       const Tensor<1, spacedim> &offset,
-                      const FullMatrix<double>  &matrix)
+                      const FullMatrix<double>  &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<double>                                     &matrix)
+    const FullMatrix<double>                                     &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;
index f03a2ac36780d13634a93de3cc88cd79a8175b8e..194530d489fc6999274edc18153e3a3cae3bd21e 100644 (file)
@@ -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<double> &);
+        const FullMatrix<double> &,
+        const double);
 
       template std::optional<types::geometric_orientation>
       orthogonal_equality<X::face_iterator>(
@@ -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<double> &);
+        const FullMatrix<double> &,
+        const double);
     \}
 #endif
   }
@@ -283,7 +285,8 @@ for (X : SEQUENTIAL_TRIANGULATION_AND_DOFHANDLERS;
         const unsigned int,
         std::vector<PeriodicFacePair<X::cell_iterator>> &,
         const Tensor<1, X::space_dimension> &,
-        const FullMatrix<double> &);
+        const FullMatrix<double> &,
+        const double);
 
       template void
       collect_periodic_faces<X>(
@@ -292,7 +295,8 @@ for (X : SEQUENTIAL_TRIANGULATION_AND_DOFHANDLERS;
         const unsigned int,
         std::vector<PeriodicFacePair<X::cell_iterator>> &,
         const Tensor<1, X::space_dimension> &,
-        const FullMatrix<double> &);
+        const FullMatrix<double> &,
+        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<double> &);
+        const FullMatrix<double> &,
+        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<double> &);
+        const FullMatrix<double> &,
+        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<double> &);
+        const FullMatrix<double> &,
+        const double);
       template void
       collect_periodic_faces<
         parallel::shared::Triangulation<deal_II_dimension,
@@ -385,7 +392,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<double> &);
+        const FullMatrix<double> &,
+        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 (file)
index 0000000..4d3a73d
--- /dev/null
@@ -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 <deal.II/base/exceptions.h>
+#include <deal.II/base/mpi.h>
+#include <deal.II/base/point.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/tria.h>
+
+#include <vector>
+
+#include "../tests.h"
+
+template <int dim>
+void
+test()
+{
+  deallog << "dim = " << dim << std::endl;
+  Triangulation<dim> tria;
+
+  std::vector<typename GridTools::PeriodicFacePair<
+    typename Triangulation<dim>::cell_iterator>>
+    periodicity_vector;
+
+  unsigned int              num_refinements = 1 << 4;
+  Point<dim>                p1;
+  Point<dim>                p2;
+  std::vector<unsigned int> 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<double> matrix        = FullMatrix<double>();
+  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 (file)
index 0000000..46c6716
--- /dev/null
@@ -0,0 +1,5 @@
+
+DEAL::dim = 2
+DEAL::16
+DEAL::dim = 3
+DEAL::256

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.