Move tolerances out of hardcoded part of code (#140)

* Move tolerances out of hardcoded part of code * Remove debug print * Remove default argument in evaluate_basis_functions * real type of tolerance * Fix typing
jorgensd · Dec 25, 2024 · 1dc49f5 · 1dc49f5
1 parent 9ea0364
commit 1dc49f5
Show file tree

Hide file tree

Showing 6 changed files with 77 additions and 47 deletions.
diff --git a/cpp/ContactConstraint.h b/cpp/ContactConstraint.h
@@ -469,7 +469,7 @@ mpc_data<T> create_contact_slip_condition(
                                                 slave_coordinates, eps2);
 
     auto [basis, basis_shape] = dolfinx_mpc::evaluate_basis_functions<U>(
-        V, slave_coordinates, local_cell_collisions);
+        V, slave_coordinates, local_cell_collisions, eps2);
     assert(basis_shape.back() == 1);
     MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
         U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>
@@ -590,7 +590,7 @@ mpc_data<T> create_contact_slip_condition(
         = dolfinx_mpc::find_local_collisions<U>(*mesh, bb_tree, recv_coords,
                                                 eps2);
     auto [recv_basis_values, shape] = dolfinx_mpc::evaluate_basis_functions<U>(
-        V, recv_coords, remote_cell_collisions);
+        V, recv_coords, remote_cell_collisions, eps2);
     MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
         U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>
         basis_span(recv_basis_values.data(), shape[0], shape[1]);
@@ -1000,7 +1000,7 @@ mpc_data<T> create_contact_inelastic_condition(
         = dolfinx_mpc::find_local_collisions<U>(*V.mesh(), bb_tree,
                                                 slave_coordinates, eps2);
     auto [basis_values, basis_shape] = dolfinx_mpc::evaluate_basis_functions<U>(
-        V, slave_coordinates, colliding_cells);
+        V, slave_coordinates, colliding_cells, eps2);
     assert(basis_shape.back() == 1);
     MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
         U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>
@@ -1185,7 +1185,7 @@ mpc_data<T> create_contact_inelastic_condition(
         = dolfinx_mpc::find_local_collisions<U>(*V.mesh(), bb_tree, recv_coords,
                                                 eps2);
     auto [basis, basis_shape] = dolfinx_mpc::evaluate_basis_functions<U>(
-        V, recv_coords, remote_cell_collisions);
+        V, recv_coords, remote_cell_collisions, eps2);
     assert(basis_shape.back() == 1);
     MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
         const U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>

diff --git a/cpp/MultiPointConstraint.h b/cpp/MultiPointConstraint.h
@@ -134,8 +134,10 @@ class MultiPointConstraint
       auto coeffs = _coeff_map->links(slave);
       assert(masters.size() == coeffs.size());
       for (std::size_t k = 0; k < masters.size(); ++k)
+      {
         vector[slave]
             += coeffs[k] * vector[masters[k]]; //+ _mpc_constants[slave];
+      }
     }
   };
 

diff --git a/cpp/PeriodicConstraint.h b/cpp/PeriodicConstraint.h
@@ -20,15 +20,20 @@ namespace impl
 /// not collapsed)
 /// @param[in] parent_space The parent space (The same space as V if not
 /// collapsed)
