diff --git a/mesh_handle/README.md b/mesh_handle/README.md index 494420fcf1..92585ea5db 100644 --- a/mesh_handle/README.md +++ b/mesh_handle/README.md @@ -10,9 +10,8 @@ The folder's most important files are: - The [mesh.hxx](mesh.hxx) defines the mesh of the handle. This is the central file of the mesh handle. - The [constructor_wrappers.hxx](constructor_wrappers.hxx) allows to define a mesh handle using a cmesh instead of a forest and provides a very small number of examples where the user needs no cmesh. - The [element.hxx](element.hxx) defines the elements (mesh or ghost elements) of the mesh handle. -- The [competences.hxx](competences.hxx) defines additional competences/functionality of an element to access additional data. +- In the folder [competences/](competences/README.md), there are files that define additional competences/functionality of an element to access additional data. - The [competence_pack.hxx](competence_pack.hxx) is needed to pack element or mesh competences to pass it as template parameter to the mesh. -- The [data_handler.hxx](data_handler.hxx) provides competences to work with element data. - The [mesh_io.hxx](mesh_io.hxx) allows to output results in vtk format. Headers in the [internal/](internal/) folder are only intended to implement details of the mesh handle, so they should not need to be included for any other purpose. \ No newline at end of file diff --git a/mesh_handle/competence_pack.hxx b/mesh_handle/competence_pack.hxx index 34f030d8b3..d983d6819c 100644 --- a/mesh_handle/competence_pack.hxx +++ b/mesh_handle/competence_pack.hxx @@ -30,8 +30,9 @@ #pragma once -#include "competences.hxx" -#include "data_handler.hxx" +#include "competences/cache_element_competences.hxx" +#include "competences/element_data_competences.hxx" +#include "competences/dg_competences.hxx" #include "internal/competence_pack_union.hxx" namespace t8_mesh_handle { @@ -96,6 +97,9 @@ using empty_mesh_competences = mesh_competence_pack<>; template using data_mesh_competences = mesh_competence_pack::template type>; +/** Predefined mesh competence pack combining all competences that are useful for discontinuous Galerkin methods. */ +using dg_mesh_competences = mesh_competence_pack; + // --- Compute union of competence packs. --- /** Compute the unique union of the competences of several competence_pack. This could be * \ref t8_mesh_handle::element_competence_pack or \ref t8_mesh_handle::mesh_competence_pack. diff --git a/mesh_handle/competences/README.md b/mesh_handle/competences/README.md new file mode 100644 index 0000000000..c556568907 --- /dev/null +++ b/mesh_handle/competences/README.md @@ -0,0 +1,19 @@ +# Competences for the mesh handle # + +Definition of additional competences/functionalities that can be used for the t8_mesh_handle::mesh. + +- [element_data_competences.hxx](element_data_competences.hxx) provides competences to work with element data. +- [cache_element_competences.hxx](cache_element_competences.hxx) defines competences to cache functionalities of elements instead of calculating them each time a function is called. +- [dg_competences.hxx](dg_competences.hxx): competences that are useful for discontinuous Galerkin methods. + +### Note for developers or users that want to provide their own competence: + +All competences have the same inheritance pattern: +We use the CRTP pattern as we may need to access members of the derived classes like t8_mesh_handle::element. +The t8_crtp_operator is used for convenience/clear code (avoid to type a static cast explicitly each time we need functionality of TUnderlying). +Especially for the competences to cache functionality, the access of members is not necessary, such that it is not obvious why we use the crtp. +For competences that extend the functionality of the element (see e.g. [dg_competences.hxx](dg_competences.hxx)), this is required. +We use it for all competences for consistency and compatibility with t8_mesh_handle::element. + +We use t8_crtp_operator instead of t8_crtp_basic to circumvent diamond shaped inheritance pattern in the case that multiple competences are used together. +The distinction of the classes is made by the second template parameter of t8_crtp_operator. diff --git a/mesh_handle/competences.hxx b/mesh_handle/competences/cache_element_competences.hxx similarity index 99% rename from mesh_handle/competences.hxx rename to mesh_handle/competences/cache_element_competences.hxx index 9782149022..aaae5ff932 100644 --- a/mesh_handle/competences.hxx +++ b/mesh_handle/competences/cache_element_competences.hxx @@ -20,7 +20,7 @@ along with t8code; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. */ -/** \file competences.hxx +/** \file cache_element_competences.hxx * Definition of the additional competences/functionalities that can be used for \ref t8_mesh_handle::mesh. * Especially, competences to cache functionalities of elements instead of calculating them each time a function * is called are provided. diff --git a/mesh_handle/competences/dg_competences.hxx b/mesh_handle/competences/dg_competences.hxx new file mode 100644 index 0000000000..2cf20858f9 --- /dev/null +++ b/mesh_handle/competences/dg_competences.hxx @@ -0,0 +1,347 @@ +/* +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 dg_competences.hxx + * Mesh competences that are useful for discontinuous Galerkin methods. This includes the storage of the remote + * ranks of the ghost elements and the face connectivity. + */ + +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +namespace t8_mesh_handle +{ +/// Dummy value for the local rank. +inline constexpr int LOCAL_RANK = -1; + +/** + * This mesh competence provides the competence to store a vector with the ranks of the according elements. + * This is always \a LOCAL_RANK for the local elements. The interesting parts are the ghost elements. + * The order of the vector is the same as in the mesh handle. + * \tparam TUnderlying Use the \ref element with specified competences as template parameter. + */ +template +struct remote_ranks_mesh_competence: public t8_crtp_operator +{ + public: + /** + * Fill \a ranks with LOCAL_RANK for each local element and with the + * corresponding remote MPI rank for each ghost element, i.e.: + * ranks[0 .. num_local_elements - 1] = LOCAL_RANK + * ranks[num_local_elements .. num_local_elements + num_ghosts - 1] + * = remote rank of that ghost + */ + void + set_rank_vector () const + { + const TUnderlying& mesh = this->underlying (); + const t8_locidx_t num_local = mesh.get_num_local_elements (); + const t8_locidx_t num_ghosts = mesh.get_num_ghosts (); + // Nothing to do if all elements are local. + if (num_ghosts == 0) { + ranks.assign (num_local, LOCAL_RANK); + return; + } + + /* --- Local elements: rank = LOCAL_RANK --- */ + ranks.resize (num_local + num_ghosts); + std::fill_n (ranks.begin (), num_local, LOCAL_RANK); + + /* --- Ghost elements: determine rank per ghost element. --- */ + /* Get array with remote ranks in ascending order. */ + const t8_forest_t forest = mesh.get_forest (); + int num_remotes = 0; + int* remotes = t8_forest_ghost_get_remotes (forest, &num_remotes); + + for (int iremote = 0; iremote < num_remotes; ++iremote) { + const int remote = remotes[iremote]; + /* First ghost-element index that belongs to this remote rank. */ + const t8_locidx_t first_elem = t8_forest_ghost_remote_first_elem (forest, remote); + /* The number of ghost elements of this remote is the difference between first_elem and the next remote's first + * element index. */ + t8_locidx_t num_elems_of_remote; + if (iremote + 1 < num_remotes) { + const int next_remote = remotes[iremote + 1]; + const t8_locidx_t next_first = t8_forest_ghost_remote_first_elem (forest, next_remote); + num_elems_of_remote = next_first - first_elem; + } + else { + num_elems_of_remote = num_ghosts - first_elem; + } + + /* Write remote rank for every ghost element of this remote. */ + for (t8_locidx_t ielem = 0; ielem < num_elems_of_remote; ++ielem) { + ranks[num_local + first_elem + ielem] = remote; + } + } + } + + /** + * Get the rank associated with a given element handle. + * \param[in] element_handle_id The ID of the element handle. + * \return The rank of the specified element if \ref set_rank_vector was called beforehand. + */ + int + get_rank (t8_locidx_t element_handle_id) const + { + T8_ASSERT (static_cast (element_handle_id) < ranks.size ()); + return ranks[element_handle_id]; + } + + protected: + mutable std::vector ranks; ///< The rank of the owner for each element. +}; + +/** Type of the face. */ +enum class face_type { + BOUNDARY, ///< Exactly 1 side, domain boundary. + CONFORMAL, ///< Exactly 2 sides, same level, all local. + MORTAR, ///< 1 large side + N small sides, all local. + MPI_CONFORMAL, ///< Exactly 2 sides, same level, exactly one remote. + MPI_MORTAR, ///< Mortar where at least one side (large or small) is remote. +}; + +/** Class for the face side of an element. One \ref face can have multiple face sides of different elements. */ +struct face_side +{ + int element_id; ///< Mesh element handle id of the side's element. + int local_face_id; ///< The face of that element pointing to this interface. + int rank = LOCAL_RANK; ///< LOCAL_RANK if owned locally, else MPI rank of owner. +}; + +/** Class for a face. A face can have multiple \ref face_side s if different elements faces share the same face. + */ +struct face +{ + face_type type; ///< The face type of the face, see \ref face_type for the possible types and their meaning. + /** The face sides of different elements adjacent to this face. + * Order conventions for the sides in the vector depend on the face type: + * - BOUNDARY: sides[0] = the single local side + * - CONFORMAL / MPI_CONFORMAL: sides[0] = primary (smaller handle id); sides[1] = secondary. + * For MPI_CONFORMAL the local side is always the primary side with the smaller handle id (local ids < ghost ids). + * - MORTAR / MPI_MORTAR: sides[0] = large side; + * sides[1..N] = small sides (in face-corner order of the large element) + */ + std::vector sides; + int orientation = 0; ///< Face orientation code for coordinate permutation. +}; + +/** + * Mesh competence to build a unique vector of faces, where each face appears exactly once. + * This is useful for DG methods, where we need to loop over faces and access the neighboring elements across the face. + * This mesh competence should be combined with the competence \ref remote_ranks_mesh_competence. + * Each vector entry is of type \ref face and we have the following uniqueness rules for the face types specified in + * \ref face_type : + * - BOUNDARY: only one local side exists, always added. + * - CONFORMAL / MPI_CONFORMAL: added by the side whose element has the smaller handle_id. + * For a local–ghost pair the local element always has the smaller id (local ids < ghost ids), + * so the local side is always the one that inserts the face. + * - MORTAR / MPI_MORTAR: the large (coarser) side owns the face and inserts it (also for ghosts). + * The small sides are specified in sides. + * Additionally, a vector is build that holds the face indices for each element. + * + * \tparam TUnderlying Use the \ref mesh with specified competences as template parameter. + */ +template +struct face_vector_mesh_competence: public t8_crtp_operator +{ + public: + /** + * Build \a m_faces so that every face touching at least one local element appears exactly once. + * Simultaneously populate \a m_element_face_vector so that m_element_face_vector[handle_id][face_id] gives the + * index of the corresponding face in m_faces. + */ + void + set_unique_face_vector () const + { + if (!m_faces.empty () || !m_element_face_vector.empty ()) { + m_faces.clear (); + m_element_face_vector.clear (); + } + const TUnderlying& mesh = this->underlying (); + T8_ASSERT (mesh.is_balanced ()); + // Ensure that rank vector is set. + if constexpr (TUnderlying::has_remote_ranks_mesh_competence ()) { + mesh.set_rank_vector (); + } + else { + SC_ABORT ("Use remote_ranks_mesh_competence together with face_vector_mesh_competence.\n"); + } + + const t8_forest_t forest = mesh.get_forest (); + SC_CHECK_ABORT (forest->incomplete_trees == 0, "This functionality is not specified for incomplete trees.\n"); + + // Vector should store one entry per face of each element, so reserve space accordingly. + m_element_face_vector.assign (mesh.get_num_local_elements () + mesh.get_num_ghosts (), {}); + + // Check each face of each element and insert into vectors. + for (const auto& elem : mesh) { + const int handle_id = elem.get_element_handle_id (); + const int num_faces = elem.get_num_faces (); + if (m_element_face_vector[handle_id].empty ()) { + m_element_face_vector[handle_id].assign (num_faces, -1); + } + for (int iface = 0; iface < num_faces; ++iface) { + // Continue if the face entry is already filled. + if (m_element_face_vector[handle_id][iface] != -1) { + continue; + } + // Get neighbors and their dual face numbers to check the face_type. + std::vector dual_faces; + auto neighs = elem.get_face_neighbors (iface, dual_faces); + const int num_neighs = static_cast (neighs.size ()); + + if (num_neighs == 0) { + /* --- BOUNDARY --- */ + face f { face_type::BOUNDARY, { { handle_id, iface, LOCAL_RANK } }, /* orientation= */ 0 }; + const int face_idx = static_cast (m_faces.size ()); + m_faces.push_back (std::move (f)); + m_element_face_vector[handle_id][iface] = face_idx; + } + else if (num_neighs == 1) { + const int neigh_id = neighs[0]->get_element_handle_id (); + const int neigh_rank = mesh.get_rank (neigh_id); + + if (elem.get_level () == neighs[0]->get_level ()) { + /* --- CONFORMAL or MPI_CONFORMAL --- */ + // Insert only from the side with the smaller handle_id so that each conformal face appears exactly once. + if (handle_id < neigh_id) { + face f; + f.type = (neigh_rank != LOCAL_RANK) ? face_type::MPI_CONFORMAL : face_type::CONFORMAL; + f.sides.push_back ({ handle_id, iface, LOCAL_RANK }); + f.sides.push_back ({ neigh_id, dual_faces[0], neigh_rank }); + f.orientation = t8_forest_leaf_face_orientation ( + forest, elem.get_local_tree_id (), t8_forest_get_scheme (forest), elem.get_forest_element (), iface); + + const int face_idx = static_cast (m_faces.size ()); + m_faces.push_back (std::move (f)); + m_element_face_vector[handle_id][iface] = face_idx; + if (m_element_face_vector[neigh_id].empty ()) { + m_element_face_vector[neigh_id].assign (neighs[0]->get_num_faces (), -1); + } + m_element_face_vector[neigh_id][dual_faces[0]] = face_idx; + } + } + else if (elem.get_level () > neighs[0]->get_level ()) { + /* --- Small side of a MORTAR or MPI_MORTAR --- */ + /* The large (coarser) neighbour owns and inserts this face. + * We only need to do something here if the neighbor is remote. + */ + if (neigh_rank == LOCAL_RANK) { + continue; // local large side will insert the face, so we can skip this. + } + + // Check if the ghost is already inserted by another small mortar. If yes, only add the face. + if (!m_element_face_vector[neigh_id].empty () && m_element_face_vector[neigh_id][dual_faces[0]] != -1) { + m_faces[m_element_face_vector[neigh_id][dual_faces[0]]].sides.push_back ( + { handle_id, iface, LOCAL_RANK }); + m_element_face_vector[handle_id][iface] = m_element_face_vector[neigh_id][dual_faces[0]]; + continue; + } + + const int orientation = t8_forest_leaf_face_orientation ( + forest, elem.get_local_tree_id (), t8_forest_get_scheme (forest), elem.get_forest_element (), iface); + face f { face_type::MPI_MORTAR, + { { neigh_id, dual_faces[0], neigh_rank }, // Large mortar first. + { handle_id, iface, LOCAL_RANK } }, + orientation }; + + const int face_idx = static_cast (m_faces.size ()); + m_faces.push_back (std::move (f)); + m_element_face_vector[handle_id][iface] = face_idx; + if (m_element_face_vector[neigh_id].empty ()) { + m_element_face_vector[neigh_id].assign (neighs[0]->get_num_faces (), -1); + } + m_element_face_vector[neigh_id][dual_faces[0]] = face_idx; + } + // No other case possible without incomplete trees. + } + else { + /* --- Large side of a MORTAR or MPI_MORTAR (num_neighs > 1) --- */ + // We are the coarse element; neighs are all small mortars. + bool any_remote = std::any_of (neighs.begin (), neighs.end (), [&] (const auto* neigh) { + return mesh.get_rank (neigh->get_element_handle_id ()) != LOCAL_RANK; + }); + // Face index for the face we are about to insert. + const int face_idx = static_cast (m_faces.size ()); + m_element_face_vector[handle_id][iface] = face_idx; + + face f; + f.type = any_remote ? face_type::MPI_MORTAR : face_type::MORTAR; + f.orientation = t8_forest_leaf_face_orientation ( + forest, elem.get_local_tree_id (), t8_forest_get_scheme (forest), elem.get_forest_element (), iface); + f.sides.push_back ({ handle_id, iface, LOCAL_RANK }); + // Add small mortars to the sides vector and record face index for the small mortars in m_element_face_vector. + for (int ineigh = 0; ineigh < num_neighs; ++ineigh) { + const int neigh_id = neighs[ineigh]->get_element_handle_id (); + const int neigh_rank = mesh.get_rank (neigh_id); + f.sides.push_back ({ neigh_id, dual_faces[ineigh], neigh_rank }); + + if (m_element_face_vector[neigh_id].empty ()) { + m_element_face_vector[neigh_id].assign (neighs[ineigh]->get_num_faces (), -1); + } + m_element_face_vector[neigh_id][dual_faces[ineigh]] = face_idx; + } + // Insert face f. + m_faces.push_back (std::move (f)); + } + } + } + } + + /** Get the vector of unique faces. \ref set_unique_face_vector must be called beforehand to populate the vector. + * \return Const reference to the vector of faces. + */ + const std::vector& + get_unique_face_vector () const + { + T8_ASSERTF (!m_faces.empty (), "m_faces has not been set. Call set_unique_face_vector() first."); + return m_faces; + } + + /** Get the vector with indices of faces for each mesh handle element. + * \ref set_unique_face_vector must be called beforehand to populate the vector. + * \return Const reference to the element-face vector. + */ + const std::vector>& + get_element_face_vector () const + { + T8_ASSERTF (!m_element_face_vector.empty (), + " m_element_face_vector has not been set. Call set_unique_face_vector() first."); + return m_element_face_vector; + } + + protected: + /** Vector with one entry per unique face. */ + mutable std::vector m_faces; + /** For each element (local + ghost), the list of indices into m_faces + * that this element participates in. One entry per face of the element. */ + mutable std::vector> m_element_face_vector; +}; +} // namespace t8_mesh_handle diff --git a/mesh_handle/data_handler.hxx b/mesh_handle/competences/element_data_competences.hxx similarity index 99% rename from mesh_handle/data_handler.hxx rename to mesh_handle/competences/element_data_competences.hxx index b5ee9ce87b..3d52fda837 100644 --- a/mesh_handle/data_handler.hxx +++ b/mesh_handle/competences/element_data_competences.hxx @@ -20,7 +20,7 @@ 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. */ -/** \file data_handler.hxx +/** \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. @@ -28,7 +28,7 @@ #pragma once #include -#include "concepts.hxx" +#include #include #include #include diff --git a/mesh_handle/element.hxx b/mesh_handle/element.hxx index aa544b2899..9bdb9e803a 100644 --- a/mesh_handle/element.hxx +++ b/mesh_handle/element.hxx @@ -27,7 +27,7 @@ along with t8code; if not, write to the Free Software Foundation, Inc., #pragma once #include -#include "data_handler.hxx" +#include "competences/element_data_competences.hxx" #include #include #include @@ -44,7 +44,7 @@ namespace t8_mesh_handle * An element without specified template parameters provides default implementations for basic functionality * as accessing the refinement level or the centroid. With this implementation, the functionality is calculated each time * the function is called. - * Use the competences defined in \ref competences.hxx as template parameter to cache the functionality instead of + * Use the competences defined in \ref cache_element_competences.hxx as template parameter to cache the functionality instead of * recalculation in every function call. * To add functionality to the element, you can also simply write your own competence and give it as a template parameter. * You can access the functions implemented in your competence via the element. @@ -349,7 +349,6 @@ class element: public TCompetences>... { get_face_neighbors (int face, std::optional>> dual_faces = std::nullopt) const { T8_ASSERT ((face >= 0) && (face < get_num_faces ())); - SC_CHECK_ABORT (!m_is_ghost_element, "get_face_neighbors is not implemented for ghost elements.\n"); if constexpr (has_face_neighbor_cache ()) { if (this->neighbor_cache_filled (face)) { if (dual_faces) { @@ -519,6 +518,16 @@ class element: public TCompetences>... { return m_tree_id; } + /** Getter for the underlying forest element pointer. + * This function is mainly relevant for writing custom competences that need to access t8code functionality. + * \return Pointer to the underlying forest element. + */ + const t8_element_t* + get_forest_element () const + { + return m_element; + } + /** Getter for the local element id in the tree 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. * For mesh specific id use \ref get_element_handle_id. diff --git a/mesh_handle/mesh.hxx b/mesh_handle/mesh.hxx index 30a7de0df5..4c88c0a4bf 100644 --- a/mesh_handle/mesh.hxx +++ b/mesh_handle/mesh.hxx @@ -30,8 +30,8 @@ #include "element.hxx" #include "competence_pack.hxx" #include "internal/adapt.hxx" +#include "competences/element_data_competences.hxx" #include "concepts.hxx" -#include "data_handler.hxx" #include #include #include @@ -177,7 +177,7 @@ class mesh: public TMeshCompetencePack::template applyincomplete_trees == 0, "The mesh handle can't adapt forests with incomplete trees.\n"); if (!m_uncommitted_forest.has_value ()) { m_uncommitted_forest.emplace (); t8_forest_init (&*m_uncommitted_forest); @@ -420,6 +421,24 @@ class mesh: public TMeshCompetencePack::template apply #include -#include +#include #include #include #include diff --git a/test/mesh_handle/t8_gtest_custom_competence.cxx b/test/mesh_handle/competences/t8_gtest_custom_competence.cxx similarity index 98% rename from test/mesh_handle/t8_gtest_custom_competence.cxx rename to test/mesh_handle/competences/t8_gtest_custom_competence.cxx index 125e4ca014..e335c72cd1 100644 --- a/test/mesh_handle/t8_gtest_custom_competence.cxx +++ b/test/mesh_handle/competences/t8_gtest_custom_competence.cxx @@ -29,7 +29,7 @@ along with t8code; if not, write to the Free Software Foundation, Inc., #include #include -#include +#include #include #include #include diff --git a/test/mesh_handle/competences/t8_gtest_dg_competences.cxx b/test/mesh_handle/competences/t8_gtest_dg_competences.cxx new file mode 100644 index 0000000000..55d7c38f69 --- /dev/null +++ b/test/mesh_handle/competences/t8_gtest_dg_competences.cxx @@ -0,0 +1,226 @@ +/* +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_dg_competences.cxx + * Checks that the competences for discontinuous Galerkin methods defined in \ref dg_competences.hxx work as expected. + */ +#include +#include + +#include +#include +#include +#include + +/** Store the rank on each element. */ +struct 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. + */ +TEST (t8_gtest_dg_competences, remote_ranks) +{ + const int level = 2; + using namespace t8_mesh_handle; + using mesh_class + = mesh, + data_mesh_competences>>; + auto mesh = handle_hypercube_hybrid_uniform_default (level, sc_MPI_COMM_WORLD, true, true, false); + mesh->set_adapt (mesh_adapt_callback_test_refine_second); + mesh->set_partition (); + mesh->set_ghost (); + mesh->commit (); + + const t8_locidx_t num_local = mesh->get_num_local_elements (); + const t8_locidx_t num_ghosts = mesh->get_num_ghosts (); + if ((mesh->get_dimension () > 1) && (num_local > 1)) { + // Ensure that we actually test with ghost elements. + ASSERT_GT (num_ghosts, 0); + } + + int mpirank; + int mpiret = sc_MPI_Comm_rank (sc_MPI_COMM_WORLD, &mpirank); + SC_CHECK_MPI (mpiret); + + // Set local rank for all local mesh elements. + 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 (); + mesh->set_rank_vector (); + for (const auto& elem : *mesh) { + EXPECT_EQ (elem.get_element_data ().rank, mpirank); + EXPECT_EQ (LOCAL_RANK, mesh->get_rank (elem.get_element_handle_id ())); + } + for (t8_locidx_t ighost = num_local; ighost < num_local + num_ghosts; ighost++) { + EXPECT_EQ ((*mesh)[ighost].get_element_data ().rank, mesh->get_rank ((*mesh)[ighost].get_element_handle_id ())); + } +} + +/** Check the competence face_vector_mesh_competence for correctness. + * The test checks that the face vector and the element-face vector are consistent with each other and with the + * neighbors of the elements. + */ +TEST (t8_gtest_dg_competences, face_vector_mesh_competence) +{ + const int level = 2; + 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, true, false); + mesh->set_adapt (mesh_adapt_callback_test_refine_second); + mesh->set_partition (); + mesh->set_ghost (); + mesh->commit (); + + // Set the face vector using the competence function. + mesh->set_unique_face_vector (); + + const auto& faces = mesh->get_unique_face_vector (); + const auto& element_face_vector = mesh->get_element_face_vector (); + + const t8_locidx_t num_local = mesh->get_num_local_elements (); + const t8_locidx_t num_ghosts = mesh->get_num_ghosts (); + ASSERT_EQ (static_cast (element_face_vector.size ()), num_local + num_ghosts); + + // Check element_face_vector first. + for (t8_locidx_t ielem = num_local; ielem < num_local + num_ghosts; ielem++) { + EXPECT_EQ (static_cast (element_face_vector[ielem].size ()), (*mesh)[ielem].get_num_faces ()); + // Check that the element_face_vector points to valid face indices in the faces vector. + for (int iface = 0; iface < (*mesh)[ielem].get_num_faces (); ++iface) { + if (element_face_vector[ielem][iface] == -1) { + EXPECT_TRUE ((*mesh)[ielem].is_ghost_element ()); + } + else { + EXPECT_GE (element_face_vector[ielem][iface], 0); + EXPECT_LE (element_face_vector[ielem][iface], static_cast (faces.size ())); + } + } + } + + // Check faces vector. + for (int iface = 0; iface < static_cast (faces.size ()); ++iface) { + const auto& face = faces[iface]; + // Check that the face has at least one side and get the element and the neighbors of the first side. + EXPECT_GE (face.sides.size (), 1); + const auto elem_first = (*mesh)[face.sides[0].element_id]; + + std::vector dual_faces; + auto neighs = elem_first.get_face_neighbors (face.sides[0].local_face_id, dual_faces); + /* --- BOUNDARY --- */ + if (face.type == face_type::BOUNDARY) { + EXPECT_EQ (face.sides.size (), 1) << "BOUNDARY face must have exactly 1 side."; + EXPECT_EQ (face.sides[0].rank, LOCAL_RANK) << "BOUNDARY side must be local."; + EXPECT_FALSE (elem_first.is_ghost_element ()) << "BOUNDARY side element must be local."; + EXPECT_EQ (neighs.size (), 0) << "BOUNDARY side must not have neighbors."; + EXPECT_EQ (face.orientation, + t8_forest_leaf_face_orientation (mesh->get_forest (), elem_first.get_local_tree_id (), + t8_forest_get_scheme (mesh->get_forest ()), + elem_first.get_forest_element (), face.sides[0].local_face_id)) + << "BOUNDARY side face orientation must match the one computed from the forest."; + } + /* --- CONFORMAL --- */ + else if (face.type == face_type::CONFORMAL || face.type == face_type::MPI_CONFORMAL) { + EXPECT_EQ (face.sides.size (), 2) << "CONFORMAL face must have exactly 2 sides."; + EXPECT_EQ (face.sides[0].rank, LOCAL_RANK) << "CONFORMAL primary side must be local."; + if (face.type == face_type::CONFORMAL) { + EXPECT_EQ (face.sides[1].rank, LOCAL_RANK) << "CONFORMAL secondary side must be local."; + } + else { + EXPECT_EQ (face.sides[1].rank, mesh->get_rank (face.sides[1].element_id)) + << "MPI_CONFORMAL secondary side must be remote."; + } + EXPECT_EQ (neighs.size (), 1) << "CONFORMAL side must have exactly one neighbor."; + EXPECT_EQ (neighs[0]->get_element_handle_id (), face.sides[1].element_id) + << "CONFORMAL side neighbor must match the second side of the face."; + EXPECT_EQ (dual_faces[0], face.sides[1].local_face_id) << "CONFORMAL side neighbor face index must match."; + EXPECT_EQ (face.orientation, + t8_forest_leaf_face_orientation (mesh->get_forest (), neighs[0]->get_local_tree_id (), + t8_forest_get_scheme (mesh->get_forest ()), + neighs[0]->get_forest_element (), dual_faces[0])) + << "CONFORMAL side face orientation must match."; + } + /* --- MORTAR --- */ + else if (face.type == face_type::MORTAR || face.type == face_type::MPI_MORTAR) { + bool has_remote = std::any_of (face.sides.begin (), face.sides.end (), + [] (const face_side& s) { return s.rank != LOCAL_RANK; }); + EXPECT_EQ (has_remote, face.type == face_type::MPI_MORTAR); + + if (face.type == face_type::MORTAR) { + // Check that first element is the large mortar. + // For MPI_MORTARs, the large mortar could have only one neighbor. + EXPECT_GT (neighs.size (), 1) << "MORTAR face must have more than 2 sides."; + } + EXPECT_EQ (neighs[0]->get_level (), elem_first.get_level () + 1) + << "MORTAR first element must be the large mortar."; + EXPECT_EQ (neighs.size (), face.sides.size () - 1) + << "MORTAR side must have as many neighbors as small sides of the face."; + + for (int ineigh = 0; ineigh < static_cast (neighs.size ()); ++ineigh) { + // The neighbors must not have the same order as the face sides. + auto neigh_face_side + = std::find_if (face.sides.begin (), face.sides.end (), [&neighs, ineigh] (const face_side& s) { + return (neighs[ineigh]->get_element_handle_id () == s.element_id); + }); + EXPECT_FALSE (neigh_face_side == face.sides.end ()) << "MORTAR side neighbor must be found in the face sides."; + EXPECT_EQ (neighs[ineigh]->get_face_neighbors (dual_faces[ineigh]).size (), 1) + << "MORTAR side neighbor must have exactly one neighbor on this face."; + EXPECT_EQ (dual_faces[ineigh], (*neigh_face_side).local_face_id) + << "MORTAR side neighbor face index must match."; + EXPECT_EQ (face.orientation, + t8_forest_leaf_face_orientation (mesh->get_forest (), neighs[ineigh]->get_local_tree_id (), + t8_forest_get_scheme (mesh->get_forest ()), + neighs[ineigh]->get_forest_element (), dual_faces[ineigh])) + << "MORTAR side face orientation must match."; + } + } + + // This ensures that the element_face_vector is filled correctly. + for (const auto& side : face.sides) { + EXPECT_EQ (element_face_vector[side.element_id][side.local_face_id], iface); + } + } +} diff --git a/test/mesh_handle/t8_gtest_handle_data.cxx b/test/mesh_handle/competences/t8_gtest_handle_data.cxx similarity index 98% rename from test/mesh_handle/t8_gtest_handle_data.cxx rename to test/mesh_handle/competences/t8_gtest_handle_data.cxx index 851a1ae3f7..fe189c5116 100644 --- a/test/mesh_handle/t8_gtest_handle_data.cxx +++ b/test/mesh_handle/competences/t8_gtest_handle_data.cxx @@ -31,7 +31,7 @@ along with t8code; if not, write to the Free Software Foundation, Inc., #include #include #include -#include +#include #include #include #include diff --git a/test/mesh_handle/t8_gtest_adapt_partition_balance.cxx b/test/mesh_handle/t8_gtest_adapt_partition_balance.cxx index fe0da53078..fc301bbb6d 100644 --- a/test/mesh_handle/t8_gtest_adapt_partition_balance.cxx +++ b/test/mesh_handle/t8_gtest_adapt_partition_balance.cxx @@ -100,7 +100,7 @@ forest_adapt_callback_example (t8_forest_t forest, t8_forest_t forest_from, t8_l //--- 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_with_userdata. + * 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. diff --git a/test/mesh_handle/t8_gtest_ghost.cxx b/test/mesh_handle/t8_gtest_ghost.cxx index 919293dc5a..0923c525b2 100644 --- a/test/mesh_handle/t8_gtest_ghost.cxx +++ b/test/mesh_handle/t8_gtest_ghost.cxx @@ -30,7 +30,7 @@ along with t8code; if not, write to the Free Software Foundation, Inc., #include #include -#include +#include #include #include #include diff --git a/test/mesh_handle/t8_gtest_mesh_handle.cxx b/test/mesh_handle/t8_gtest_mesh_handle.cxx index 02abb6de0f..32e0de7ecc 100644 --- a/test/mesh_handle/t8_gtest_mesh_handle.cxx +++ b/test/mesh_handle/t8_gtest_mesh_handle.cxx @@ -31,7 +31,7 @@ along with t8code; if not, write to the Free Software Foundation, Inc., #include #include -#include +#include #include #include