diff --git a/mesh_handle/adapt.hxx b/mesh_handle/adapt.hxx new file mode 100644 index 0000000000..b23f689624 --- /dev/null +++ b/mesh_handle/adapt.hxx @@ -0,0 +1,283 @@ +/* + This file is part of t8code. + t8code is a C library to manage a collection (a forest) of multiple + connected adaptive space-trees of general element classes in parallel. + + Copyright (C) 2025 the developers + + t8code is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + t8code is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with t8code; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +*/ + +/** \file adapt.hxx + * This file provides the functionality to adapt a \ref t8_mesh_handle::mesh + * according to user defined callbacks. + */ + +#ifndef T8_ADAPT_HXX +#define T8_ADAPT_HXX + +#include +#include +#include "mesh.hxx" +#include +#include + +namespace t8_mesh_handle +{ +/** Callback function prototype to decide for coarsening of a family of elements in a mesh handle. + * This callback works with a family of elements, i.e., multiple elements that can be coarsened + * into one parent element. + * \param [in] mesh The mesh that should be adapted. + * \param [in] elements The element family considered to be coarsened. + * \tparam TMeshClass The mesh handle class. + * \return true if the family should be coarsened, false otherwise. + */ +template +using coarsen_element_family + = std::function& elements)>; + +/** Callback function prototype to decide for the refinement of an element of a mesh handle. + * \param [in] mesh The mesh that should be adapted. + * \param [in] element The element to consider for refinement. + * \tparam TMeshClass The mesh handle class. + * \return true if the element should be refined, false otherwise. + */ +template +using refine_element = std::function; + +/** Namespace detail to hide implementation details from the user. */ +namespace detail +{ + +/** Virtual base class for mesh adaptation contexts. + * We need this base class and not only \ref MeshAdaptContext for the \ref AdaptRegistry. + * AdaptRegistry should not be templated because we need to access registered contexts in \ref mesh_adapt_callback_wrapper, + * where we do not know the type of the mesh. Therefore, we work with a map of forests to instances of this base class to remain template free. + */ +struct MeshAdaptContextBase +{ + /** Virtual destructor for safe polymorphic deletion. + */ + virtual ~MeshAdaptContextBase () = default; + + /** Pure virtual callback for mesh adaptation. + * \param[in] lelement_handle_id Local element ID in the mesh handle. + * \param [in] is_family If 1, the entries in \a elements form a family. If 0, they do not. + * \param [in] num_elements The number of entries in \a elements. + * \param [in] elements Pointers to a family or, if \a is_family is zero, pointer to one element. + * \return 1 if the first entry in \a elements should be refined, + * -1 if the family \a elements shall be coarsened, + * 0 else. + */ + virtual int + adapt_callback (const t8_locidx_t lelement_handle_id, const int is_family, const int num_elements, + t8_element_t* elements[]) + = 0; +}; + +/** Mesh adaptation context holding the mesh handle and the user defined callbacks. + * Class inherits from \ref MeshAdaptContextBase and implements the virtual adapt callback using the mesh and the callbacks. + * \tparam TMesh The mesh handle class. + */ +template +struct MeshAdaptContext final: MeshAdaptContextBase +{ + using Element = typename TMesh::element_class; /**< Type alias for the element class. */ + + /** Constructor of the context with the mesh handle and the user defined callbacks. + * \param [in] mesh_handle The mesh handle to adapt. + * \param [in] refine_callback The refinement callback. + * \param [in] coarsen_callback The coarsening callback. + */ + MeshAdaptContext (TMesh& mesh_handle, refine_element refine_callback, + coarsen_element_family coarsen_callback) + : m_mesh_handle (mesh_handle), m_refine_callback (std::move (refine_callback)), + m_coarsen_callback (std::move (coarsen_callback)) + { + } + + /** Callback for mesh adaptation using user defined callbacks. + * \param [in] lelement_handle_id Local flat element ID in the mesh handle. + * \param [in] is_family If 1, the entries in \a elements form a family. If 0, they do not. + * \param [in] num_elements The number of entries in \a elements. + * \param [in] elements Pointers to a family or, if \a is_family is zero, pointer to one element. + * \return 1 if the first entry in \a elements should be refined, + * -1 if the family \a elements shall be coarsened, + * 0 else. + */ + int + adapt_callback (const t8_locidx_t lelement_handle_id, const int is_family, const int num_elements, + t8_element_t* elements[]) override + { + // Check if refine callback is set and call it using the correct mesh handle function arguments. + if (m_refine_callback) { + Element elem = m_mesh_handle[lelement_handle_id]; + if (m_refine_callback (m_mesh_handle, elem)) { + return 1; + } + } + + // Check if a family is provided for adaption, if coarsen callback is set and call it using the correct mesh handle function arguments. + if (is_family && m_coarsen_callback) { + std::vector element_family; + for (int i = 0; i < num_elements; i++) { + element_family.push_back (m_mesh_handle[lelement_handle_id + i]); + } + if (m_coarsen_callback (m_mesh_handle, element_family)) { + return -1; + } + } + // Do nothing if the callbacks no refinement or coarsening should be done according to the callbacks. + return 0; + } + + private: + TMesh& m_mesh_handle; /**< The mesh handle to adapt. */ + refine_element m_refine_callback; /**< The refinement callback. */ + coarsen_element_family m_coarsen_callback; /**< The coarsening callback. */ +}; + +/** Registry pattern is used to register contexts, which provides access to the callbacks and the mesh handle. + * This globally accessible static class is required to get the handle and the callbacks in the forest callback, + * as the predefined header permits to give these as function arguments. + */ +class AdaptRegistry { + public: + /** Static function to register \a context using \a forest as identifier. + * This makes the context publicly available using the Registry. + * \param [in] forest The forest identifier. + * \param [in] context The context to register. + */ + static void + register_context (t8_forest_t forest, MeshAdaptContextBase& context) + { + auto& map = get_map (); + auto [it, inserted] = map.emplace (forest, context); + if (!inserted) { + throw std::logic_error ("Context already registered"); + } + } + + /** Static function to unregister a context using \a forest as identifier. + * \param [in] forest The forest identifier. + */ + static void + unregister_context (t8_forest_t forest) + { + auto& map = get_map (); + const auto erased = map.erase (forest); + T8_ASSERT (erased == 1); + } + + /** Getter for a context using \a forest as identifier. + * \param [in] forest The forest identifier. + * \return Pointer to the context registered with the id \a forest if found, nullptr otherwise. + */ + static MeshAdaptContextBase* + get (t8_forest_t forest) + { + auto& map = get_map (); + auto it = map.find (forest); + return it != map.end () ? &it->second.get () : nullptr; + } + + private: + /** Get the static map associating t8_forest_t with MeshAdaptContextBase references. + * We use a getter instead of private member variable to ensure single initialization + * \return Reference to the static unordered map of t8_forest_t to MeshAdaptContextBase references. + */ + static std::unordered_map>& + get_map () + { + static std::unordered_map> map; + return map; + } +}; + +/** Wrapper around the mesh handle adapt functionality to be able to pass the callbacks to the classic adapt routine of a forest. + * The function header fits the definition of \ref t8_forest_adapt_t. + * \param [in] forest Unused; forest to which the new elements belong. + * \param [in] forest_from Forest that is adapted. + * \param [in] which_tree Local tree containing \a elements. + * \param [in] tree_class Unused; eclass of \a which_tree. + * \param [in] lelement_id The local element id in the tree of the first element in elements. + * \param [in] scheme Unused; scheme of the forest. + * \param [in] is_family If 1, the entries in \a elements form a family. If 0, they do not. + * \param [in] num_elements The number of entries in \a elements. + * \param [in] elements Pointers to a family or, if \a is_family is zero, pointer to one element. + * \return 1 if the first entry in \a elements should be refined, + * -1 if the family \a elements shall be coarsened, + * 0 else. + */ +int +mesh_adapt_callback_wrapper ([[maybe_unused]] t8_forest_t forest, t8_forest_t forest_from, t8_locidx_t which_tree, + [[maybe_unused]] t8_eclass_t tree_class, t8_locidx_t lelement_id, + [[maybe_unused]] const t8_scheme* scheme, const int is_family, const int num_elements, + t8_element_t* elements[]) +{ + // Get static adapt context from the registry. + // Via this, we can access the mesh handle and the user defined callbacks that are using mesh handle functionality. + auto* context = AdaptRegistry::get (forest_from); + if (!context) { + t8_global_infof ( + "Something went wrong while registering the adaption callbacks. Please check your implementation."); + return 0; // No adaption as default. + } + // Convert to index used in the mesh handle. + t8_locidx_t mesh_index = t8_forest_get_tree_element_offset (forest_from, which_tree) + lelement_id; + // Call the actual adapt callback stored in the context. + return context->adapt_callback (mesh_index, is_family, num_elements, elements); +} + +} // namespace detail + +/** Adapt a mesh handle according to the provided callbacks. + * The forest used to define the mesh handle is replaced in this function. + * If you want to additionally keep the current forest, call \ref t8_forest_ref before this function. + * \param [in,out] mesh_handle The mesh handle to adapt. + * \param [in] refine_callback Callback to decide if a single element should be refined. + * \param [in] coarsen_callback Callback to decide if a family of elements should be coarsened. + * \param [in] recursive Specifying whether adaptation is to be done recursively or not. + * \tparam TMesh The mesh handle class. + */ +template +void +adapt_mesh (TMesh& mesh_handle, refine_element refine_callback, coarsen_element_family coarsen_callback, + bool recursive) +{ + auto forest_from = mesh_handle.get_forest (); + // Initialize forest for the adapted mesh. + t8_forest_t forest; + t8_forest_init (&forest); + + // Create and register adaptation context holding the mesh handle and the user defined callbacks. + auto context + = detail::MeshAdaptContext (mesh_handle, std::move (refine_callback), std::move (coarsen_callback)); + detail::AdaptRegistry::register_context (forest_from, context); + + // Set up the forest for adaptation using the wrapper callback. + t8_forest_set_adapt (forest, forest_from, detail::mesh_adapt_callback_wrapper, recursive); + t8_forest_set_ghost (forest, 1, T8_GHOST_FACES); + t8_forest_set_user_data (forest, t8_forest_get_user_data (forest_from)); + t8_forest_commit (forest); + + // Replace the forest in the mesh handle with the adapted forest. + mesh_handle.set_forest (forest); + // Clean up. + detail::AdaptRegistry::unregister_context (forest_from); +} + +} // namespace t8_mesh_handle +#endif /* !T8_ADAPT_HXX */ diff --git a/mesh_handle/mesh.hxx b/mesh_handle/mesh.hxx index 2dd8f52ed5..e30459121b 100644 --- a/mesh_handle/mesh.hxx +++ b/mesh_handle/mesh.hxx @@ -73,6 +73,16 @@ class mesh { update_elements (); } + /** + * Destructor for a mesh of the handle. + * The forest in use will be unreferenced. + * Call \ref t8_forest_ref before if you want to keep it alive. + */ + ~mesh () + { + t8_forest_unref (&m_forest); + } + /** * Getter for the number of local elements in the mesh. * \return Number of local elements in the mesh. diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 3a7ad0f181..0ca7623b19 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -217,6 +217,7 @@ if( T8CODE_BUILD_MESH_HANDLE ) add_t8_cpp_test( NAME t8_gtest_custom_competence_serial SOURCES mesh_handle/t8_gtest_custom_competence.cxx ) add_t8_cpp_test( NAME t8_gtest_cache_competence_serial SOURCES mesh_handle/t8_gtest_cache_competence.cxx ) add_t8_cpp_test( NAME t8_gtest_handle_data_parallel SOURCES mesh_handle/t8_gtest_handle_data.cxx ) + add_t8_cpp_test( NAME t8_gtest_handle_adapt_serial SOURCES mesh_handle/t8_gtest_adapt.cxx ) endif() copy_test_file( test_cube_unstructured_1.inp ) diff --git a/test/mesh_handle/t8_gtest_adapt.cxx b/test/mesh_handle/t8_gtest_adapt.cxx new file mode 100644 index 0000000000..344415576c --- /dev/null +++ b/test/mesh_handle/t8_gtest_adapt.cxx @@ -0,0 +1,142 @@ +/* +This file is part of t8code. +t8code is a C library to manage a collection (a forest) of multiple +connected adaptive space-trees of general element classes in parallel. + +Copyright (C) 2025 the developers + +t8code is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +t8code is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with t8code; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +*/ + +/** + * \file t8_gtest_adapt.cxx + * Tests for the adapt routines of mesh handles. + * This tests uses the callbacks and user data of tutorial step 3 as example. + * The adaptation criterion is to look at the midpoint coordinates of the current element and if + * they are inside a sphere around a given midpoint we refine, if they are outside, we coarsen. + */ + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +/** Dummy user data taken from tutorial for test purposes. */ +struct dummy_user_data +{ + t8_3D_point midpoint; /**< The midpoint of our sphere. */ + double refine_if_inside_radius; /**< if an element's center is smaller than this value, we refine the element. */ + double coarsen_if_outside_radius; /**< if an element's center is larger this value, we coarsen its family. */ +}; + +/** Callback function implementation to decide for the refinement of an element of a mesh handle. + * \param [in] mesh The mesh that should be adapted. + * \param [in] element The element to consider for refinement. + * \tparam TMeshClass The mesh handle class. + * \return true if the element should be refined, false otherwise. + */ +template +bool +refine_element_test (const TMeshClass &mesh, const typename TMeshClass::element_class &element) +{ + typename TMeshClass::UserDataType user_data = mesh.get_user_data (); + auto element_centroid = element.get_centroid (); + double dist = t8_dist (element_centroid, user_data.midpoint); + return (dist < user_data.refine_if_inside_radius); +} + +/** Callback function implementation to decide for coarsening. + * \param [in] mesh The mesh that should be adapted. + * \param [in] elements The element family considered to be coarsened. + * \tparam TMeshClass The mesh handle class. + * \return true if the family should be coarsened, false otherwise. + */ +template +bool +coarsen_element_family_test (const TMeshClass &mesh, const std::vector &elements) +{ + typename TMeshClass::UserDataType user_data = mesh.get_user_data (); + auto element_centroid = elements[0].get_centroid (); + double dist = t8_dist (element_centroid, user_data.midpoint); + return (dist > user_data.coarsen_if_outside_radius); +} + +/** Adapt callback implementation for a forest. + * This callback defines the same adaptation rules as \ref coarsen_element_family_test + * together with \ref refine_element_test. Callback is not used for the mesh handle + * but for the forest compared to the adapted mesh handle. + */ +int +forest_adapt_callback_example (t8_forest_t forest, t8_forest_t forest_from, t8_locidx_t which_tree, + [[maybe_unused]] t8_eclass_t tree_class, [[maybe_unused]] t8_locidx_t lelement_id, + [[maybe_unused]] const t8_scheme *scheme, const int is_family, + [[maybe_unused]] const int num_elements, t8_element_t *elements[]) +{ + const struct dummy_user_data *adapt_data = (const struct dummy_user_data *) t8_forest_get_user_data (forest); + t8_3D_point centroid; + t8_forest_element_centroid (forest_from, which_tree, elements[0], centroid.data ()); + double dist = t8_dist (centroid, adapt_data->midpoint); + if (dist < adapt_data->refine_if_inside_radius) { + /* Refine this element. */ + return 1; + } + else if (is_family && dist > adapt_data->coarsen_if_outside_radius) { + /* Coarsen this family. */ + return -1; + } + /* Do not change this element. */ + return 0; +} + +/** Test the adapt routine of a mesh handle. + * We compare the result to a classically adapted forest with similar callback. + */ +TEST (t8_gtest_handle_adapt, compare_adapt_with_forest) +{ + // Define forest, a mesh handle and user data. + const int level = 3; + t8_cmesh_t cmesh = t8_cmesh_new_hypercube_hybrid (sc_MPI_COMM_WORLD, 0, 0); + const t8_scheme *init_scheme = t8_scheme_new_default (); + t8_forest_t forest = t8_forest_new_uniform (cmesh, init_scheme, level, 0, sc_MPI_COMM_WORLD); + using mesh_class = t8_mesh_handle::mesh, dummy_user_data>; + mesh_class mesh_handle = mesh_class (forest); + struct dummy_user_data user_data = { + t8_3D_point ({ 0.5, 0.5, 1 }), /* Midpoints of the sphere. */ + 0.2, /* Refine if inside this radius. */ + 0.4 /* Coarsen if outside this radius. */ + }; + mesh_handle.set_user_data (&user_data); + + // Ref the forest as we want to keep using it after the adapt call to compare results. + t8_forest_ref (forest); + + // Adapt mesh handle and the forest with similar callbacks. + t8_mesh_handle::adapt_mesh (mesh_handle, refine_element_test, + coarsen_element_family_test, false); + forest = t8_forest_new_adapt (forest, forest_adapt_callback_example, 0, 1, &user_data); + + // Compare results. + EXPECT_TRUE (t8_forest_is_equal (mesh_handle.get_forest (), forest)); + // Clean up. + t8_forest_unref (&forest); +} diff --git a/test/mesh_handle/t8_gtest_cache_competence.cxx b/test/mesh_handle/t8_gtest_cache_competence.cxx index 15ae2e2038..d037af9138 100644 --- a/test/mesh_handle/t8_gtest_cache_competence.cxx +++ b/test/mesh_handle/t8_gtest_cache_competence.cxx @@ -98,7 +98,9 @@ class t8_gtest_cache_competence: public testing::Test { void TearDown () override { - t8_forest_unref (&forest); + if (forest->rc.refcount > 0) { + t8_forest_unref (&forest); + } } t8_forest_t forest; diff --git a/test/mesh_handle/t8_gtest_compare_handle_to_forest.cxx b/test/mesh_handle/t8_gtest_compare_handle_to_forest.cxx index 9a962fcc23..9bffb63ee5 100644 --- a/test/mesh_handle/t8_gtest_compare_handle_to_forest.cxx +++ b/test/mesh_handle/t8_gtest_compare_handle_to_forest.cxx @@ -78,7 +78,4 @@ TEST (t8_gtest_compare_handle_to_forest, compare_handle_to_forest) mesh_iterator++; } } - - // Unref the forest. - t8_forest_unref (&forest); } diff --git a/test/mesh_handle/t8_gtest_custom_competence.cxx b/test/mesh_handle/t8_gtest_custom_competence.cxx index bb73005e2a..6206d1cdc7 100644 --- a/test/mesh_handle/t8_gtest_custom_competence.cxx +++ b/test/mesh_handle/t8_gtest_custom_competence.cxx @@ -98,6 +98,7 @@ TEST (t8_gtest_custom_competence, custom_competence) EXPECT_EQ (level, it->get_level_dummy ()); } + t8_forest_ref (forest); // Test with two custom competences and a predefined competence. using competences = t8_mesh_handle::competence_pack; using mesh_class = t8_mesh_handle::mesh; @@ -114,7 +115,4 @@ TEST (t8_gtest_custom_competence, custom_competence) } EXPECT_TRUE (it->centroid_cache_filled ()); } - - // Unref the forest. - t8_forest_unref (&forest); } diff --git a/test/mesh_handle/t8_gtest_ghost.cxx b/test/mesh_handle/t8_gtest_ghost.cxx index 561e3b77a1..681afea2f7 100644 --- a/test/mesh_handle/t8_gtest_ghost.cxx +++ b/test/mesh_handle/t8_gtest_ghost.cxx @@ -55,7 +55,9 @@ class t8_mesh_ghost_test: public testing::TestWithParamrc.refcount > 0) { + t8_forest_unref (&forest); + } } t8_forest_t forest; int level; diff --git a/test/mesh_handle/t8_gtest_handle_data.cxx b/test/mesh_handle/t8_gtest_handle_data.cxx index 17929ccce9..82ef2e99de 100644 --- a/test/mesh_handle/t8_gtest_handle_data.cxx +++ b/test/mesh_handle/t8_gtest_handle_data.cxx @@ -38,12 +38,12 @@ along with t8code; if not, write to the Free Software Foundation, Inc., #include // --- Test for user data. --- -// Dummy user data taken from a tutorial for test purposes. +/** Dummy user data taken from a tutorial for test purposes. */ struct dummy_user_data { - t8_3D_point midpoint; /* The midpoint of our sphere. */ - double refine_if_inside_radius; /* if an element's center is smaller than this value, we refine the element. */ - double coarsen_if_outside_radius; /* if an element's center is larger this value, we coarsen its family. */ + t8_3D_point midpoint; /**< The midpoint of our sphere. */ + double refine_if_inside_radius; /**< if an element's center is smaller than this value, we refine the element. */ + double coarsen_if_outside_radius; /**< if an element's center is larger this value, we coarsen its family. */ }; /** Check that user data can be set and accesses for the handle. @@ -71,9 +71,6 @@ TEST (t8_gtest_handle_data, set_and_get_user_data) EXPECT_EQ (mesh_user_data.midpoint, user_data.midpoint); EXPECT_EQ (mesh_user_data.refine_if_inside_radius, user_data.refine_if_inside_radius); EXPECT_EQ (mesh_user_data.coarsen_if_outside_radius, user_data.coarsen_if_outside_radius); - - // Unref the forest. - t8_forest_unref (&forest); } // --- Test for element data. --- @@ -146,7 +143,4 @@ TEST (t8_gtest_handle_data, set_and_get_element_data) EXPECT_EQ (mesh[ighost].get_element_data ().volume, mesh[ighost].get_volume ()); } } - - // Unref the forest. - t8_forest_unref (&forest); } diff --git a/test/mesh_handle/t8_gtest_mesh_handle.cxx b/test/mesh_handle/t8_gtest_mesh_handle.cxx index 542dcc684a..12370ddb6e 100644 --- a/test/mesh_handle/t8_gtest_mesh_handle.cxx +++ b/test/mesh_handle/t8_gtest_mesh_handle.cxx @@ -52,7 +52,9 @@ class t8_mesh_handle_test: public testing::TestWithParamrc.refcount > 0) { + t8_forest_unref (&forest); + } } t8_forest_t forest; @@ -135,6 +137,7 @@ TEST_P (t8_mesh_handle_test, test_competences) } // --- Version with cached centroid variable. --- + t8_forest_ref (forest); using competence_centroid = t8_mesh_handle::competence_pack; using mesh_class_centroid = t8_mesh_handle::mesh; using element_class_centroid = typename mesh_class_centroid::element_class;