From 1dc49f5f2b2c209a056a7eed6821838afe278333 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Schartum=20Dokken?= Date: Wed, 25 Dec 2024 21:28:55 +0100 Subject: [PATCH] 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 --- cpp/ContactConstraint.h | 8 +-- cpp/MultiPointConstraint.h | 2 + cpp/PeriodicConstraint.h | 75 +++++++++++++--------- cpp/utils.h | 8 ++- python/dolfinx_mpc/mpc.cpp | 13 ++-- python/dolfinx_mpc/multipointconstraint.py | 18 ++++-- 6 files changed, 77 insertions(+), 47 deletions(-) diff --git a/cpp/ContactConstraint.h b/cpp/ContactConstraint.h index 09e8b471..15d4102c 100644 --- a/cpp/ContactConstraint.h +++ b/cpp/ContactConstraint.h @@ -469,7 +469,7 @@ mpc_data create_contact_slip_condition( slave_coordinates, eps2); auto [basis, basis_shape] = dolfinx_mpc::evaluate_basis_functions( - 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> @@ -590,7 +590,7 @@ mpc_data create_contact_slip_condition( = dolfinx_mpc::find_local_collisions(*mesh, bb_tree, recv_coords, eps2); auto [recv_basis_values, shape] = dolfinx_mpc::evaluate_basis_functions( - V, recv_coords, remote_cell_collisions); + V, recv_coords, remote_cell_collisions, eps2); MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> basis_span(recv_basis_values.data(), shape[0], shape[1]); @@ -1000,7 +1000,7 @@ mpc_data create_contact_inelastic_condition( = dolfinx_mpc::find_local_collisions(*V.mesh(), bb_tree, slave_coordinates, eps2); auto [basis_values, basis_shape] = dolfinx_mpc::evaluate_basis_functions( - 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> @@ -1185,7 +1185,7 @@ mpc_data create_contact_inelastic_condition( = dolfinx_mpc::find_local_collisions(*V.mesh(), bb_tree, recv_coords, eps2); auto [basis, basis_shape] = dolfinx_mpc::evaluate_basis_functions( - 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> diff --git a/cpp/MultiPointConstraint.h b/cpp/MultiPointConstraint.h index 34606a2e..977cdcbc 100644 --- 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 index 5f6f9e9e..bc3daabb 100644 --- a/cpp/PeriodicConstraint.h +++ b/cpp/PeriodicConstraint.h @@ -20,6 +20,10 @@ 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 dolfinx_mpc::mpc_data _create_periodic_condition( @@ -27,8 +31,9 @@ dolfinx_mpc::mpc_data _create_periodic_condition( std::span slave_blocks, const std::function(std::span)>& relation, T scale, const std::function& parent_map, - const dolfinx::fem::FunctionSpace& parent_space) + const dolfinx::fem::FunctionSpace& 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& sub_dofs) { @@ -71,10 +76,6 @@ dolfinx_mpc::mpc_data _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::epsilon(); - auto mesh = V.mesh(); auto dofmap = V.dofmap(); auto imap = dofmap->index_map; @@ -138,7 +139,7 @@ dolfinx_mpc::mpc_data _create_periodic_condition( = dolfinx_mpc::find_local_collisions(*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( - 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> tabulated_basis_values(basis_values.data(), basis_shape); @@ -358,7 +359,7 @@ dolfinx_mpc::mpc_data _create_periodic_condition( = dolfinx_mpc::find_local_collisions(*mesh, tree, coords_recvb, tol); auto [remote_basis_valuesb, r_basis_shape] = dolfinx_mpc::evaluate_basis_functions(V, coords_recvb, - remote_cell_collisions); + remote_cell_collisions, tol); MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< const U, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents> remote_basis_values(remote_basis_valuesb.data(), r_basis_shape); @@ -490,6 +491,10 @@ dolfinx_mpc::mpc_data _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 dolfinx_mpc::mpc_data geometrical_condition( @@ -502,7 +507,7 @@ dolfinx_mpc::mpc_data geometrical_condition( indicator, const std::function(std::span)>& relation, const std::vector>>& bcs, - T scale, bool collapse) + T scale, bool collapse, const U tol) { std::vector reduced_blocks; if (collapse) @@ -526,7 +531,7 @@ dolfinx_mpc::mpc_data geometrical_condition( auto sub_map = [&parent_map](const std::int32_t& i) { return parent_map[i]; }; return _create_periodic_condition(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 geometrical_condition( reduced_blocks.push_back(slave_blocks[i]); auto sub_map = [](const std::int32_t& dof) { return dof; }; return _create_periodic_condition(*V, std::span(reduced_blocks), - relation, scale, sub_map, *V); + relation, scale, sub_map, *V, tol); } } @@ -558,6 +563,10 @@ dolfinx_mpc::mpc_data 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 dolfinx_mpc::mpc_data topological_condition( @@ -566,7 +575,7 @@ dolfinx_mpc::mpc_data topological_condition( const std::int32_t tag, const std::function(std::span)>& relation, const std::vector>>& bcs, - T scale, bool collapse) + T scale, bool collapse, const U tol) { std::vector entities = meshtag->find(tag); V->mesh()->topology_mutable()->create_connectivity( @@ -594,7 +603,7 @@ dolfinx_mpc::mpc_data topological_condition( = [&parent_map](const std::int32_t& i) { return parent_map[i]; }; // Create mpc on sub space dolfinx_mpc::mpc_data sub_data = _create_periodic_condition( - 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 topological_condition( const auto sub_map = [](const std::int32_t& dof) { return dof; }; return _create_periodic_condition(*V, std::span(reduced_blocks), - relation, scale, sub_map, *V); + relation, scale, sub_map, *V, tol); } }; @@ -633,10 +642,11 @@ mpc_data create_periodic_condition_geometrical( const std::function(std::span)>& relation, const std::vector>>& bcs, - double scale, bool collapse) + double scale, bool collapse, + const double tol = 500 * std::numeric_limits::epsilon()) { return impl::geometrical_condition(V, indicator, relation, - bcs, scale, collapse); + bcs, scale, collapse, tol); } mpc_data> create_periodic_condition_geometrical( @@ -651,10 +661,11 @@ mpc_data> create_periodic_condition_geometrical( const std::vector< std::shared_ptr>>>& bcs, - std::complex scale, bool collapse) + std::complex scale, bool collapse, + const double tol = 500 * std::numeric_limits::epsilon()) { return impl::geometrical_condition, double>( - V, indicator, relation, bcs, scale, collapse); + V, indicator, relation, bcs, scale, collapse, tol); } mpc_data create_periodic_condition_topological( @@ -664,10 +675,11 @@ mpc_data create_periodic_condition_topological( const std::function(std::span)>& relation, const std::vector>>& bcs, - double scale, bool collapse) + double scale, bool collapse, + const double tol = 500 * std::numeric_limits::epsilon()) { return impl::topological_condition(V, meshtag, tag, relation, - bcs, scale, collapse); + bcs, scale, collapse, tol); } mpc_data> create_periodic_condition_topological( @@ -678,10 +690,11 @@ mpc_data> create_periodic_condition_topological( const std::vector< std::shared_ptr>>>& bcs, - std::complex scale, bool collapse) + std::complex scale, bool collapse, + const double tol = 500 * std::numeric_limits::epsilon()) { return impl::topological_condition, double>( - V, meshtag, tag, relation, bcs, scale, collapse); + V, meshtag, tag, relation, bcs, scale, collapse, tol); } mpc_data create_periodic_condition_geometrical( @@ -695,10 +708,11 @@ mpc_data create_periodic_condition_geometrical( const std::function(std::span)>& relation, const std::vector>>& bcs, - float scale, bool collapse) + float scale, bool collapse, + const float tol = 500 * std::numeric_limits::epsilon()) { return impl::geometrical_condition(V, indicator, relation, bcs, - scale, collapse); + scale, collapse, tol); } mpc_data> create_periodic_condition_geometrical( @@ -713,10 +727,11 @@ mpc_data> create_periodic_condition_geometrical( const std::vector< std::shared_ptr>>>& bcs, - std::complex scale, bool collapse) + std::complex scale, bool collapse, + const float tol = 500 * std::numeric_limits::epsilon()) { return impl::geometrical_condition, float>( - V, indicator, relation, bcs, scale, collapse); + V, indicator, relation, bcs, scale, collapse, tol); } mpc_data create_periodic_condition_topological( @@ -726,10 +741,11 @@ mpc_data create_periodic_condition_topological( const std::function(std::span)>& relation, const std::vector>>& bcs, - float scale, bool collapse) + float scale, bool collapse, + const float tol = 500 * std::numeric_limits::epsilon()) { return impl::topological_condition(V, meshtag, tag, relation, - bcs, scale, collapse); + bcs, scale, collapse, tol); } mpc_data> create_periodic_condition_topological( @@ -740,10 +756,11 @@ mpc_data> create_periodic_condition_topological( const std::vector< std::shared_ptr>>>& bcs, - std::complex scale, bool collapse) + std::complex scale, bool collapse, + const float tol = 500 * std::numeric_limits::epsilon()) { return impl::topological_condition, float>( - V, meshtag, tag, relation, bcs, scale, collapse); + V, meshtag, tag, relation, bcs, scale, collapse, tol); } } // namespace dolfinx_mpc \ No newline at end of file diff --git a/cpp/utils.h b/cpp/utils.h index 3d9be9bb..a8aba51d 100644 --- a/cpp/utils.h +++ b/cpp/utils.h @@ -962,13 +962,16 @@ dolfinx_mpc::mpc_data 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::pair, std::array> evaluate_basis_functions(const dolfinx::fem::FunctionSpace& V, std::span x, - std::span cells) + std::span 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& V, else { // Pull-back physical point xp to reference coordinate Xp - cmap.pull_back_nonaffine(Xp, xp, coord_dofs, - 5000 * std::numeric_limits::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::compute_jacobian(dphi, coord_dofs, diff --git a/python/dolfinx_mpc/mpc.cpp b/python/dolfinx_mpc/mpc.cpp index 269bd48c..e5f70353 100644 --- a/python/dolfinx_mpc/mpc.cpp +++ b/python/dolfinx_mpc/mpc.cpp @@ -152,7 +152,7 @@ void declare_functions(nb::module_& m) nb::ndarray, nb::numpy>&)>& relation, const std::vector>>& 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>& V, @@ -193,7 +193,7 @@ void declare_functions(nb::module_& m) nb::ndarray, nb::numpy>&)>& relation, const std::vector>>& bcs, - T scale, bool collapse) + T scale, bool collapse, const U tol) { auto _relation = [&relation](std::span x) -> std::vector { @@ -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 diff --git a/python/dolfinx_mpc/multipointconstraint.py b/python/dolfinx_mpc/multipointconstraint.py index 8ef83eb6..45fa32a3 100644 --- a/python/dolfinx_mpc/multipointconstraint.py +++ b/python/dolfinx_mpc/multipointconstraint.py @@ -14,7 +14,7 @@ import dolfinx.mesh as _mesh import numpy import numpy.typing as npt -from dolfinx import default_scalar_type +from dolfinx import default_real_type, default_scalar_type import dolfinx_mpc.cpp @@ -230,6 +230,7 @@ def create_periodic_constraint_topological( relation: Callable[[numpy.ndarray], numpy.ndarray], bcs: List[_fem.DirichletBC], scale: _float_classes = default_scalar_type(1.0), + tol: _float_classes = 500 * numpy.finfo(default_real_type).eps, ): """ Create periodic condition for all closure dofs of on all entities in `meshtag` with value `tag`. @@ -242,17 +243,20 @@ def create_periodic_constraint_topological( relation: Lambda-function describing the geometrical relation bcs: Dirichlet boundary conditions for the problem (Periodic constraints will be ignored for these dofs) scale: Float for scaling bc + 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. """ bcs_ = [bc._cpp_object for bc in bcs] if isinstance(scale, numpy.generic): # nanobind conversion of numpy dtypes to general Python types scale = scale.item() # type: ignore if V is self.V: mpc_data = dolfinx_mpc.cpp.mpc.create_periodic_constraint_topological( - self.V._cpp_object, meshtag._cpp_object, tag, relation, bcs_, scale, False + self.V._cpp_object, meshtag._cpp_object, tag, relation, bcs_, scale, False, float(tol) ) elif self.V.contains(V): mpc_data = dolfinx_mpc.cpp.mpc.create_periodic_constraint_topological( - V._cpp_object, meshtag._cpp_object, tag, relation, bcs_, scale, True + V._cpp_object, meshtag._cpp_object, tag, relation, bcs_, scale, True, float(tol) ) else: raise RuntimeError("The input space has to be a sub space (or the full space) of the MPC") @@ -265,6 +269,7 @@ def create_periodic_constraint_geometrical( relation: Callable[[numpy.ndarray], numpy.ndarray], bcs: List[_fem.DirichletBC], scale: _float_classes = default_scalar_type(1.0), + tol: _float_classes = 500 * numpy.finfo(default_real_type).eps, ): """ Create a periodic condition for all degrees of freedom whose physical location satisfies @@ -278,17 +283,20 @@ def create_periodic_constraint_geometrical( bcs: Dirichlet boundary conditions for the problem (Periodic constraints will be ignored for these dofs) scale: Float for scaling bc + 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. """ if isinstance(scale, numpy.generic): # nanobind conversion of numpy dtypes to general Python types scale = scale.item() # type: ignore bcs = [] if bcs is None else [bc._cpp_object for bc in bcs] if V is self.V: mpc_data = dolfinx_mpc.cpp.mpc.create_periodic_constraint_geometrical( - self.V._cpp_object, indicator, relation, bcs, scale, False + self.V._cpp_object, indicator, relation, bcs, scale, False, float(tol) ) elif self.V.contains(V): mpc_data = dolfinx_mpc.cpp.mpc.create_periodic_constraint_geometrical( - V._cpp_object, indicator, relation, bcs, scale, True + V._cpp_object, indicator, relation, bcs, scale, True, float(tol) ) else: raise RuntimeError("The input space has to be a sub space (or the full space) of the MPC")