+/// @param[in] tol Tolerance for adding scaled basis values to MPC. Any
+/// contribution that is less than this value is ignored. The tolerance is also
+/// added as padding for the bounding box trees and corresponding collision
+/// searches to determine periodic degrees of freedom.
 /// @returns The multi point constraint
 template <typename T, std::floating_point U>
 dolfinx_mpc::mpc_data<T> _create_periodic_condition(
     const dolfinx::fem::FunctionSpace<U>& V,
     std::span<std::int32_t> slave_blocks,
     const std::function<std::vector<U>(std::span<const U>)>& relation, T scale,
     const std::function<const std::int32_t(const std::int32_t&)>& parent_map,
-    const dolfinx::fem::FunctionSpace<U>& parent_space)
+    const dolfinx::fem::FunctionSpace<U>& parent_space, const U tol)
 {
+
   // Map a list of indices in collapsed space back to the parent space
   auto sub_to_parent = [&parent_map](const std::vector<std::int32_t>& sub_dofs)
   {
@@ -71,10 +76,6 @@ dolfinx_mpc::mpc_data<T> _create_periodic_condition(
         "Periodic conditions for vector valued spaces are not "
         "implemented");
 
-  // Tolerance for adding scaled basis values to MPC. Any scaled basis
-  // value with lower absolute value than the tolerance is ignored
-  const U tol = 500 * std::numeric_limits<U>::epsilon();
-
   auto mesh = V.mesh();
   auto dofmap = V.dofmap();
   auto imap = dofmap->index_map;
@@ -138,7 +139,7 @@ dolfinx_mpc::mpc_data<T> _create_periodic_condition(
       = dolfinx_mpc::find_local_collisions<U>(*mesh, tree, mapped_T_b, tol);
   dolfinx::common::Timer t0("~~Periodic: Local cell and eval basis");
   auto [basis_values, basis_shape] = dolfinx_mpc::evaluate_basis_functions<U>(
-      V, mapped_T_b, local_cell_collisions);
+      V, mapped_T_b, local_cell_collisions, tol);
   MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
       const U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 3>>
       tabulated_basis_values(basis_values.data(), basis_shape);
@@ -358,7 +359,7 @@ dolfinx_mpc::mpc_data<T> _create_periodic_condition(
       = dolfinx_mpc::find_local_collisions<U>(*mesh, tree, coords_recvb, tol);
   auto [remote_basis_valuesb, r_basis_shape]
       = dolfinx_mpc::evaluate_basis_functions<U>(V, coords_recvb,
-                                                 remote_cell_collisions);
+                                                 remote_cell_collisions, tol);
   MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
       const U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 3>>
       remote_basis_values(remote_basis_valuesb.data(), r_basis_shape);
@@ -490,6 +491,10 @@ dolfinx_mpc::mpc_data<T> _create_periodic_condition(
 /// @param[in] scale Scaling of the periodic condition
 /// @param[in] collapse If true, the list of marked dofs is in the collapsed
 /// input space
+/// @param[in] tol Tolerance for adding scaled basis values to MPC. Any
+/// contribution that is less than this value is ignored. The tolerance is also
+/// added as padding for the bounding box trees and corresponding collision
+/// searches to determine periodic degrees of freedom.
 /// @returns The multi point constraint
 template <typename T, std::floating_point U>
 dolfinx_mpc::mpc_data<T> geometrical_condition(
@@ -502,7 +507,7 @@ dolfinx_mpc::mpc_data<T> geometrical_condition(
         indicator,
     const std::function<std::vector<U>(std::span<const U>)>& relation,
     const std::vector<std::shared_ptr<const dolfinx::fem::DirichletBC<T>>>& bcs,
-    T scale, bool collapse)
+    T scale, bool collapse, const U tol)
 {
   std::vector<std::int32_t> reduced_blocks;
   if (collapse)
@@ -526,7 +531,7 @@ dolfinx_mpc::mpc_data<T> geometrical_condition(
     auto sub_map
         = [&parent_map](const std::int32_t& i) { return parent_map[i]; };
     return _create_periodic_condition<T>(V_sub, std::span(reduced_blocks),
-                                         relation, scale, sub_map, *V);
+                                         relation, scale, sub_map, *V, tol);
   }
   else
   {
@@ -542,7 +547,7 @@ dolfinx_mpc::mpc_data<T> geometrical_condition(
         reduced_blocks.push_back(slave_blocks[i]);
     auto sub_map = [](const std::int32_t& dof) { return dof; };
     return _create_periodic_condition<T>(*V, std::span(reduced_blocks),
-                                         relation, scale, sub_map, *V);
+                                         relation, scale, sub_map, *V, tol);
   }
 }
 
@@ -558,6 +563,10 @@ dolfinx_mpc::mpc_data<T> geometrical_condition(
 /// @param[in] scale Scaling of the periodic condition
 /// @param[in] collapse If true, the list of marked dofs is in the collapsed
 /// input space
+/// @param[in] tol Tolerance for adding scaled basis values to MPC. Any
+/// contribution that is less than this value is ignored. The tolerance is also
+/// added as padding for the bounding box trees and corresponding collision
+/// searches to determine periodic degrees of freedom.
 /// @returns The multi point constraint
 template <typename T, std::floating_point U>
 dolfinx_mpc::mpc_data<T> topological_condition(
@@ -566,7 +575,7 @@ dolfinx_mpc::mpc_data<T> topological_condition(
     const std::int32_t tag,
     const std::function<std::vector<U>(std::span<const U>)>& relation,
     const std::vector<std::shared_ptr<const dolfinx::fem::DirichletBC<T>>>& bcs,
-    T scale, bool collapse)
+    T scale, bool collapse, const U tol)
 {
   std::vector<std::int32_t> entities = meshtag->find(tag);
   V->mesh()->topology_mutable()->create_connectivity(
@@ -594,7 +603,7 @@ dolfinx_mpc::mpc_data<T> topological_condition(
         = [&parent_map](const std::int32_t& i) { return parent_map[i]; };
     // Create mpc on sub space
     dolfinx_mpc::mpc_data<T> sub_data = _create_periodic_condition<T>(
-        V_sub, std::span(reduced_blocks), relation, scale, sub_map, *V);
+        V_sub, std::span(reduced_blocks), relation, scale, sub_map, *V, tol);
     return sub_data;
   }
   else
@@ -613,7 +622,7 @@ dolfinx_mpc::mpc_data<T> topological_condition(
     const auto sub_map = [](const std::int32_t& dof) { return dof; };
 
     return _create_periodic_condition<T, U>(*V, std::span(reduced_blocks),
-                                            relation, scale, sub_map, *V);
+                                            relation, scale, sub_map, *V, tol);
   }
 };
 
@@ -633,10 +642,11 @@ mpc_data<double> create_periodic_condition_geometrical(
     const std::function<std::vector<double>(std::span<const double>)>& relation,
     const std::vector<std::shared_ptr<const dolfinx::fem::DirichletBC<double>>>&
         bcs,
-    double scale, bool collapse)
+    double scale, bool collapse,
+    const double tol = 500 * std::numeric_limits<double>::epsilon())
 {
   return impl::geometrical_condition<double, double>(V, indicator, relation,
-                                                     bcs, scale, collapse);
+                                                     bcs, scale, collapse, tol);
 }
 
 mpc_data<std::complex<double>> create_periodic_condition_geometrical(
@@ -651,10 +661,11 @@ mpc_data<std::complex<double>> create_periodic_condition_geometrical(
     const std::vector<
         std::shared_ptr<const dolfinx::fem::DirichletBC<std::complex<double>>>>&
         bcs,
-    std::complex<double> scale, bool collapse)
+    std::complex<double> scale, bool collapse,
+    const double tol = 500 * std::numeric_limits<double>::epsilon())
 {
   return impl::geometrical_condition<std::complex<double>, double>(
-      V, indicator, relation, bcs, scale, collapse);
+      V, indicator, relation, bcs, scale, collapse, tol);
 }
 
 mpc_data<double> create_periodic_condition_topological(
@@ -664,10 +675,11 @@ mpc_data<double> create_periodic_condition_topological(
     const std::function<std::vector<double>(std::span<const double>)>& relation,
     const std::vector<std::shared_ptr<const dolfinx::fem::DirichletBC<double>>>&
         bcs,
-    double scale, bool collapse)
+    double scale, bool collapse,
+    const double tol = 500 * std::numeric_limits<double>::epsilon())
 {
   return impl::topological_condition<double, double>(V, meshtag, tag, relation,
-                                                     bcs, scale, collapse);
+                                                     bcs, scale, collapse, tol);
 }
 
 mpc_data<std::complex<double>> create_periodic_condition_topological(
@@ -678,10 +690,11 @@ mpc_data<std::complex<double>> create_periodic_condition_topological(
     const std::vector<
         std::shared_ptr<const dolfinx::fem::DirichletBC<std::complex<double>>>>&
         bcs,
-    std::complex<double> scale, bool collapse)
+    std::complex<double> scale, bool collapse,
+    const double tol = 500 * std::numeric_limits<double>::epsilon())
 {
   return impl::topological_condition<std::complex<double>, double>(
-      V, meshtag, tag, relation, bcs, scale, collapse);
+      V, meshtag, tag, relation, bcs, scale, collapse, tol);
 }
 
 mpc_data<float> create_periodic_condition_geometrical(
@@ -695,10 +708,11 @@ mpc_data<float> create_periodic_condition_geometrical(
     const std::function<std::vector<float>(std::span<const float>)>& relation,
     const std::vector<std::shared_ptr<const dolfinx::fem::DirichletBC<float>>>&
         bcs,
-    float scale, bool collapse)
+    float scale, bool collapse,
+    const float tol = 500 * std::numeric_limits<float>::epsilon())
 {
   return impl::geometrical_condition<float, float>(V, indicator, relation, bcs,
-                                                   scale, collapse);
+                                                   scale, collapse, tol);
 }
 
 mpc_data<std::complex<float>> create_periodic_condition_geometrical(
@@ -713,10 +727,11 @@ mpc_data<std::complex<float>> create_periodic_condition_geometrical(
     const std::vector<
         std::shared_ptr<const dolfinx::fem::DirichletBC<std::complex<float>>>>&
         bcs,
-    std::complex<float> scale, bool collapse)
+    std::complex<float> scale, bool collapse,
+    const float tol = 500 * std::numeric_limits<float>::epsilon())
 {
   return impl::geometrical_condition<std::complex<float>, float>(
-      V, indicator, relation, bcs, scale, collapse);
+      V, indicator, relation, bcs, scale, collapse, tol);
 }
 
 mpc_data<float> create_periodic_condition_topological(
@@ -726,10 +741,11 @@ mpc_data<float> create_periodic_condition_topological(
     const std::function<std::vector<float>(std::span<const float>)>& relation,
     const std::vector<std::shared_ptr<const dolfinx::fem::DirichletBC<float>>>&
         bcs,
-    float scale, bool collapse)
+    float scale, bool collapse,
+    const float tol = 500 * std::numeric_limits<float>::epsilon())
 {
   return impl::topological_condition<float, float>(V, meshtag, tag, relation,
-                                                   bcs, scale, collapse);
+                                                   bcs, scale, collapse, tol);
 }
 
 mpc_data<std::complex<float>> create_periodic_condition_topological(
@@ -740,10 +756,11 @@ mpc_data<std::complex<float>> create_periodic_condition_topological(
     const std::vector<
         std::shared_ptr<const dolfinx::fem::DirichletBC<std::complex<float>>>>&
         bcs,
-    std::complex<float> scale, bool collapse)
+    std::complex<float> scale, bool collapse,
+    const float tol = 500 * std::numeric_limits<float>::epsilon())
 {
   return impl::topological_condition<std::complex<float>, float>(
-      V, meshtag, tag, relation, bcs, scale, collapse);
+      V, meshtag, tag, relation, bcs, scale, collapse, tol);
 }
 
 } // namespace dolfinx_mpc
