* 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
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
*
* 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
&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);
/**
&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);
/** @} */
/**
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;
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
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;
// 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())
{
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;
// 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);
}
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);
if (i == direction)
continue;
- if (std::abs(distance[i]) > 1.e-10)
+ if (std::abs(distance[i]) > abs_tol)
return false;
}
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"));
face2->vertex(*it),
direction,
offset,
- matrix))
+ matrix,
+ abs_tol))
{
face1_vertices[i] = *it;
face2_vertices[i] = i;
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>(
const X::face_iterator &,
const unsigned int,
const Tensor<1, deal_II_space_dimension> &,
- const FullMatrix<double> &);
+ const FullMatrix<double> &,
+ const double);
\}
#endif
}
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>(
const unsigned int,
std::vector<PeriodicFacePair<X::cell_iterator>> &,
const Tensor<1, X::space_dimension> &,
- const FullMatrix<double> &);
+ const FullMatrix<double> &,
+ const double);
\}
#endif
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<
parallel::distributed::Triangulation<
deal_II_dimension,
deal_II_space_dimension>::space_dimension> &,
- const FullMatrix<double> &);
+ const FullMatrix<double> &,
+ const double);
\}
# endif
#endif
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,
parallel::shared::Triangulation<
deal_II_dimension,
deal_II_space_dimension>::space_dimension> &,
- const FullMatrix<double> &);
+ const FullMatrix<double> &,
+ const double);
\}
#endif
}
--- /dev/null
+// ------------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+DEAL::dim = 2
+DEAL::16
+DEAL::dim = 3
+DEAL::256