diff --git a/mesh_handle/competence_pack.hxx b/mesh_handle/competence_pack.hxx index d983d6819c..75950daec4 100644 --- a/mesh_handle/competence_pack.hxx +++ b/mesh_handle/competence_pack.hxx @@ -65,9 +65,9 @@ using all_cache_element_competences using cache_face_element_competences = element_competence_pack; -/** Predefined element competence pack combining all competences related to data. - * Please note that you must combine this with \ref t8_mesh_handle::data_mesh_competences. */ -using data_element_competences = element_competence_pack; +/** Predefined element data competence pack. + * Please note that you must combine this with \ref t8_mesh_handle::data_mesh_competences_basic. */ +using data_element_competences_basic = element_competence_pack; // --- Mesh competence pack. --- /** Class to pack different mesh competences into one template parameter for the \ref mesh class. @@ -91,11 +91,19 @@ struct mesh_competence_pack /** Empty competence pack. */ using empty_mesh_competences = mesh_competence_pack<>; -/** Predefined mesh competence pack combining all competences related to data. - * If you want to access the data also via the elements, combine this with \ref t8_mesh_handle::data_element_competences. +/** Predefined mesh competence pack to handle element data. + * If you want to access the data also via the elements, combine this with \ref t8_mesh_handle::data_element_competences_basic. */ template -using data_mesh_competences = mesh_competence_pack::template type>; +using data_mesh_competences_basic = mesh_competence_pack::template type>; + +/** Predefined mesh competence pack with the functionality to work with element data and to interpolate element data after adaptation. + * If you want to access the data also via the elements, combine this with \ref t8_mesh_handle::data_element_competences_basic. + */ +template +using interpolate_data_mesh_competence_pack + = mesh_competence_pack::template type, + interpolate_element_data_mesh_competence>; /** Predefined mesh competence pack combining all competences that are useful for discontinuous Galerkin methods. */ using dg_mesh_competences = mesh_competence_pack; diff --git a/mesh_handle/competences/element_data_competences.hxx b/mesh_handle/competences/element_data_competences.hxx index 3d52fda837..5ed7bf8ae7 100644 --- a/mesh_handle/competences/element_data_competences.hxx +++ b/mesh_handle/competences/element_data_competences.hxx @@ -22,26 +22,62 @@ /** \file element_data_competences.hxx * Handler for the element data of a \ref t8_mesh_handle::mesh. - * The file defines a mesh and an element competence for element data handling. - * Use both competences together if you want to manage element data for the elements of the mesh and access it directly for each element. + * The file defines mesh and element competences for element data handling. + * The mesh competences make it possible to manage element data and exchange it for ghost elements between processes. + * The element competences makes it possible to access these element data directly for each element of the mesh. + * A competence to interpolate data after adaptation using a user defined callback is provided. */ #pragma once #include #include #include +#include #include +#include #include #include +#include +#include +#include namespace t8_mesh_handle { + +/** Namespace detail to hide implementation details from the user. */ +namespace detail +{ +/** Helper function to wrap a span based interpolation callback (see \ref mesh::interpolate_callback_type) into + * the element-index based \ref interpolate_element_data_mesh_competence::internal_interpolate_callback_type. + * The returned wrapper receives the index/count pairs, builds the element spans, and forwards them to \a callback. + * The spans are built here and not in \ref interpolate_element_data_mesh_competence because \a TMesh is complete, + * so element_class is nameable — which it is not inside the competence (see note on + * \ref interpolate_element_data_mesh_competence::internal_interpolate_callback_type). + * This is used in \ref interpolate_element_data_mesh_competence::set_interpolate_callback + * \tparam TMesh The (complete) mesh handle type. + * \param [in] callback The span based user callback of type \ref mesh::interpolate_callback_type. Taken by value and + * moved into the returned wrapper, which owns it. + * \return Callback of type \ref interpolate_element_data_mesh_competence::internal_interpolate_callback_type. + */ +template +auto +to_replace_callback (typename TMesh::interpolate_callback_type callback) +{ + return + [callback = std::move (callback)] (const TMesh& mesh_old, TMesh& mesh_new, const int refine, const int num_old, + const t8_locidx_t first_old, const int num_new, const t8_locidx_t first_new) { + callback (mesh_old, mesh_new, refine, std::span (&mesh_old[first_old], num_old), + std::span (&mesh_new[first_new], num_new)); + }; +} +} // namespace detail + // --- Mesh competence for element data management. --- /** Handler for the element data of a \ref mesh. * Use this competence if you want to manage element data for the elements of the mesh. - * Use the helper \ref element_data_mesh_competence to get this competence with the correct template parameters form for the mesh. + * Use the helper \ref element_data_mesh_competence to get this competence with the correct template parameters form. * If you want to access the data not only in vector form but also directly for each element, - * you can combine this competence with the \ref element_data_element_competence competence. + * you can combine this competence with \ref element_data_element_competence. * In summary you can use the competences like this: * mesh, * mesh_competence_pack::template type>>; @@ -79,18 +115,30 @@ class element_data_mesh_competence_impl: public t8_crtp_basic { return m_element_data; } + /** Get the element data vector by moving it out of the competence. + * In contrast to \ref get_element_data, this transfers ownership of the data instead of + * returning a reference. After this call the internal element data vector is left empty (moved-from), + * so \ref set_element_data should be used before accessing the data again. + * \return Element data vector with data of Type TElementDataType, moved out of the competence. + */ + std::vector + take_element_data () + { + return std::move (m_element_data); + } + /** Exchange the element data for ghost elements between processes. * This routine has to be called on each process after setting the element data for all local elements. */ void exchange_ghost_data () { - // t8_forest_ghost_exchange_data expects an sc_array, so we need to wrap our data array to one. - sc_array* sc_array_wrapper; + // Extend element data array to hold also the ghost elements. const auto num_local_elements = this->underlying ().get_num_local_elements (); const auto num_ghosts = this->underlying ().get_num_ghosts (); m_element_data.resize (num_local_elements + num_ghosts); - sc_array_wrapper + // t8_forest_ghost_exchange_data expects an sc_array, so we need to wrap our data array to one. + sc_array* sc_array_wrapper = sc_array_new_data (m_element_data.data (), sizeof (ElementDataType), num_local_elements + num_ghosts); // Data exchange: entries with indices > num_local_elements will get overwritten. @@ -127,10 +175,9 @@ struct element_data_mesh_competence * \tparam TUnderlying Use the \ref element with specified competences as template parameter. */ template -struct element_data_element_competence: public t8_crtp_basic +struct element_data_element_competence: public t8_crtp_operator { public: - // --- Getter and setter for element data. --- /** Set the element data for the element. * \note You can only set element data for non-ghost elements. * \param [in] element_data The element data to be set of Type TMeshClass::ElementDataType. @@ -167,4 +214,124 @@ struct element_data_element_competence: public t8_crtp_basic } }; +// --- Mesh competence to interpolate data. --- +/** Mesh competence to interpolate the element data after an adaptation step. + * The \ref element_data_mesh_competence stores a vector of element data, but that data has to be updated if + * the mesh is adapted, since the elements it refers to are refined, coarsened or reordered. This competence adds the + * ability to interpolate the data after the adaptation via a user defined callback set using \ref set_interpolate_callback. + * The next \ref mesh::commit applies it. + * \note It therefore only makes sense in combination with the element data competence + * (see \ref interpolate_data_mesh_competence_pack, which bundles the two). + * \tparam TUnderlying Use the \ref mesh class here. + */ +template +class interpolate_element_data_mesh_competence: + public t8_crtp_operator { + public: + /** Mesh internal, element-index based storage for the interpolation callback. + * Users should use the easier span based callback type \ref mesh::interpolate_callback_type. + * \see set_interpolate_callback for registering a callback. + * The span based callback is automatically wrapped in \ref set_interpolate_callback to match this type and stored as + * \ref m_interpolate_callback to be used in the next \ref mesh::commit. + * \note We can not store or the span based \ref mesh::interpolate_callback_type directly. This competence uses the + * CRTP pattern, so while the competence is instantiated, the mesh (\a TUnderlying) is still an incomplete type. + * A data member of type \ref mesh::interpolate_callback_type, would require \c TUnderlying::element_class, + * which is not available for the incomplete type. Therefore we use this index based callback type for storage + * without the need for element_class. Using \ref set_interpolate_callback, we move the element_class lookup + * to the call site. + * \param [in] mesh_old The old mesh that is adapted from. + * \param [in,out] mesh_new The new mesh constructed from \a mesh_old. + * \param [in] refine -1 if a family in the old mesh got coarsened, 0 if the element was not touched, + * 1 if the element got refined. + * \param [in] num_old The number of outgoing elements. + * \param [in] first_old The local mesh handle index of the first outgoing element in the old mesh. + * \param [in] num_new The number of incoming elements. + * \param [in] first_new The local mesh handle index of the first incoming element in the new mesh. + + */ + using internal_interpolate_callback_type + = std::function; + + /** Register a user callback to interpolate the element data after adaptation. + * Note that data can only be interpolated for a level difference of at most one. + * Please use the type \ref mesh::interpolate_callback_type for the callback. + * \see mesh::interpolate_callback_type for the expected callback shape and the meaning of its arguments. + * \note This function is templated on purpose: it is only instantiated at the call site, where the mesh is a + * complete type and \ref mesh::element_class is nameable. The element class is needed in the definition + * of the interpolate_callback_type. You do not have to provide this template, it is normally auto deduced. + * \note This function is templated on purpose: it is only instantiated at the call site, where the mesh is a + * complete type and \ref mesh::element_class (needed to name \ref mesh::interpolate_callback_type) is + * nameable, which it is not inside this competence. + * The template parameter is deduced from the passed callback, so you do not have to provide it explicitly! + * \tparam TInterpolateCallback The user callback type \ref mesh::interpolate_callback_type. + * \param [in] interpolate_callback The span based interpolation callback. + * + */ + template + void + set_interpolate_callback (TInterpolateCallback&& interpolate_callback) + { + /* We wrap the user defined, span based callback using \ref detail::to_replace_callback to the index based callback + * type \ref internal_interpolate_callback_type to be able to store the callback without the need of knowing + * \ref mesh::element_class.*/ + m_interpolate_callback + = detail::to_replace_callback (std::forward (interpolate_callback)); + } + + protected: + /** Decide whether \ref mesh::set_partition has been requested for the upcoming \ref mesh::commit. + * With the interpolation competence the partition step is postponed so that it runs after the element data has been + * interpolated onto the new mesh; \ref mesh::set_partition therefore records its choice in + * \ref m_partition_for_coarsening instead of building the partitioned forest directly. + * \return true if \ref mesh::set_partition has been called (i.e. \ref m_partition_for_coarsening holds a value), + * false otherwise. + */ + bool + set_partition_called () + { + return m_partition_for_coarsening.has_value (); + } + + /** Repartition the element data so it follows a newly partitioned forest. + * The element data currently belongs to \a forest_from; this moves it to the layout of \a forest_to, which must + * have been created by partitioning \a forest_from. + * Analogous to \ref element_data_mesh_competence_impl::exchange_ghost_data, but using \ref t8_forest_partition_data. + * This function is called from \ref mesh::commit after the interpolated data has been produced and the partitioned + * forest has been committed. + * \param [in] forest_from The (committed) forest the current element data belongs to. + * \param [in] forest_to The committed forest that was partitioned from \a forest_from. + * \note Both forests could also be accessed directly (by this->underlying()) but this requires that the function is + * called on the exact right states of m_forest and m_uncommitted_forest. + * Providing the variables is the saver implementation. + */ + void + repartition_element_data (t8_forest_t forest_from, t8_forest_t forest_to) + { + using element_data_type = typename TUnderlying::ElementDataType; + // Take ownership of old data and wrap into sc_array. This is because the forest functions expect an sc_array. + std::vector old_data = this->underlying ().take_element_data (); + sc_array* data_in = sc_array_new_data (old_data.data (), sizeof (element_data_type), + t8_forest_get_local_num_leaf_elements (forest_from)); + const t8_locidx_t num_new_local = t8_forest_get_local_num_leaf_elements (forest_to); + // Define vector for the new data and wrap it. + std::vector partitioned_data (num_new_local); + sc_array* data_out = sc_array_new_data (partitioned_data.data (), sizeof (element_data_type), num_new_local); + // Partition magic using the forest function. + t8_forest_partition_data (forest_from, forest_to, data_in, data_out); + // Clean up. + sc_array_destroy (data_in); + sc_array_destroy (data_out); + // Set the partitioned element data to the mesh. + this->underlying ().set_element_data (std::move (partitioned_data)); + } + + internal_interpolate_callback_type + m_interpolate_callback; /**< The wrapped element-index based interpolation callback, + * applied on the next \ref mesh::commit. */ + std::optional + m_partition_for_coarsening; /**< Postponed \ref mesh::set_partition request: a value means partition on the next + * commit (with value passed to \ref mesh::set_partition); no value means do not partition. */ +}; + } // namespace t8_mesh_handle diff --git a/mesh_handle/element.hxx b/mesh_handle/element.hxx index 3b1854b452..5528086c42 100644 --- a/mesh_handle/element.hxx +++ b/mesh_handle/element.hxx @@ -55,6 +55,7 @@ namespace t8_mesh_handle * 2.) for the cached options to keep the number of member variables of the default element to a minimum to save memory. * The choice between calculate and cache is a tradeoff between runtime and memory usage. * + * \tparam TMeshClass The class of the mesh the element belongs to. * \tparam TCompetences The competences you want to add to the default functionality of the element. */ @@ -65,8 +66,7 @@ class element: public TCompetences>... { parameters specified. */ friend TMeshClass; /**< Define TMeshClass as friend to be able to access e.g. the constructor. */ friend struct element_data_element_competence< - SelfType>; /**< Define the competence to access element data as friend to - be able to access e.g. the mesh. */ + SelfType>; /**< Define the competence as friend to be able to access e.g. the mesh from competence. */ /** Private constructor for an element of a mesh. This could be a simple mesh element or a ghost element. * This constructor should only be called by the TMeshClass (and invisible for the user). diff --git a/mesh_handle/internal/adapt.hxx b/mesh_handle/internal/adapt.hxx index 1ca88e1fbb..4cb863ccaa 100644 --- a/mesh_handle/internal/adapt.hxx +++ b/mesh_handle/internal/adapt.hxx @@ -30,6 +30,7 @@ #include #include #include +#include #include #include @@ -68,7 +69,7 @@ struct mesh_adapt_context_base * Struct inherits from \ref mesh_adapt_context_base and implements the virtual adapt callback using the mesh and the callback. * \tparam TMeshClass The mesh handle class. */ -template +template struct mesh_adapt_context final: mesh_adapt_context_base { /** Constructor of the context with the mesh handle and the user defined callback. @@ -98,8 +99,8 @@ struct mesh_adapt_context final: mesh_adapt_context_base } private: - TMeshClass& m_mesh_handle; /**< The mesh handle to adapt. */ - typename TMeshClass::adapt_callback_type m_adapt_callback; /**< The adapt callback. */ + TMeshClass& m_mesh_handle; /**< The mesh handle to adapt. */ + const typename TMeshClass::adapt_callback_type m_adapt_callback; /**< The adapt callback. */ }; /** Registry pattern is used to register contexts, which provides access to the adapt callback and the mesh handle. diff --git a/mesh_handle/internal/interpolate.hxx b/mesh_handle/internal/interpolate.hxx new file mode 100644 index 0000000000..a206af8698 --- /dev/null +++ b/mesh_handle/internal/interpolate.hxx @@ -0,0 +1,223 @@ +/* + 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) 2026 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 interpolate.hxx + * This file provides helper functionality to interpolate the element data of a \ref t8_mesh_handle::mesh + * onto a newly adapted mesh according to a user defined callback. + * During adaptation the element data attached to the old elements has to be transferred to the new (refined, + * coarsened or unchanged) elements. t8code reports this old-to-new correspondence through + * \ref t8_forest_iterate_replace, which expects a plain C callback of type \ref t8_forest_replace_t. + * The helpers in this file bridge the gap between the mesh handle patterns and these functions, + * following the same registry pattern as \ref adapt.hxx. + */ + +#pragma once + +#include +#include +#include +#include +#include + +namespace t8_mesh_handle +{ + +/** Namespace detail to hide implementation details from the user. */ +namespace detail +{ + +/** Virtual base class for mesh interpolation contexts. + * We need this base class and not only \ref mesh_interpolate_context for the \ref interpolate_registry. + * interpolate_registry should not be templated because we need to access registered contexts in + * \ref mesh_replace_callback_wrapper, where we do not know the type of the mesh. Therefore, we work with a map of + * forests to instances of this (type erased) base class to remain template free. + */ +struct mesh_interpolate_context_base +{ + /** Virtual destructor for safe polymorphic deletion. + */ + virtual ~mesh_interpolate_context_base () = default; + + /** Pure virtual callback to interpolate the element data of one set of old elements onto a set of new elements. + * The indices refer to the flat, process local element numbering of the mesh handle (mesh handle id), + * not to the tree local numbering used by the forest. The conversion is done in \ref mesh_replace_callback_wrapper. + * \param [in] refine -1 if a family in the old mesh got coarsened, 0 if the element was not touched, + * 1 if the element got refined. + * \param [in] num_old The number of outgoing (old) elements. + * \param [in] first_old The local mesh handle id of the first outgoing element in the old mesh. + * \param [in] num_new The number of incoming (new) elements. + * \param [in] first_new The local mesh handle id of the first incoming element in the new mesh. + */ + virtual void + interpolate (const int refine, const int num_old, const t8_locidx_t first_old, const int num_new, + const t8_locidx_t first_new) + = 0; +}; + +/** Templated mesh interpolation context holding the old and new mesh handle and the user defined callback. + * Struct inherits from \ref mesh_interpolate_context_base and implements the virtual interpolate callback using the + * two mesh handles and the callback. + * Type erasure via the base class lets \ref interpolate_registry store this context without being templated on the + * mesh type. + * \tparam TMesh The mesh handle class. + */ +template +struct mesh_interpolate_context final: mesh_interpolate_context_base +{ + using callback_type = + typename TMesh::internal_interpolate_callback_type; /**< The user defined interpolate callback type. */ + + /** Constructor of the context with the old and new mesh handle and the user defined callback. + * \param [in] mesh_old The old mesh that is being adapted. Only read from during interpolation. + * \param [in, out] mesh_new The new mesh constructed from \a mesh_old. Written to during interpolation. + * \param [in] interpolate_callback The interpolate callback. Moved into the context. + */ + mesh_interpolate_context (const TMesh& mesh_old, TMesh& mesh_new, callback_type&& interpolate_callback) + : m_mesh_old (mesh_old), m_mesh_new (mesh_new), m_callback (std::forward (interpolate_callback)) + { + } + + /** Interpolation of one group of elements using the old and the new mesh and the user defined callback. + * This function is called by \ref mesh_replace_callback_wrapper for each group. + * \param [in] refine -1 if a family got coarsened, 0 if the element was not touched, 1 if it got refined. + * \param [in] num_old The number of outgoing (old) elements. + * \param [in] first_old The local mesh handle id of the first outgoing element in the old mesh. + * \param [in] num_new The number of incoming (new) elements. + * \param [in] first_new The local mesh handle id of the first incoming element in the new mesh. + */ + void + interpolate (const int refine, const int num_old, const t8_locidx_t first_old, const int num_new, + const t8_locidx_t first_new) override + { + // Check if the interpolate callback is set and call it using the old and new mesh handle. + T8_ASSERTF (m_callback, "No interpolate callback set."); + m_callback (m_mesh_old, m_mesh_new, refine, num_old, first_old, num_new, first_new); + } + + private: + const TMesh& m_mesh_old; /**< The old mesh to read the element data from. */ + TMesh& m_mesh_new; /**< The new mesh to write the interpolated element data to. */ + const callback_type m_callback; /**< The user defined interpolate callback. */ +}; + +/** Registry pattern is used to register contexts, which provide access to the interpolate callback and the mesh + * handles. This globally accessible static class is required to get the meshes and the callback in the forest + * replace callback \ref mesh_replace_callback_wrapper, as the predefined \ref t8_forest_replace_t header does not + * permit to pass these as function arguments. It uses the same idea as \ref adapt_registry. + */ +class interpolate_registry { + public: + /** Static function to register \a context using \a forest as identifier. + * This makes the context publicly available through the registry. + * \param [in] forest The forest identifier. In our case, this is the old forest we interpolate data from. + * ( Consistent to the adapt_registry.) + * \param [in] context The context to register. Use unique pointer to ensure proper memory management and ownership. + * \note We need the forest and not the mesh as key because this class must not be templated on the mesh type. + */ + static void + register_context (t8_forest_t forest, std::unique_ptr context) + { + auto& map = get_map (); + auto [it, inserted] = map.emplace (forest, std::move (context)); + if (!inserted) { + t8_global_errorf ("ERROR: Context already registered!"); + } + } + + /** Static function to unregister a context using \a forest as identifier. + * \param [in] forest The forest identifier. In our case, this is the old forest we interpolate from. + */ + static void + unregister_context (t8_forest_t forest) + { + auto& map = get_map (); + [[maybe_unused]] 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. In our case, this is the old forest we interpolate from. + * \return Pointer to the context registered with the id \a forest if found, nullptr otherwise. + */ + static mesh_interpolate_context_base* + 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 mesh_interpolate_context_base 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 mesh_interpolate_context_base references. + */ + static std::unordered_map>& + get_map () + { + static std::unordered_map> map; + return map; + } +}; + +/** Wrapper around the mesh handle interpolation functionality to be able to pass the callback to the classic replace + * routine \ref t8_forest_iterate_replace of a forest. The function header fits the definition of + * \ref t8_forest_replace_t. + * \param [in] forest_old The forest that is adapted fitting the old element data. + * \param [in, out] forest_new The forest that is newly constructed from \a forest_old. + * Data will be interpolated to match this forest. + * \param [in] which_tree The local tree containing \a first_outgoing and \a first_incoming. + * \param [in] tree_class Unused; The eclass of the local tree containing \a first_outgoing and \a first_incoming. + * \param [in] scheme Unused; The scheme of the forest. + * \param [in] refine -1 if a family got coarsened, 0 if the element was not touched, 1 if it got refined. + * \param [in] num_outgoing The number of outgoing elements. + * \param [in] first_outgoing The tree local index of the first outgoing element. + * 0 <= first_outgoing < which_tree->num_elements + * \param [in] num_incoming The number of incoming elements. + * \param [in] first_incoming The tree local index of the first incoming element. + * 0 <= first_incom < new_which_tree->num_elements + */ +inline void +mesh_replace_callback_wrapper (t8_forest_t forest_old, t8_forest_t forest_new, t8_locidx_t which_tree, + [[maybe_unused]] const t8_eclass_t tree_class, + [[maybe_unused]] const t8_scheme_c* scheme, const int refine, const int num_outgoing, + const t8_locidx_t first_outgoing, const int num_incoming, + const t8_locidx_t first_incoming) +{ + // Get the static interpolate context from the registry. + // Via this, we can access the old and new mesh handle and the user defined interpolate callback. + auto* context = interpolate_registry::get (forest_old); + if (!context) { + t8_global_productionf ("Interpolate context not found. Did you forget to register it?"); + return; + } + + // Convert the tree local indices reported by the forest to the flat, process local indices used in the mesh handle + // (the mesh handle id). + const t8_locidx_t first_old_global = t8_forest_get_tree_element_offset (forest_old, which_tree) + first_outgoing; + const t8_locidx_t first_new_global = t8_forest_get_tree_element_offset (forest_new, which_tree) + first_incoming; + // Call the actual interpolate callback stored in the context. + context->interpolate (refine, num_outgoing, first_old_global, num_incoming, first_new_global); +} + +} // namespace detail +} // namespace t8_mesh_handle diff --git a/mesh_handle/mesh.hxx b/mesh_handle/mesh.hxx index b08fba5e78..e9fbfd6c04 100644 --- a/mesh_handle/mesh.hxx +++ b/mesh_handle/mesh.hxx @@ -30,12 +30,14 @@ #include "element.hxx" #include "competence_pack.hxx" #include "internal/adapt.hxx" +#include "internal/interpolate.hxx" #include "competences/element_data_competences.hxx" #include "concepts.hxx" #include #include #include #include +#include #include #include #include @@ -77,6 +79,7 @@ class mesh: public TMeshCompetencePack::template apply::iterator; /**< Non-const iterator type for the mesh elements. */ friend struct element_data_element_competence; /**< Friend struct to access its element data vector. */ + // --- Definition of callback types. --- /** Callback function prototype to decide for refining and coarsening of a family of elements * or one element in a mesh handle. * If \a elements contains more than one element, they must form a family and we decide whether this family should be @@ -110,6 +113,49 @@ class mesh: public TMeshCompetencePack::template apply elements, TUserDataType user_data)>; + /** Callback function prototype to interpolate the element data after refining or coarsening. + * \note You need to include \ref interpolate_element_data_mesh_competence to you competences to be able to + * interpolate. The best way to do this is via the predefined pack \ref interpolate_data_mesh_competence_pack + * defined in \ref competence_pack.hxx. + * + * For each group of elements that changed during adaption, the outgoing elements of the old mesh are passed in + * \a old_elements and the incoming elements of the new mesh in \a new_elements; the callback reads the old data + * and writes the interpolated data onto the new elements. \a refine is the value \ref adapt_callback_type returned + * for this group. + * \see interpolate_element_data_mesh_competence::set_interpolate_callback for the usage of this callback. + * \param [in] mesh_old The old mesh that is adapted from. + * \param [in,out] mesh_new The new mesh constructed from \a mesh_old. + * \param [in] refine -1 if the family \a old_elements got coarsened, 0 if the element was not touched, + * 1 if the element got refined. Same convention as the return of \ref adapt_callback_type. + * \param [in] old_elements Span over the outgoing elements: the whole family on coarsening, + * a single element if refined or untouched. + * \param [in,out] new_elements Span over the incoming elements to write the interpolated data to: the children on + * refinement, a single element if coarsened or untouched. + */ + using interpolate_callback_type + = std::function old_elements, std::span new_elements)>; + + /** Templated callback function prototype to interpolate the element data after refining or coarsening, + * including user data. + * See the version without user_data \ref interpolate_callback_type for more details! + * Use \ref mesh_interpolate_callback_wrapper to convert this type into \ref interpolate_callback_type + * to be able to pass the callback to \ref interpolate_element_data_mesh_competence::set_interpolate_callback + * (see \ref element_data_competences.hxx). + * \tparam TUserDataType The type of the user data to be passed to the callback. + * \param [in] mesh_old The old mesh that is adapted from. + * \param [in,out] mesh_new The new mesh constructed from \a mesh_old. + * \param [in] refine -1 if the family got coarsened, 0 if the element was not touched, 1 if it got refined. + * \param [in] old_elements Span over the outgoing elements from \a mesh_old. + * \param [in,out] new_elements Span over the incoming elements to write the interpolated data to from \a mesh_new. + * \param [in] user_data The user data to be used during the interpolation. + */ + template + using interpolate_callback_type_with_userdata = std::function old_elements, + std::span new_elements, TUserDataType user_data)>; + + // --- Constructor and destructor. --- /** * Constructor for a mesh of the handle. * \param [in] forest The forest from which the mesh should be created. @@ -141,6 +187,16 @@ class mesh: public TMeshCompetencePack::template apply + static interpolate_callback_type + mesh_interpolate_callback_wrapper ( + interpolate_callback_type_with_userdata interpolate_callback_with_userdata, + const TUserDataType& user_data) + { + return [=] (const SelfType& mesh_old, SelfType& mesh_new, const int refine, + std::span old_elements, std::span new_elements) { + return interpolate_callback_with_userdata (mesh_old, mesh_new, refine, old_elements, new_elements, user_data); + }; + } + /** Wrapper to convert an adapt callback with user data of type \ref adapt_callback_type_with_userdata * into a callback without user data of type \ref adapt_callback_type using the defined user data \a user_data. * This is required to pass an adapt callback with user data to \ref set_adapt. @@ -298,7 +376,7 @@ class mesh: public TMeshCompetencePack::template applym_partition_for_coarsening = set_for_coarsening; + return; + } if (!m_uncommitted_forest.has_value ()) { t8_forest_t new_forest; t8_forest_init (&new_forest); @@ -338,7 +422,7 @@ class mesh: public TMeshCompetencePack::template applym_interpolate_callback) { + // Create new intermediate mesh to interpolate the data from the current mesh to the new mesh. + SelfType new_mesh (m_uncommitted_forest.value ()); + t8_forest_ref (m_uncommitted_forest.value ()); + // Register the interpolate context with the callback for the new mesh. With this, the standard + // iterate replace can be called. + detail::interpolate_registry::register_context ( + m_forest, std::make_unique> ( + *this, new_mesh, std::move (this->m_interpolate_callback))); + t8_forest_iterate_replace (m_uncommitted_forest.value (), m_forest, detail::mesh_replace_callback_wrapper); + detail::interpolate_registry::unregister_context (m_forest); + // Override the element data of the current mesh with the interpolated data from the "new mesh". + this->m_element_data = new_mesh.take_element_data (); + // Now we update the forest of the current mesh with the new forest and partition it if required. + t8_forest_unref (&m_forest); + if (this->set_partition_called ()) { + t8_forest_init (&m_forest); + t8_forest_set_partition (m_forest, m_uncommitted_forest.value (), + this->m_partition_for_coarsening.value ()); + const t8_ghost_type_t ghost_type = m_uncommitted_forest.value ()->ghost_type; + if (ghost_type != T8_GHOST_NONE) { + t8_forest_set_ghost (m_forest, true, ghost_type); + } + t8_forest_commit (m_forest); + + /* Now we repartition also the data: The interpolated data follows m_uncommitted_forest. + * We align it now with the partitioned m_forest. */ + this->repartition_element_data (m_uncommitted_forest.value (), m_forest); + } + else { + // Update underlying forest of the mesh for the case where we do not repartition. + m_forest = m_uncommitted_forest.value (); + } + // Cleanup and update the elements of the mesh. + m_uncommitted_forest.reset (); + update_elements (); + return; + } + else { + t8_global_infof ("No interpolation context set.\n"); + } + } + else { + t8_global_infof ("The element data was not interpolated during adaptation. Use set_element_data() to provide " + "new data or use the mesh competence interpolate_element_data_mesh_competence.\n"); + } } } t8_forest_unref (&m_forest); @@ -421,6 +554,15 @@ class mesh: public TMeshCompetencePack::template apply #include +#include #include #include @@ -33,32 +34,11 @@ along with t8code; if not, write to the Free Software Foundation, Inc., #include /** Store the rank on each element. */ -struct data_per_element +struct rank_data_per_element { int rank; ///< Rank of the element. }; -/** Callback function for the mesh handle to decide for refining or coarsening of (a family of) elements. - * The adaptation criterion is to refine every element with even id. - * The function header fits the definition of \ref TMesh::adapt_callback_type. - * \tparam TMeshClass The mesh handle class. - * \param [in] mesh The mesh that should be adapted. - * \param [in] elements One element or a family of elements to consider for adaptation. - * \return 1 if the first entry in \a elements should be refined, - * -1 if the family \a elements shall be coarsened, - * 0 else. - */ -template -int -mesh_adapt_callback_test_refine_second ([[maybe_unused]] const TMeshClass& mesh, - std::span elements) -{ - if ((elements[0].get_element_handle_id ()) % 2 == 0) { - return 1; - } - return 0; -} - /** Check the competence remote_ranks_mesh_competence for correctness. The ranks are set as data first, exchanged for * ghost elements and then checked against the competence functionality. */ @@ -66,11 +46,11 @@ TEST (t8_gtest_dg_competences, remote_ranks) { const int level = 2; using namespace t8_mesh_handle; - using mesh_class - = mesh, - data_mesh_competences>>; + using mesh_class = mesh, + data_mesh_competences_basic>>; auto mesh = handle_hypercube_hybrid_uniform_default (level, sc_MPI_COMM_WORLD, true, false); - mesh->set_adapt (mesh_adapt_callback_test_refine_second); + mesh->set_adapt (adapt_callback_refine_second); mesh->set_partition (); mesh->set_ghost (); mesh->commit (); @@ -87,7 +67,7 @@ TEST (t8_gtest_dg_competences, remote_ranks) SC_CHECK_MPI (mpiret); // Set local rank for all local mesh elements. - std::vector element_data (num_local, { mpirank }); + std::vector element_data (num_local, { mpirank }); mesh->set_element_data (std::move (element_data)); // Get element data and check that the remote ranks competence works as expected. mesh->exchange_ghost_data (); @@ -111,7 +91,7 @@ TEST (t8_gtest_dg_competences, face_vector_mesh_competence) using namespace t8_mesh_handle; using mesh_class = mesh, dg_mesh_competences>; auto mesh = handle_hypercube_hybrid_uniform_default (level, sc_MPI_COMM_WORLD, true, false); - mesh->set_adapt (mesh_adapt_callback_test_refine_second); + mesh->set_adapt (adapt_callback_refine_second); mesh->set_partition (); mesh->set_ghost (); mesh->commit (); diff --git a/test/mesh_handle/competences/t8_gtest_handle_data.cxx b/test/mesh_handle/competences/t8_gtest_handle_data.cxx index cfb8b60b44..c9f66e2ea0 100644 --- a/test/mesh_handle/competences/t8_gtest_handle_data.cxx +++ b/test/mesh_handle/competences/t8_gtest_handle_data.cxx @@ -3,7 +3,7 @@ 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 +Copyright (C) 2026 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 @@ -27,6 +27,7 @@ along with t8code; if not, write to the Free Software Foundation, Inc., #include #include +#include #include #include @@ -36,19 +37,12 @@ along with t8code; if not, write to the Free Software Foundation, Inc., #include #include -/** Dummy element data taken from a tutorial for test purposes. */ -struct data_per_element -{ - int level; - double volume; -}; - /** Check that element data can be set for the handle and that exchanging data for the ghosts works. */ TEST (t8_gtest_handle_data, set_and_get_element_data) { const int level = 2; - using mesh_class = t8_mesh_handle::mesh>; + using mesh_class = t8_mesh_handle::mesh>; auto mesh = t8_mesh_handle::handle_hypercube_hybrid_uniform_default (level, sc_MPI_COMM_WORLD, true, false); @@ -56,14 +50,15 @@ TEST (t8_gtest_handle_data, set_and_get_element_data) // Ensure that we actually test with ghost elements. EXPECT_GT (mesh->get_num_ghosts (), 0); } + EXPECT_TRUE (mesh->has_element_data_handler_competence ()); - // Create element data for all local mesh elements. + // Create element data for all local mesh elements and set via mesh competence. std::vector element_data; for (const auto &elem : *mesh) { element_data.push_back ({ elem.get_level (), elem.get_volume () }); } mesh->set_element_data (std::move (element_data)); - // Get element data and check that the data for all elements (including ghosts) is correct. + // Exchange element data for ghosts and check that the data for all elements (including ghosts) is correct. mesh->exchange_ghost_data (); auto mesh_element_data = mesh->get_element_data (); for (t8_locidx_t ielem = 0; ielem < mesh->get_num_local_elements () + mesh->get_num_ghosts (); ielem++) { @@ -71,8 +66,7 @@ TEST (t8_gtest_handle_data, set_and_get_element_data) EXPECT_EQ (mesh_element_data[ielem].volume, (*mesh)[ielem].get_volume ()) << "ielem = " << ielem; } - // Modify element data for elements that are in the first half of the global trees. - EXPECT_TRUE (mesh->has_element_data_handler_competence ()); + // Modify element data via the element competence for elements that are in the first half of the global trees. auto forest = mesh->get_forest (); t8_gloidx_t barrier = t8_forest_get_num_global_trees (forest) / 2.0; const int newlevel = 42; @@ -83,6 +77,7 @@ TEST (t8_gtest_handle_data, set_and_get_element_data) elem.set_element_data (elem_data); } } + // Exchange data for ghosts and check that the data for all elements (including ghosts) is correct. mesh->exchange_ghost_data (); for (const auto &elem : *mesh) { if (t8_forest_global_tree_id (forest, elem.get_local_tree_id ()) < barrier) { @@ -117,10 +112,11 @@ TEST (t8_gtest_handle_data, set_and_get_element_data) TEST (t8_gtest_handle_data, test_union_mesh_competence_pack) { using namespace t8_mesh_handle; - using mesh_class = mesh< - union_competence_packs_type, - union_competence_packs_type, data_mesh_competences, - empty_mesh_competences>>; + using mesh_class + = mesh, + union_competence_packs_type, + data_mesh_competences_basic, empty_mesh_competences>>; EXPECT_TRUE (mesh_class::has_element_data_handler_competence ()); using element_class = typename mesh_class::element_class; diff --git a/test/mesh_handle/t8_gtest_adapt_partition_balance.cxx b/test/mesh_handle/t8_gtest_adapt_partition_balance.cxx index f173729259..3edf845561 100644 --- a/test/mesh_handle/t8_gtest_adapt_partition_balance.cxx +++ b/test/mesh_handle/t8_gtest_adapt_partition_balance.cxx @@ -24,117 +24,15 @@ along with t8code; if not, write to the Free Software Foundation, Inc., * \file t8_gtest_adapt_partition_balance.cxx * Tests for the adapt, partition and balance routines of mesh handle. */ -#include "t8_types/t8_vec.h" #include #include +#include "t8_gtest_common.hxx" #include -#include -#include #include #include #include #include -#include - -/** Dummy user data taken from tutorial for test purposes. */ -struct dummy_user_data -{ - t8_3D_vec 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 for the mesh handle to decide for refining or coarsening of (a family of) elements. - * 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. - * The function header fits the definition of \ref TMesh::adapt_callback_type_with_userdata. - * \tparam TMeshClass The mesh handle class. - * \param [in] mesh The mesh that should be adapted. - * \param [in] elements One element or a family of elements to consider for adaptation. - * \param [in] user_data The user data to be used during the adaptation process. - * \return 1 if the first entry in \a elements should be refined, - * -1 if the family \a elements shall be coarsened, - * 0 else. - */ -template -int -adapt_callback_test ([[maybe_unused]] const TMeshClass &mesh, - std::span elements, const dummy_user_data &user_data) -{ - const auto element_centroid = elements[0].get_centroid (); - const double dist = t8_dist (element_centroid, user_data.midpoint); - if (dist < user_data.refine_if_inside_radius) { - return 1; - } - // Check if we got a family and if yes, if we should coarsen. - if ((elements.size () > 1) && (dist > user_data.coarsen_if_outside_radius)) { - return -1; - } - return 0; -} - -/** Adapt callback implementation for a forest. - * This callback defines the same adaptation rules as \ref adapt_callback_test, - * but it is used for the forest instead of the 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_vec 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) { - return 1; - } - else if (is_family && dist > adapt_data->coarsen_if_outside_radius) { - return -1; - } - return 0; -} - -//--- Second callback type for testing purpose: Refine every second element. --- -/** Callback function for the mesh handle to decide for refining or coarsening of (a family of) elements. - * The adaptation criterion is to refine every element with even id. - * The function header fits the definition of \ref TMesh::adapt_callback_type. - * \tparam TMeshClass The mesh handle class. - * \param [in] mesh The mesh that should be adapted. - * \param [in] elements One element or a family of elements to consider for adaptation. - * \return 1 if the first entry in \a elements should be refined, - * -1 if the family \a elements shall be coarsened, - * 0 else. - */ -template -int -mesh_adapt_callback_test_refine_second ([[maybe_unused]] const TMeshClass &mesh, - std::span elements) -{ - if ((elements[0].get_element_handle_id ()) % 2 == 0) { - return 1; - } - return 0; -} - -/** Adapt callback implementation for a forest. The adaptation criterion is to refine every element with even id. - * This callback defines the same adaptation rules as \ref mesh_adapt_callback_test_refine_second, - * but it is used for the forest instead of the mesh handle. - */ -int -forest_adapt_callback_refine_second ([[maybe_unused]] t8_forest_t forest, [[maybe_unused]] 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, - [[maybe_unused]] const int is_family, [[maybe_unused]] const int num_elements, - [[maybe_unused]] t8_element_t *elements[]) -{ - if ((t8_forest_get_tree_element_offset (forest_from, which_tree) + lelement_id) % 2 == 0) { - return 1; - } - return 0; -} /** Test the adapt, partition and balance routines of a mesh handle. * The test compares the results of the mesh handle to a forest adapted with the same criterion and balanced and partitioned similarly. @@ -151,10 +49,10 @@ TEST (t8_gtest_handle_adapt, compare_with_forest) t8_forest_t forest = t8_forest_new_uniform (cmesh, init_scheme, level, 0, sc_MPI_COMM_WORLD); using mesh_class = t8_mesh_handle::mesh<>; mesh_class mesh_handle = mesh_class (forest); - struct dummy_user_data user_data = { - t8_3D_vec ({ 0.5, 0.5, 1 }), /**< Midpoints of the sphere. */ - 0.2, /**< Refine if inside this radius. */ - 0.4 /**< Coarsen if outside this radius. */ + dummy_user_data user_data { + t8_3D_vec { 0.5, 0.5, 1 }, /**< Midpoints of the sphere. */ + 0.2, /**< Refine if inside this radius. */ + 0.4 /**< Coarsen if outside this radius. */ }; // Ref the forest as we want to keep using it after the adapt call to compare results. @@ -162,10 +60,10 @@ TEST (t8_gtest_handle_adapt, compare_with_forest) // Adapt mesh handle. mesh_handle.set_adapt ( - mesh_class::mesh_adapt_callback_wrapper (adapt_callback_test, user_data)); + mesh_class::mesh_adapt_callback_wrapper (adapt_callback_sphere, user_data)); mesh_handle.commit (); // Adapt forest classically. - forest = t8_forest_new_adapt (forest, forest_adapt_callback_example, 0, 0, &user_data); + forest = t8_forest_new_adapt (forest, forest_adapt_callback_sphere, 0, 0, &user_data); // Compare results. EXPECT_TRUE (t8_forest_is_equal (mesh_handle.get_forest (), forest)); @@ -174,7 +72,7 @@ TEST (t8_gtest_handle_adapt, compare_with_forest) mesh_handle.set_balance (); mesh_handle.set_partition (); mesh_handle.set_adapt ( - mesh_class::mesh_adapt_callback_wrapper (adapt_callback_test, user_data)); + mesh_class::mesh_adapt_callback_wrapper (adapt_callback_sphere, user_data)); mesh_handle.commit (); EXPECT_TRUE (mesh_handle.is_balanced ()); @@ -182,14 +80,14 @@ TEST (t8_gtest_handle_adapt, compare_with_forest) t8_forest_t forest_compare; t8_forest_init (&forest_compare); t8_forest_set_user_data (forest_compare, &user_data); - t8_forest_set_adapt (forest_compare, forest, forest_adapt_callback_example, false); + t8_forest_set_adapt (forest_compare, forest, forest_adapt_callback_sphere, false); t8_forest_set_partition (forest_compare, NULL, false); t8_forest_set_balance (forest_compare, NULL, false); t8_forest_commit (forest_compare); EXPECT_TRUE (t8_forest_is_equal (mesh_handle.get_forest (), forest_compare)); // Adapt again with the second callback. - mesh_handle.set_adapt (mesh_adapt_callback_test_refine_second); + mesh_handle.set_adapt (adapt_callback_refine_second); mesh_handle.commit (); t8_forest_t forest_refine; diff --git a/test/mesh_handle/t8_gtest_common.hxx b/test/mesh_handle/t8_gtest_common.hxx new file mode 100644 index 0000000000..5405cf0fdd --- /dev/null +++ b/test/mesh_handle/t8_gtest_common.hxx @@ -0,0 +1,172 @@ +/* +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) 2026 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_common.hxx + * Collection of data types and callbacks for mesh handle meshes and forests to be used in the mesh handle test cases. + */ +#pragma once + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +/** Dummy element data taken from a tutorial for test purposes. */ +struct data_per_element +{ + int level; /**< The level of the element. */ + double volume; /**< The volume of the element. */ + + /** Comparison operator to check if two data entries are the same. */ + bool + operator== (const data_per_element &) const + = default; +}; + +/** Adapt callback for a mesh handle hypercube that exercises all possible operations (coarsen, refine, nothing). + * Adaption with this callback produces an unbalanced mesh (if the initial refinement level fits). + * The callback coarsens the left half (x < 0.5), refines a band 0.5 <= x < 0.75 and leaves the rest untouched. + * The centroid is used to determine the x-coordinate. + * \tparam TMeshClass The mesh handle class. + * \param [in] mesh Unused. + * \param [in] elements One element, or a family of elements to consider for adaptation. + * \return -1 to coarsen the family, 1 to refine the first element, 0 otherwise. + */ +template +int +adapt_callback_coarsen_left_refine_middle ([[maybe_unused]] const TMeshClass &mesh, + std::span elements) +{ + const double x_centroid = elements[0].get_centroid ()[0]; + if ((elements.size () > 1) && (x_centroid < 0.5)) { + return -1; // Coarsen the left half. + } + if (x_centroid >= 0.5 && x_centroid < 0.75) { + return 1; // Refine band 0.5 <= x < 0.75. + } + return 0; // Untouched right band. +} + +//--- Callback to refine every second element. --- +/** Callback function for the mesh handle to decide for refining or coarsening of (a family of) elements. + * The adaptation criterion is to refine every element with even id. + * The function header fits the definition of \ref t8_mesh_handle::mesh::adapt_callback_type. + * \tparam TMeshClass The mesh handle class. + * \param [in] mesh The mesh that should be adapted. + * \param [in] elements One element or a family of elements to consider for adaptation. + * \return 1 if the first entry in \a elements should be refined, + * -1 if the family \a elements shall be coarsened, + * 0 else. + */ +template +int +adapt_callback_refine_second ([[maybe_unused]] const TMeshClass &mesh, + std::span elements) +{ + if ((elements[0].get_element_handle_id ()) % 2 == 0) { + return 1; + } + return 0; +} + +/** Adapt callback implementation for a forest. The adaptation criterion is to refine every element with even id. + * This callback defines the same adaptation rules as \ref adapt_callback_refine_second, + * but it is used for the forest instead of the mesh handle. + */ +int +forest_adapt_callback_refine_second ([[maybe_unused]] t8_forest_t forest, [[maybe_unused]] 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, + [[maybe_unused]] const int is_family, [[maybe_unused]] const int num_elements, + [[maybe_unused]] t8_element_t *elements[]) +{ + if ((t8_forest_get_tree_element_offset (forest_from, which_tree) + lelement_id) % 2 == 0) { + return 1; + } + return 0; +} + +//--- Define callbacks to adapt according to a sphere-based criterion. --- +/** Dummy user data taken from tutorial for test purposes. This user data type is used in the sphere adapt callbacks. */ +struct dummy_user_data +{ + t8_3D_vec 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 for the mesh handle to decide for refining or coarsening of (a family of) elements. + * The function header fits the definition of \ref t8_mesh_handle::mesh::adapt_callback_type_with_userdata. + * \tparam TMeshClass The mesh handle class. + * \param [in] mesh The mesh that should be adapted. + * \param [in] elements One element or a family of elements to consider for adaptation. + * \param [in] user_data The user data to be used during the adaptation process. + * \return 1 if the first entry in \a elements should be refined, + * -1 if the family \a elements shall be coarsened, + * 0 else. + */ +template +int +adapt_callback_sphere ([[maybe_unused]] const TMeshClass &mesh, + std::span elements, const dummy_user_data &user_data) +{ + auto element_centroid = elements[0].get_centroid (); + double dist = t8_dist (element_centroid, user_data.midpoint); + if (dist < user_data.refine_if_inside_radius) { + return 1; + } + // Check if we got a family and if yes, if we should coarsen. + if ((elements.size () > 1) && (dist > user_data.coarsen_if_outside_radius)) { + return -1; + } + return 0; +} + +/** Adapt callback implementation for a forest. + * This callback defines the same adaptation rules as adapt_callback_sphere defined above, + * but it is used for the forest instead of the mesh handle. + */ +int +forest_adapt_callback_sphere (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_vec 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) { + return 1; + } + else if (is_family && dist > adapt_data->coarsen_if_outside_radius) { + return -1; + } + return 0; +} diff --git a/test/mesh_handle/t8_gtest_interpolate.cxx b/test/mesh_handle/t8_gtest_interpolate.cxx new file mode 100644 index 0000000000..1d2e46ac7c --- /dev/null +++ b/test/mesh_handle/t8_gtest_interpolate.cxx @@ -0,0 +1,150 @@ +/* +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) 2026 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_interpolate.cxx + * Tests to check that the data interpolation for the mesh handle works as intended. + */ + +#include +#include +#include "t8_gtest_common.hxx" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +/** Dummy user data for the interpolation. */ +struct interpolate_user_data +{ + int level_step; /**< Levels added when refining and subtracted when coarsening (1 in the standard case). + This is to be applied to the level entry of the dummy user data. */ +}; + +/** Interpolation callback for the mesh handle, using user data. + * The function header fits the definition of \ref TMesh::interpolate_callback_type_with_userdata. + * Copies the data of untouched elements, averages the parent volume over the children on refinement, and sums the + * children's volume onto the parent on coarsening. The level changes by \a user_data.level_step per step. + * \tparam TMeshClass The mesh handle class. + * \param [in] mesh_old The old mesh that is adapted from. + * \param [in,out] mesh_new The new mesh constructed from \a mesh_old. + * \param [in] refine -1 if the family got coarsened, 0 if the element was not touched, 1 if it got refined. + * \param [in] old_elements Span over the outgoing elements from \a mesh_old. + * \param [in,out] new_elements Span over the incoming elements to write the interpolated data to from \a mesh_new. + * \param [in] user_data The user data to be used during the interpolation. + */ +template +void +interpolate_callback ([[maybe_unused]] const TMeshClass& mesh_old, [[maybe_unused]] TMeshClass& mesh_new, + const int refine, std::span old_elements, + std::span new_elements, + const interpolate_user_data& user_data) +{ + /* Untouched: copy data. */ + if (refine == 0) { + new_elements[0].set_element_data (old_elements[0].get_element_data ()); + } + /* Refined: children share the parent volume equally, level increases by user_data.level_step. */ + else if (refine == 1) { + const auto& parent_data = old_elements[0].get_element_data (); + for (auto& child : new_elements) { + child.set_element_data ( + data_per_element { parent_data.level + user_data.level_step, parent_data.volume / new_elements.size () }); + } + } + /* Coarsened: parent volume is the sum of the children, level decreases by user_data.level_step. */ + else if (refine == -1) { + double tmp_volume = 0; + for (const auto& child : old_elements) { + tmp_volume += child.get_element_data ().volume; + } + new_elements[0].set_element_data ( + data_per_element { old_elements[0].get_element_data ().level - user_data.level_step, tmp_volume }); + } +} + +/** Test for the data interpolation of the mesh handle using the interpolate_callback. */ +TEST (t8_gtest_handle_interpolate, test_interpolate_data) +{ + const int level + = 3; ///< 3 is the minimum level to have each case of refine, coarsen, copy but still include balance. + using mesh_class = t8_mesh_handle::mesh>; + auto mesh = t8_mesh_handle::handle_hypercube_uniform_default (T8_ECLASS_HEX, level, sc_MPI_COMM_WORLD); + + interpolate_user_data dummy_interpolate_user_data { 1 }; ///< Level changes by one per adaptation step. + + // Create element data for all local mesh elements and set via mesh competence. + std::vector element_data; + for (const auto& elem : *mesh) { + element_data.push_back ({ elem.get_level (), elem.get_volume () }); + } + mesh->set_element_data (std::move (element_data)); + + // Adapt the mesh and set all options. + mesh->set_adapt (adapt_callback_coarsen_left_refine_middle); + mesh->set_balance (); + mesh->set_partition (); + mesh->set_ghost (); + mesh->set_interpolate_callback (mesh_class::mesh_interpolate_callback_wrapper ( + interpolate_callback, dummy_interpolate_user_data)); + mesh->commit (); + + // Check basics. + EXPECT_TRUE (mesh->is_balanced ()); + EXPECT_TRUE (mesh->get_num_ghosts () > 0); + // Ensure partitioned. + int mpi_size = 0; + int mpiret = sc_MPI_Comm_size (sc_MPI_COMM_WORLD, &mpi_size); + SC_CHECK_MPI (mpiret); + int num_global_elements_averaged = (int) (mesh->get_num_global_elements () / mpi_size); + EXPECT_LE (mesh->get_num_local_elements (), num_global_elements_averaged + 1); + EXPECT_GE (mesh->get_num_local_elements (), num_global_elements_averaged - 1); + + // Test interpolation. + // Variables to demonstrate that we tested interpolation for coarsening and for refinement. + int found_refined = false; + int found_coarsened = false; + for (auto& elem : *mesh) { + if (!found_refined && (elem.get_level () > level)) { + found_refined = true; + } + if (!found_coarsened && (elem.get_level () < level)) { + found_coarsened = true; + } + // Check that element data and actual geometric data match. + EXPECT_EQ (elem.get_level (), elem.get_element_data ().level); + // For hexes, our interpolation method for volume is accurate. + EXPECT_NEAR (elem.get_volume (), elem.get_element_data ().volume, T8_PRECISION_SQRT_EPS); + } + int global_refined = 0, global_coarsened = 0; + sc_MPI_Allreduce (&found_refined, &global_refined, 1, sc_MPI_INT, sc_MPI_LOR, sc_MPI_COMM_WORLD); + sc_MPI_Allreduce (&found_coarsened, &global_coarsened, 1, sc_MPI_INT, sc_MPI_LOR, sc_MPI_COMM_WORLD); + EXPECT_TRUE (global_refined); + EXPECT_TRUE (global_coarsened); +} diff --git a/test/mesh_handle/t8_gtest_mesh_handle.cxx b/test/mesh_handle/t8_gtest_mesh_handle.cxx index 32e0de7ecc..8bb9870b4b 100644 --- a/test/mesh_handle/t8_gtest_mesh_handle.cxx +++ b/test/mesh_handle/t8_gtest_mesh_handle.cxx @@ -29,6 +29,7 @@ along with t8code; if not, write to the Free Software Foundation, Inc., #include #include #include +#include "t8_gtest_common.hxx" #include #include @@ -196,4 +197,21 @@ TEST (t8_mesh_handle_test, test_union_element_competence_pack) EXPECT_FALSE (element_class::has_face_neighbor_cache ()); } +/** Check that the unique union of multiple mesh competence packs works as intended. */ +TEST (t8_mesh_handle_test, test_union_mesh_competence_pack) +{ + using namespace t8_mesh_handle; + using mesh_class + = mesh, + union_competence_packs_type, + data_mesh_competences_basic, empty_mesh_competences>>; + EXPECT_TRUE (mesh_class::has_element_data_handler_competence ()); + using element_class = typename mesh_class::element_class; + + EXPECT_TRUE (element_class::has_element_data_handler_competence ()); + EXPECT_TRUE (element_class::has_volume_cache ()); + EXPECT_TRUE (element_class::has_diameter_cache ()); +} + INSTANTIATE_TEST_SUITE_P (t8_gtest_mesh, t8_mesh_handle_test, testing::Combine (AllEclasses, testing::Range (2, 3))); diff --git a/tutorials/mesh_handle/t8_mesh_element_data.cxx b/tutorials/mesh_handle/t8_mesh_element_data.cxx index 03249f0002..7422d160e0 100644 --- a/tutorials/mesh_handle/t8_mesh_element_data.cxx +++ b/tutorials/mesh_handle/t8_mesh_element_data.cxx @@ -204,8 +204,8 @@ main (int argc, char **argv) { /* We put the mesh in its own scope so that it is automatically destroyed at the end of the scope. * This is only necessary because sc_finalize checks if there are leftover references. * This unique pointer would have been destroyed automatically at the end of the programme. */ - using mesh_class = t8_mesh_handle::mesh>; + using mesh_class = t8_mesh_handle::mesh>; auto mesh = build_mesh (comm, level); t8_mesh_handle::write_mesh_to_vtk (*mesh, prefix_mesh);