diff --git a/cpp/utils.h b/cpp/utils.h
@@ -962,13 +962,16 @@ dolfinx_mpc::mpc_data<T> distribute_ghost_data(
 /// @param[in,out] u The values at the points. Values are not computed
 /// for points with a negative cell index. This argument must be
 /// passed with the correct size.
+/// @param[in] tol The stopping tolerance for the l2 norm of the Newton update
+/// when pulling back a point on a non-affine cell. For affine cells this value
+/// is not used.
 /// @returns basis values (not unrolled for block size) for each point. shape
 /// (num_points, number_of_dofs, value_size). Flattened row major
 template <std::floating_point U>
 std::pair<std::vector<U>, std::array<std::size_t, 3>>
 evaluate_basis_functions(const dolfinx::fem::FunctionSpace<U>& V,
                          std::span<const U> x,
-                         std::span<const std::int32_t> cells)
+                         std::span<const std::int32_t> cells, const U tol)
 {
   assert(x.size() % 3 == 0);
   const std::size_t num_points = x.size() / 3;
@@ -1132,8 +1135,7 @@ evaluate_basis_functions(const dolfinx::fem::FunctionSpace<U>& V,
     else
     {
       // Pull-back physical point xp to reference coordinate Xp
-      cmap.pull_back_nonaffine(Xp, xp, coord_dofs,
-                               5000 * std::numeric_limits<U>::epsilon(), 15);
+      cmap.pull_back_nonaffine(Xp, xp, coord_dofs, tol, 15);
 
       cmap.tabulate(1, std::span(Xpb.data(), tdim), {1, tdim}, phi_b);
       dolfinx::fem::CoordinateElement<U>::compute_jacobian(dphi, coord_dofs,

diff --git a/python/dolfinx_mpc/mpc.cpp b/python/dolfinx_mpc/mpc.cpp
@@ -152,7 +152,7 @@ void declare_functions(nb::module_& m)
              nb::ndarray<const U, nb::ndim<2>, nb::numpy>&)>& relation,
          const std::vector<std::shared_ptr<const dolfinx::fem::DirichletBC<T>>>&
              bcs,
-         T scale, bool collapse)
+         T scale, bool collapse, const U tol)
       {
         auto _indicator
             = [&indicator](MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
@@ -180,10 +180,10 @@ void declare_functions(nb::module_& m)
           return output;
         };
         return dolfinx_mpc::create_periodic_condition_geometrical(
-            V, _indicator, _relation, bcs, scale, collapse);
+            V, _indicator, _relation, bcs, scale, collapse, tol);
       },
       "V"_a, "indicator"_a, "relation"_a, "bcs"_a, nb::arg("scale").noconvert(),
