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..8046d407ce 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 this element data directly for each element of the mesh. + * Also, 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,120 @@ 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 cannot store 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 (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 safer 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..4fa044789e 100644 --- a/mesh_handle/element.hxx +++ b/mesh_handle/element.hxx @@ -28,6 +28,7 @@ along with t8code; if not, write to the Free Software Foundation, Inc., #include #include "competences/element_data_competences.hxx" +#include "competences/cache_element_competences.hxx" #include #include #include @@ -36,6 +37,7 @@ along with t8code; if not, write to the Free Software Foundation, Inc., #include #include #include +#include namespace t8_mesh_handle { @@ -55,6 +57,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 +68,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). @@ -117,7 +119,7 @@ class element: public TCompetences>... { static constexpr bool has_volume_cache () { - return requires (SelfType& element) { element.volume_cache_filled (); }; + return has_competence (); } /** Function that checks if a cache for the element's diameter exists. @@ -126,7 +128,7 @@ class element: public TCompetences>... { static constexpr bool has_diameter_cache () { - return requires (SelfType& element) { element.diameter_cache_filled (); }; + return has_competence (); } /** Function that checks if a cache for the vertex coordinates exists. @@ -135,7 +137,7 @@ class element: public TCompetences>... { static constexpr bool has_vertex_cache () { - return requires (SelfType& element) { element.vertex_cache_filled (); }; + return has_competence (); } /** Function that checks if a cache for the centroid exists. @@ -144,7 +146,7 @@ class element: public TCompetences>... { static constexpr bool has_centroid_cache () { - return requires (SelfType& element) { element.centroid_cache_filled (); }; + return has_competence (); } /** Function that checks if a cache for the face neighbors exists. @@ -153,7 +155,7 @@ class element: public TCompetences>... { static constexpr bool has_face_neighbor_cache () { - return requires (SelfType& element) { element.neighbor_cache_filled (0); }; + return has_competence (); } /** Function that checks if a cache for the element's face areas exists. * \return true if a cache exists, false otherwise. @@ -161,7 +163,7 @@ class element: public TCompetences>... { static constexpr bool has_face_areas_cache () { - return requires (SelfType& element) { element.face_area_cache_filled (0); }; + return has_competence (); } /** Function that checks if a cache for the element's face centroids exists. @@ -170,7 +172,7 @@ class element: public TCompetences>... { static constexpr bool has_face_centroids_cache () { - return requires (SelfType& element) { element.face_centroid_cache_filled (0); }; + return has_competence (); } /** Function that checks if a cache for the element's face normals exists. @@ -179,16 +181,16 @@ class element: public TCompetences>... { static constexpr bool has_face_normals_cache () { - return requires (SelfType& element) { element.face_normal_cache_filled (0); }; + return has_competence (); } - /** Function that checks if a competence for element data handling is given to the element. + /** Function that checks if a competence for element-data handling is given to the element. * \return true if element has a data handler, false otherwise. */ static constexpr bool has_element_data_handler_competence () { - return requires (SelfType& element) { element.get_element_data (); }; + return has_competence (); } // --- Functionality of the element. In each function, it is checked if a cached version exists (and is used then). --- @@ -232,6 +234,14 @@ class element: public TCompetences>... { return t8_forest_get_scheme (m_mesh->m_forest)->element_get_shape (get_tree_class (), m_element); } + /** The number of children of the element when it is refined during adaption. + */ + int + get_num_children_per_refinement () const + { + return t8_forest_get_scheme (m_mesh->m_forest)->element_get_num_children (get_tree_class (), m_element); + } + /** Getter for the element's volume. * This is only an approximation. * This function uses or sets the cached version defined in TCompetence if available and calculates if not. @@ -516,7 +526,7 @@ class element: public TCompetences>... { #endif // --- Function to access mesh specific id. --- - /** Getter for the index of the element in the mesh to which the element belongs. + /** Getter for the local index of the element in the mesh to which the element belongs. * \return The local element id of the element in the mesh. */ t8_locidx_t @@ -531,6 +541,17 @@ class element: public TCompetences>... { return t8_forest_get_tree_element_offset (m_mesh->m_forest, m_tree_id) + m_element_id; } + /** Getter for the global index of the element in the mesh to which the element belongs. + * \return The global element id of the element in the mesh. + */ + t8_locidx_t + get_global_element_handle_id () const + { + SC_CHECK_ABORT (!m_is_ghost_element, "ERROR: Ghost elements do not have a global handle id.\n"); + return t8_forest_get_first_local_leaf_element_id (m_mesh->m_forest) + + t8_forest_get_tree_element_offset (m_mesh->m_forest, m_tree_id) + m_element_id; + } + //--- Getter for the member variables. --- /** Getter for the tree id of the element in the forest related to the mesh. * \warning This is related to t8code's tree structure and should not be confused with \ref mesh specific ids. @@ -600,6 +621,17 @@ class element: public TCompetences>... { { return t8_forest_get_tree_class (m_mesh->m_forest, m_tree_id); } + + /** Check whether the element was instantiated with a given competence. + * \tparam TCompetence The competence template. + * \return true if TCompetence is a base class of the element, false otherwise. + */ + template