-      nb::arg("collapse").noconvert());
+      nb::arg("collapse").noconvert(), nb::arg("tol").noconvert());
   m.def(
       "create_periodic_constraint_topological",
       [](std::shared_ptr<const dolfinx::fem::FunctionSpace<U>>& V,
@@ -193,7 +193,7 @@ void declare_functions(nb::module_& m)
              nb::ndarray<const U, nb::ndim<2>, nb::numpy>&)>& relation,
          const std::vector<std::shared_ptr<const dolfinx::fem::DirichletBC<T>>>&
              bcs,
-         T scale, bool collapse)
+         T scale, bool collapse, const U tol)
       {
         auto _relation = [&relation](std::span<const U> x) -> std::vector<U>
         {
@@ -204,10 +204,11 @@ void declare_functions(nb::module_& m)
           return output;
         };
         return dolfinx_mpc::create_periodic_condition_topological(
-            V, meshtags, dim, _relation, bcs, scale, collapse);
+            V, meshtags, dim, _relation, bcs, scale, collapse, tol);
       },
       "V"_a, "meshtags"_a, "dim"_a, "relation"_a, "bcs"_a,
-      nb::arg("scale").noconvert(), nb::arg("collapse").noconvert());
+      nb::arg("scale").noconvert(), nb::arg("collapse").noconvert(),
+      nb::arg("tol").noconvert());
 }
 
 template <typename T, std::floating_point U>