From 8efb4e10fae6c38861954d838e90bcf36c11c74d Mon Sep 17 00:00:00 2001 From: Lena Ploetzke Date: Mon, 20 Apr 2026 10:27:38 +0200 Subject: [PATCH 01/12] Rename file --- mesh_handle/README.md | 2 +- mesh_handle/competence_pack.hxx | 2 +- .../cache_element_competences.hxx} | 2 +- mesh_handle/element.hxx | 2 +- test/mesh_handle/t8_gtest_cache_competence.cxx | 2 +- test/mesh_handle/t8_gtest_custom_competence.cxx | 2 +- test/mesh_handle/t8_gtest_ghost.cxx | 2 +- test/mesh_handle/t8_gtest_mesh_handle.cxx | 2 +- 8 files changed, 8 insertions(+), 8 deletions(-) rename mesh_handle/{competences.hxx => competences/cache_element_competences.hxx} (99%) diff --git a/mesh_handle/README.md b/mesh_handle/README.md index fc578bf32e..641fc47f93 100644 --- a/mesh_handle/README.md +++ b/mesh_handle/README.md @@ -10,7 +10,7 @@ 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/), 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 [mesh_io.hxx](mesh_io.hxx) allows to output results in vtk format. diff --git a/mesh_handle/competence_pack.hxx b/mesh_handle/competence_pack.hxx index d7406b5112..8b132e4520 100644 --- a/mesh_handle/competence_pack.hxx +++ b/mesh_handle/competence_pack.hxx @@ -26,7 +26,7 @@ #pragma once -#include "competences.hxx" +#include "competences/cache_element_competences.hxx" #include "internal/competence_pack_union.hxx" namespace t8_mesh_handle { 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/element.hxx b/mesh_handle/element.hxx index 78cb844f22..a715035542 100644 --- a/mesh_handle/element.hxx +++ b/mesh_handle/element.hxx @@ -43,7 +43,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. diff --git a/test/mesh_handle/t8_gtest_cache_competence.cxx b/test/mesh_handle/t8_gtest_cache_competence.cxx index 6c460e01b2..8aac637b34 100644 --- a/test/mesh_handle/t8_gtest_cache_competence.cxx +++ b/test/mesh_handle/t8_gtest_cache_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/t8_gtest_custom_competence.cxx b/test/mesh_handle/t8_gtest_custom_competence.cxx index 5db5f1f8a7..8bf1b65d0c 100644 --- a/test/mesh_handle/t8_gtest_custom_competence.cxx +++ b/test/mesh_handle/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/t8_gtest_ghost.cxx b/test/mesh_handle/t8_gtest_ghost.cxx index 7195a5b43e..04ad61dc97 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 d33eff9de5..426a3b2b81 100644 --- a/test/mesh_handle/t8_gtest_mesh_handle.cxx +++ b/test/mesh_handle/t8_gtest_mesh_handle.cxx @@ -30,7 +30,7 @@ along with t8code; if not, write to the Free Software Foundation, Inc., #include #include -#include +#include #include #include From 83395d138832655657b2cdb1c42e6df7628ab225 Mon Sep 17 00:00:00 2001 From: Lena Ploetzke Date: Mon, 20 Apr 2026 13:11:14 +0200 Subject: [PATCH 02/12] rename data handler --- mesh_handle/README.md | 3 +- mesh_handle/competence_pack.hxx | 2 +- mesh_handle/competences/README.md | 18 ++++++ mesh_handle/competences/dg_competences.hxx | 61 +++++++++++++++++++ .../element_data_competences.hxx} | 10 ++- mesh_handle/element.hxx | 2 +- mesh_handle/mesh.hxx | 2 +- .../t8_gtest_adapt_partition_balance.cxx | 2 +- test/mesh_handle/t8_gtest_handle_data.cxx | 2 +- .../mesh_handle/t8_mesh_element_data.cxx | 6 +- 10 files changed, 96 insertions(+), 12 deletions(-) create mode 100644 mesh_handle/competences/README.md create mode 100644 mesh_handle/competences/dg_competences.hxx rename mesh_handle/{data_handler.hxx => competences/element_data_competences.hxx} (93%) diff --git a/mesh_handle/README.md b/mesh_handle/README.md index 97d3d15b4d..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. -- In the folder [competences/](competences/), there are files that define 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 4dad3a89fe..310f555a97 100644 --- a/mesh_handle/competence_pack.hxx +++ b/mesh_handle/competence_pack.hxx @@ -31,7 +31,7 @@ #pragma once #include "competences/cache_element_competences.hxx" -#include "data_handler.hxx" +#include "competences/element_data_competences.hxx" #include "internal/competence_pack_union.hxx" namespace t8_mesh_handle { diff --git a/mesh_handle/competences/README.md b/mesh_handle/competences/README.md new file mode 100644 index 0000000000..39c668a5ad --- /dev/null +++ b/mesh_handle/competences/README.md @@ -0,0 +1,18 @@ +# 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. + +### 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, 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/dg_competences.hxx b/mesh_handle/competences/dg_competences.hxx new file mode 100644 index 0000000000..716aafd70e --- /dev/null +++ b/mesh_handle/competences/dg_competences.hxx @@ -0,0 +1,61 @@ +/* +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_element_competences.hxx + * TODO + */ + +// #pragma once + +// #include +// #include +// #include +// #include +// #include + +// namespace t8_mesh_handle +// { + +// /** +// * TODO +// * \tparam TUnderlying Use the \ref element with specified competences as template parameter. +// */ +// template +// struct face_vect: public t8_crtp_operator +// { +// public: +// /** +// * Function that checks if the cache for the volume has been filled. +// * \return true if the cache has been filled, false otherwise. +// */ +// bool +// volume_cache_filled () const +// { +// return m_volume.has_value (); +// } + +// protected: +// mutable std::optional +// m_volume; /**< Cache for the volume. Use optional to allow no value if cache is not filled. */ +// }; + +// } diff --git a/mesh_handle/data_handler.hxx b/mesh_handle/competences/element_data_competences.hxx similarity index 93% rename from mesh_handle/data_handler.hxx rename to mesh_handle/competences/element_data_competences.hxx index e8e18ef54a..66f01704d5 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. @@ -146,8 +146,12 @@ struct element_data_element_competence: public t8_crtp_basic { T8_ASSERT (this->underlying ().m_mesh->has_element_data_handler_competence ()); SC_CHECK_ABORT (!this->underlying ().is_ghost_element (), "Element data cannot be set for ghost elements.\n"); - this->underlying ().m_mesh->m_element_data.reserve (this->underlying ().m_mesh->get_num_local_elements () - + this->underlying ().m_mesh->get_num_ghosts ()); + if ((t8_locidx_t) this->underlying ().m_mesh->m_element_data.size () + <= this->underlying ().get_element_handle_id ()) { + this->underlying ().m_mesh->m_element_data.reserve (this->underlying ().m_mesh->get_num_local_elements () + + this->underlying ().m_mesh->get_num_ghosts ()); + this->underlying ().m_mesh->m_element_data.resize (this->underlying ().m_mesh->get_num_local_elements ()); + } this->underlying ().m_mesh->m_element_data[this->underlying ().get_element_handle_id ()] = std::move (element_data); } diff --git a/mesh_handle/element.hxx b/mesh_handle/element.hxx index 2d7f0a80df..fce2936187 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 diff --git a/mesh_handle/mesh.hxx b/mesh_handle/mesh.hxx index 65b96daf70..c063d965ed 100644 --- a/mesh_handle/mesh.hxx +++ b/mesh_handle/mesh.hxx @@ -30,7 +30,7 @@ #include "element.hxx" #include "competence_pack.hxx" #include "internal/adapt.hxx" -#include "data_handler.hxx" +#include "competences/element_data_competences.hxx" #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 fdab97fffa..b965662cb2 100644 --- a/test/mesh_handle/t8_gtest_adapt_partition_balance.cxx +++ b/test/mesh_handle/t8_gtest_adapt_partition_balance.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/t8_gtest_handle_data.cxx b/test/mesh_handle/t8_gtest_handle_data.cxx index 851a1ae3f7..fe189c5116 100644 --- a/test/mesh_handle/t8_gtest_handle_data.cxx +++ b/test/mesh_handle/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/tutorials/mesh_handle/t8_mesh_element_data.cxx b/tutorials/mesh_handle/t8_mesh_element_data.cxx index 44d1b5132b..fa0afc8d13 100644 --- a/tutorials/mesh_handle/t8_mesh_element_data.cxx +++ b/tutorials/mesh_handle/t8_mesh_element_data.cxx @@ -114,7 +114,8 @@ void set_element_data_mesh (TMeshClass &mesh) { for (auto &elem : mesh) { - elem.set_element_data ({ elem.get_level (), elem.get_volume () }); + data_per_element_type data ({ elem.get_level (), elem.get_volume () }); + elem.set_element_data (data); } } @@ -202,7 +203,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, data_per_element_type>; + using mesh_class = t8_mesh_handle::mesh>; auto mesh = build_mesh (comm, level); t8_mesh_handle::write_mesh_to_vtk (*mesh, prefix_mesh); From 6ad032d0f8c069fee5521a195ecd335b96c52eba Mon Sep 17 00:00:00 2001 From: Lena Ploetzke Date: Mon, 20 Apr 2026 16:00:22 +0200 Subject: [PATCH 03/12] ideas --- mesh_handle/competences/dg_competences.hxx | 113 ++++++++++++++------- 1 file changed, 78 insertions(+), 35 deletions(-) diff --git a/mesh_handle/competences/dg_competences.hxx b/mesh_handle/competences/dg_competences.hxx index 716aafd70e..8eccb2fb10 100644 --- a/mesh_handle/competences/dg_competences.hxx +++ b/mesh_handle/competences/dg_competences.hxx @@ -24,38 +24,81 @@ along with t8code; if not, write to the Free Software Foundation, Inc., * TODO */ -// #pragma once - -// #include -// #include -// #include -// #include -// #include - -// namespace t8_mesh_handle -// { - -// /** -// * TODO -// * \tparam TUnderlying Use the \ref element with specified competences as template parameter. -// */ -// template -// struct face_vect: public t8_crtp_operator -// { -// public: -// /** -// * Function that checks if the cache for the volume has been filled. -// * \return true if the cache has been filled, false otherwise. -// */ -// bool -// volume_cache_filled () const -// { -// return m_volume.has_value (); -// } - -// protected: -// mutable std::optional -// m_volume; /**< Cache for the volume. Use optional to allow no value if cache is not filled. */ -// }; - -// } +#pragma once + +#include +#include +#include +#include +#include + +namespace t8_mesh_handle +{ + +struct mortar_type +{ + static const int CONFORMAL = 0; + static const int SMALLMORTAR = 1; + static const int BIGMORTAR = 2; + static const int BOUNDARY = -1; +} struct face +{ + mortar_type m_mortartype; + std::vector m_element_ids; + std::optional ranks; +}; + +/** + * TODO + * \tparam TUnderlying Use the \ref element with specified competences as template parameter. + */ +template +struct face_vector_mesh_competence: public t8_crtp_operator +{ + public: + /** + * TODO + */ + void + set_unique_face_vector () const + { + for (const auto &elem : this->underlying ()) { + for (int iface = 0; iface < elem.get_num_faces (); ++iface) { + auto neighs = elem.get_face_neighbors (iface); + auto num_neighs = neighs.size (); + if (num_neighs == 0) { + // Face is a geometry boundary. + m_faces.push_back ({ mortar_type::BOUNDARY, { elem.get_element_handle_id () } }); + } + else if (num_neighs > 1) { + // Face is a big mortar. + std::vector temp; + temp.push_back (elem.get_element_handle_id ()); + for (int ineigh = 0; ineigh < num_neighs; ++ineigh) { + temp.push_back (neighs[ineigh]->get_element_handle_id ()); + } + m_faces.push_back ({ mortar_type::BIGMORTAR, temp }); + } + else { //num_neighs==1 + if (neighs[0]->get_level () == elem.get_level ()) { + // For conformal faces, the face should only be added once to the unique faces vector. + if (elem.get_element_handle_id () < neighs[0]->get_element_handle_id ()) { + m_faces.push_back ( + { mortar_type::CONFORMAL, { elem.get_element_handle_id (), neighs[0]->get_element_handle_id () } }); + } + } + else { + // Face is the small mortar. + m_faces.push_back ( + { mortar_type::SMALLMORTAR, { elem.get_element_handle_id (), neighs[0]->get_element_handle_id () } }); + } + } + } + } + } + + protected: + mutable std::vector m_faces; + mutable std::vector> m_elements_to_faces; +}; +} // namespace t8_mesh_handle From 5c209af475e1101ff92ea48b3377dd4f3b4e92fe Mon Sep 17 00:00:00 2001 From: Lena Ploetzke Date: Tue, 21 Apr 2026 16:55:20 +0200 Subject: [PATCH 04/12] changes to have vec with faces --- mesh_handle/competences/dg_competences.hxx | 51 +++++++++++++++++----- 1 file changed, 40 insertions(+), 11 deletions(-) diff --git a/mesh_handle/competences/dg_competences.hxx b/mesh_handle/competences/dg_competences.hxx index 8eccb2fb10..f754142904 100644 --- a/mesh_handle/competences/dg_competences.hxx +++ b/mesh_handle/competences/dg_competences.hxx @@ -34,18 +34,44 @@ along with t8code; if not, write to the Free Software Foundation, Inc., namespace t8_mesh_handle { +// Sentinel value for "this side is local" +inline constexpr int LOCAL_RANK = -1; -struct mortar_type +struct face_side { - static const int CONFORMAL = 0; - static const int SMALLMORTAR = 1; - static const int BIGMORTAR = 2; - static const int BOUNDARY = -1; -} struct face + int m_element_id; // local element id, or ghost layer index if remote + int m_local_face_id; // which face of that element points to this interface + int m_orientation; // face orientation code for coordinate permutation + int m_rank = LOCAL_RANK; // LOCAL_RANK if owned locally, else MPI rank of owner + bool + is_remote () const + { + return m_rank != LOCAL_RANK; + } +}; + +enum class face_type { + CONFORMAL, // exactly 2 sides, same level, all local + MORTAR, // 1 large side + N small sides (N=2 in 2D, up to 4 in 3D hex) + BOUNDARY, // exactly 1 side, domain boundary + MPI_CONFORMAL, // exactly 2 sides, same level, at least one remote + MPI_MORTAR, // mortar where at least one side is remote +}; + +struct face { - mortar_type m_mortartype; - std::vector m_element_ids; - std::optional ranks; + face_type m_type; + std::vector m_sides; + // Convention for MORTAR / MPI_MORTAR: + // m_sides[0] = large side + // m_sides[1..N] = small sides (in face-corner order of the large element) + // Convention for CONFORMAL / MPI_CONFORMAL: + // m_sides[0] = primary (lower global id or local owner) + // m_sides[1] = secondary + // Convention for BOUNDARY: + // m_sides[0] = the single local side + /** \note Currently unused but may be filled with t8_cmesh_set_attribute. */ + int m_boundary_tag = 0; // meaningful only for BOUNDARY // \note }; /** @@ -64,7 +90,8 @@ struct face_vector_mesh_competence: public t8_crtp_operatorunderlying ()) { for (int iface = 0; iface < elem.get_num_faces (); ++iface) { - auto neighs = elem.get_face_neighbors (iface); + std::vector dual_faces; + auto neighs = elem.get_face_neighbors (iface, dual_faces); auto num_neighs = neighs.size (); if (num_neighs == 0) { // Face is a geometry boundary. @@ -99,6 +126,8 @@ struct face_vector_mesh_competence: public t8_crtp_operator m_faces; - mutable std::vector> m_elements_to_faces; + // Inverse map: for each element, which face entries does it participate in? + // Store as (type, index_into_vector, side_role) for each local face + std::vector, MAX_FACES_PER_ELEM>> m_element_face_map; }; } // namespace t8_mesh_handle From cef3be589a018bbb92e9c741f222d89ff2242ae5 Mon Sep 17 00:00:00 2001 From: Lena Ploetzke Date: Tue, 5 May 2026 16:56:49 +0200 Subject: [PATCH 05/12] add remote rank functionality --- mesh_handle/competences/dg_competences.hxx | 93 ++++++++++++++++++++-- 1 file changed, 88 insertions(+), 5 deletions(-) diff --git a/mesh_handle/competences/dg_competences.hxx b/mesh_handle/competences/dg_competences.hxx index f754142904..c117231968 100644 --- a/mesh_handle/competences/dg_competences.hxx +++ b/mesh_handle/competences/dg_competences.hxx @@ -34,15 +34,99 @@ along with t8code; if not, write to the Free Software Foundation, Inc., namespace t8_mesh_handle { -// Sentinel value for "this side is local" +/// 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 m_ranks with LOCAL_RANK for each local element and with the + * corresponding remote MPI rank for each ghost element, i.e.: + * m_ranks[0 .. num_local_elements - 1] = LOCAL_RANK + * m_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 (); + t8_forest_t forest = mesh.get_forest (); + + const t8_locidx_t num_local = mesh.get_num_local_elements (); + const t8_locidx_t num_ghosts = mesh.get_num_ghosts (); + + m_ranks.resize (num_local + num_ghosts); + + /* --- local elements: rank = LOCAL_RANK --- */ + for (t8_locidx_t ilocal = 0; ilocal < num_local; ++ilocal) { + m_ranks[ilocal] = LOCAL_RANK; + } + + if (num_ghosts == 0) { + return; + } + + /* --- ghost elements: determine rank per ghost element --- */ + + /* t8_forest_ghost_get_remotes returns a plain int array of length + * num_remotes in ascending order, and fills num_remotes. */ + 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 (0-based within ghost elements only) + * 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 its first element index and the next remote's first + * element index (or the total ghost count for the last remote). */ + 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 the remote rank for every ghost element of this remote. + * The handle id of a ghost element is num_local + (ghost-only index). */ + for (t8_locidx_t ielem = 0; ielem < num_elems_of_remote; ++ielem) { + m_ranks[num_local + first_elem + ielem] = remote; + } + } + } + + int + get_rank (t8_locidx_t element_handle_id) + { + T8_ASSERT (element_handle_id < m_ranks.size ()); + return m_ranks[element_handle_id]; + } + + protected: + mutable std::vector m_ranks; +}; + struct face_side { int m_element_id; // local element id, or ghost layer index if remote int m_local_face_id; // which face of that element points to this interface int m_orientation; // face orientation code for coordinate permutation int m_rank = LOCAL_RANK; // LOCAL_RANK if owned locally, else MPI rank of owner + bool is_remote () const { @@ -71,7 +155,7 @@ struct face // Convention for BOUNDARY: // m_sides[0] = the single local side /** \note Currently unused but may be filled with t8_cmesh_set_attribute. */ - int m_boundary_tag = 0; // meaningful only for BOUNDARY // \note + int m_boundary_tag = 0; // meaningful only for BOUNDARY }; /** @@ -88,7 +172,7 @@ struct face_vector_mesh_competence: public t8_crtp_operatorunderlying ()) { + for (const auto& elem : this->underlying ()) { for (int iface = 0; iface < elem.get_num_faces (); ++iface) { std::vector dual_faces; auto neighs = elem.get_face_neighbors (iface, dual_faces); @@ -127,7 +211,6 @@ struct face_vector_mesh_competence: public t8_crtp_operator m_faces; // Inverse map: for each element, which face entries does it participate in? - // Store as (type, index_into_vector, side_role) for each local face - std::vector, MAX_FACES_PER_ELEM>> m_element_face_map; + std::vector> m_element_face_vector; }; } // namespace t8_mesh_handle From 7ede1dea9526acd704956f2024ab357e8cee8b54 Mon Sep 17 00:00:00 2001 From: Lena Ploetzke Date: Thu, 7 May 2026 14:48:33 +0200 Subject: [PATCH 06/12] new ideas --- mesh_handle/competences/dg_competences.hxx | 183 +++++++++++++++++---- 1 file changed, 149 insertions(+), 34 deletions(-) diff --git a/mesh_handle/competences/dg_competences.hxx b/mesh_handle/competences/dg_competences.hxx index c117231968..4fa4797605 100644 --- a/mesh_handle/competences/dg_competences.hxx +++ b/mesh_handle/competences/dg_competences.hxx @@ -21,7 +21,7 @@ along with t8code; if not, write to the Free Software Foundation, Inc., */ /** \file dg_element_competences.hxx - * TODO + * TODO: IDEA for test: write rank on element data and then check if this is correct for the ghosts */ #pragma once @@ -57,6 +57,7 @@ struct remote_ranks_mesh_competence: public t8_crtp_operatorunderlying (); t8_forest_t forest = mesh.get_forest (); @@ -160,57 +161,171 @@ struct face /** * TODO - * \tparam TUnderlying Use the \ref element with specified competences as template parameter. + * \tparam TUnderlying Use the \ref mesh with specified competences as template parameter. */ template struct face_vector_mesh_competence: public t8_crtp_operator { public: /** - * TODO + * 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] holds one index into m_faces + * per face of that element. + * + * Uniqueness rules + * ---------------- + * 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. + * The small sides are added to m_faces[*].m_sides by the + * large-side owner; they do not trigger their own insertion. */ void set_unique_face_vector () const { - for (const auto& elem : this->underlying ()) { - for (int iface = 0; iface < elem.get_num_faces (); ++iface) { - std::vector dual_faces; - auto neighs = elem.get_face_neighbors (iface, dual_faces); - auto num_neighs = neighs.size (); - if (num_neighs == 0) { - // Face is a geometry boundary. - m_faces.push_back ({ mortar_type::BOUNDARY, { elem.get_element_handle_id () } }); - } - else if (num_neighs > 1) { - // Face is a big mortar. - std::vector temp; - temp.push_back (elem.get_element_handle_id ()); - for (int ineigh = 0; ineigh < num_neighs; ++ineigh) { - temp.push_back (neighs[ineigh]->get_element_handle_id ()); + const TUnderlying& mesh = this->underlying (); + if constexpr (requires (TUnderlying& meshvar) { meshvar.set_rank_vector (); }) + mesh.set_rank_vector (); + } + else + { + SC_ABORT ("Use remote_ranks_mesh_competence together with face_vector_mesh_competence.\n"); + } + + /* One entry per element (local + ghost), each entry is a list of + * indices into m_faces. */ + m_element_face_vector.assign (mesh.get_num_local_elements () + mesh.get_num_ghosts (), {}); + + for (const auto& elem : mesh) { + const int handle_id = elem.get_element_handle_id (); + + for (int iface = 0; iface < elem.get_num_faces (); ++iface) { + 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; + f.m_type = face_type::BOUNDARY; + f.m_sides.push_back ({ handle_id, iface, /*orientation=*/0, LOCAL_RANK }); + + const int face_idx = static_cast (m_faces.size ()); + m_faces.push_back (std::move (f)); + m_element_face_vector[handle_id].push_back (face_idx); + } + else if (num_neighs == 1) { + const int neigh_handle = neighs[0]->get_element_handle_id (); + const int neigh_rank = mesh.get_rank (neigh_handle); + const bool neigh_is_remote = (neigh_rank != LOCAL_RANK); + + 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. For a local–ghost + * pair the local element always wins (local ids < ghost ids). */ + if (handle_id < neigh_handle) { + face f; + f.m_type = neigh_is_remote ? face_type::MPI_CONFORMAL : face_type::CONFORMAL; + f.m_sides.push_back ({ handle_id, iface, /*orientation=*/0, LOCAL_RANK }); + f.m_sides.push_back ({ neigh_handle, dual_faces[0], /*orientation=*/0, neigh_rank }); + + const int face_idx = static_cast (m_faces.size ()); + m_faces.push_back (std::move (f)); + m_element_face_vector[handle_id].push_back (face_idx); + m_element_face_vector[neigh_handle].push_back (face_idx); } - m_faces.push_back ({ mortar_type::BIGMORTAR, temp }); - } - else { //num_neighs==1 - if (neighs[0]->get_level () == elem.get_level ()) { - // For conformal faces, the face should only be added once to the unique faces vector. - if (elem.get_element_handle_id () < neighs[0]->get_element_handle_id ()) { - m_faces.push_back ( - { mortar_type::CONFORMAL, { elem.get_element_handle_id (), neighs[0]->get_element_handle_id () } }); + else { + /* The neighbour (local) will insert this face; just record + * that we will participate so the index can be filled in when + * the neighbour processes its face. Because we iterate in + * handle_id order the neighbour has already done so. */ + /* Find the face index already inserted by the neighbour. */ + for (int fi : m_element_face_vector[neigh_handle]) { + const face& candidate = m_faces[fi]; + if (candidate.m_sides.size () >= 2 + && ((candidate.m_sides[0].m_element_id == neigh_handle + && candidate.m_sides[1].m_element_id == handle_id) + || (candidate.m_sides[1].m_element_id == neigh_handle + && candidate.m_sides[0].m_element_id == handle_id))) { + m_element_face_vector[handle_id].push_back (fi); + break; + } } } - else { - // Face is the small mortar. - m_faces.push_back ( - { mortar_type::SMALLMORTAR, { elem.get_element_handle_id (), neighs[0]->get_element_handle_id () } }); + } + 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 are on a small side; do not insert, but record the + * face index once the large side has inserted it. + * Because we may reach the small side before the large side + * has been processed (it is a ghost), we defer: the large + * side's insertion loop will also update our entry. */ + /* Nothing to insert here; the large side handles it. */ + } + else { + /* ---- Large side of a MORTAR or MPI_MORTAR (num_neighs==1) ---- + * Our level < neighbour level: we are the coarse side. + * This branch is reached when a coarse local element borders a + * single finer element (possible in 2D with triangles). */ + face f; + f.m_type = neigh_is_remote ? face_type::MPI_MORTAR : face_type::MORTAR; + f.m_sides.push_back ({ handle_id, iface, /*orientation=*/0, LOCAL_RANK }); + f.m_sides.push_back ({ neigh_handle, dual_faces[0], /*orientation=*/0, neigh_rank }); + + const int face_idx = static_cast (m_faces.size ()); + m_faces.push_back (std::move (f)); + m_element_face_vector[handle_id].push_back (face_idx); + m_element_face_vector[neigh_handle].push_back (face_idx); + } + } + else { + /* ---- Large side of a MORTAR or MPI_MORTAR (num_neighs > 1) ---- + * We are the coarse element; neighs are all finer small sides. */ + bool any_remote = false; + for (int in = 0; in < num_neighs; ++in) { + if (neighs[in]->get_rank () != LOCAL_RANK) { + any_remote = true; + break; } } + + face f; + f.m_type = any_remote ? face_type::MPI_MORTAR : face_type::MORTAR; + /* m_sides[0] = large side (us) */ + f.m_sides.push_back ({ handle_id, iface, /*orientation=*/0, LOCAL_RANK }); + /* m_sides[1..N] = small sides in face-corner order */ + for (int in = 0; in < num_neighs; ++in) { + const int nh = neighs[in]->get_element_handle_id (); + const int nr = neighs[in]->get_rank (); + f.m_sides.push_back ({ nh, dual_faces[in], /*orientation=*/0, nr }); + } + + const int face_idx = static_cast (m_faces.size ()); + m_faces.push_back (std::move (f)); + + /* Record this face index for the large side and all small sides. */ + m_element_face_vector[handle_id].push_back (face_idx); + for (int in = 0; in < num_neighs; ++in) { + m_element_face_vector[neighs[in]->get_element_handle_id ()].push_back (face_idx); + } } } } +} - protected: - mutable std::vector m_faces; - // Inverse map: for each element, which face entries does it participate in? - std::vector> m_element_face_vector; +protected: + 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 From a8beefeee6e02f7cbe71012fb14ba84549d9db9d Mon Sep 17 00:00:00 2001 From: Lena Ploetzke Date: Mon, 11 May 2026 13:42:01 +0200 Subject: [PATCH 07/12] engineer function set unique face vector --- mesh_handle/competences/dg_competences.hxx | 359 ++++++++++----------- mesh_handle/mesh.hxx | 19 ++ 2 files changed, 189 insertions(+), 189 deletions(-) diff --git a/mesh_handle/competences/dg_competences.hxx b/mesh_handle/competences/dg_competences.hxx index 4fa4797605..0906f7422a 100644 --- a/mesh_handle/competences/dg_competences.hxx +++ b/mesh_handle/competences/dg_competences.hxx @@ -26,6 +26,8 @@ along with t8code; if not, write to the Free Software Foundation, Inc., #pragma once +#include +#include #include #include #include @@ -57,41 +59,31 @@ struct remote_ranks_mesh_competence: public t8_crtp_operatorunderlying (); - t8_forest_t forest = mesh.get_forest (); - const t8_locidx_t num_local = mesh.get_num_local_elements (); const t8_locidx_t num_ghosts = mesh.get_num_ghosts (); - - m_ranks.resize (num_local + num_ghosts); - - /* --- local elements: rank = LOCAL_RANK --- */ - for (t8_locidx_t ilocal = 0; ilocal < num_local; ++ilocal) { - m_ranks[ilocal] = LOCAL_RANK; - } - if (num_ghosts == 0) { + m_ranks.assign (num_local, LOCAL_RANK); return; } - /* --- ghost elements: determine rank per ghost element --- */ + /* --- local elements: rank = LOCAL_RANK --- */ + m_ranks.resize (num_local + num_ghosts); + std::fill_n (m_ranks.begin (), num_local, LOCAL_RANK); - /* t8_forest_ghost_get_remotes returns a plain int array of length - * num_remotes in ascending order, and fills num_remotes. */ + /* --- 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 (0-based within ghost elements only) - * that belongs to this remote rank. */ + /* 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 its first element index and the next remote's first - * element index (or the total ghost count for the last 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]; @@ -102,14 +94,17 @@ struct remote_ranks_mesh_competence: public t8_crtp_operator m_ranks; + mutable std::vector m_ranks; ///< The rank of the owner for each element. }; -struct face_side -{ - int m_element_id; // local element id, or ghost layer index if remote - int m_local_face_id; // which face of that element points to this interface - int m_orientation; // face orientation code for coordinate permutation - int m_rank = LOCAL_RANK; // LOCAL_RANK if owned locally, else MPI rank of owner - - bool - is_remote () const - { - return m_rank != LOCAL_RANK; - } +/** 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. }; -enum class face_type { - CONFORMAL, // exactly 2 sides, same level, all local - MORTAR, // 1 large side + N small sides (N=2 in 2D, up to 4 in 3D hex) - BOUNDARY, // exactly 1 side, domain boundary - MPI_CONFORMAL, // exactly 2 sides, same level, at least one remote - MPI_MORTAR, // mortar where at least one side is remote +/** Class for the face side of an element. One face can have multiple face sides of different elements. */ +struct face_side +{ + int m_element_id; ///< Mesh element handle id of the side's element. + int m_local_face_id; ///< The face of that element pointing to this interface. + int m_orientation; ///< Face orientation code for coordinate permutation. + int m_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 m_type; + face_type m_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: m_sides[0] = the single local side + * - MORTAR / MPI_MORTAR: m_sides[0] = large side; + * m_sides[1..N] = small sides (in face-corner order of the large element) + * - CONFORMAL / MPI_CONFORMAL: m_sides[0] = primary (smaller handle id); m_sides[1] = secondary + */ std::vector m_sides; - // Convention for MORTAR / MPI_MORTAR: - // m_sides[0] = large side - // m_sides[1..N] = small sides (in face-corner order of the large element) - // Convention for CONFORMAL / MPI_CONFORMAL: - // m_sides[0] = primary (lower global id or local owner) - // m_sides[1] = secondary - // Convention for BOUNDARY: - // m_sides[0] = the single local side - /** \note Currently unused but may be filled with t8_cmesh_set_attribute. */ - int m_boundary_tag = 0; // meaningful only for BOUNDARY + /** Boundary tag; meaningful only for BOUNDARY. + * \note Currently unused but may be filled with t8_cmesh_set_attribute. + */ + int m_boundary_tag = 0; }; /** - * TODO + * 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 m_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 @@ -168,164 +174,139 @@ struct face_vector_mesh_competence: public t8_crtp_operatorunderlying (); - if constexpr (requires (TUnderlying& meshvar) { meshvar.set_rank_vector (); }) + // Ensure that rank vector is set. + if constexpr (mesh.has_remote_ranks_mesh_competence ()) { mesh.set_rank_vector (); - } - else - { - SC_ABORT ("Use remote_ranks_mesh_competence together with face_vector_mesh_competence.\n"); - } - - /* One entry per element (local + ghost), each entry is a list of - * indices into m_faces. */ - m_element_face_vector.assign (mesh.get_num_local_elements () + mesh.get_num_ghosts (), {}); - - for (const auto& elem : mesh) { - const int handle_id = elem.get_element_handle_id (); - - for (int iface = 0; iface < elem.get_num_faces (); ++iface) { - 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; - f.m_type = face_type::BOUNDARY; - f.m_sides.push_back ({ handle_id, iface, /*orientation=*/0, LOCAL_RANK }); - - const int face_idx = static_cast (m_faces.size ()); - m_faces.push_back (std::move (f)); - m_element_face_vector[handle_id].push_back (face_idx); - } - else if (num_neighs == 1) { - const int neigh_handle = neighs[0]->get_element_handle_id (); - const int neigh_rank = mesh.get_rank (neigh_handle); - const bool neigh_is_remote = (neigh_rank != LOCAL_RANK); + } + 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, "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 (); + for (int iface = 0; iface < elem.get_num_faces (); ++iface) { + // 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; + f.m_type = face_type::BOUNDARY; + f.m_sides.push_back ({ handle_id, iface, /*orientation=*/0, LOCAL_RANK }); - 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. For a local–ghost - * pair the local element always wins (local ids < ghost ids). */ - if (handle_id < neigh_handle) { + const int face_idx = static_cast (m_faces.size ()); + m_faces.push_back (std::move (f)); + m_element_face_vector[handle_id].push_back (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.m_type = (neigh_rank != LOCAL_RANK) ? face_type::MPI_CONFORMAL : face_type::CONFORMAL; + const int orientation = t8_forest_leaf_face_orientation ( + forest, elem.get_tree_id (), t8_forest_get_scheme (forest), elem.get_element (), iface); + + f.m_sides.push_back ({ handle_id, iface, orientation, LOCAL_RANK }); + const int orientation_neigh + = t8_forest_leaf_face_orientation (forest, neighs[0]->get_tree_id (), t8_forest_get_scheme (forest), + neighs[0]->get_element (), dual_faces[0]); + f.m_sides.push_back ({ neigh_id, dual_faces[0], orientation_neigh, neigh_rank }); + + const int face_idx = static_cast (m_faces.size ()); + m_faces.push_back (std::move (f)); + m_element_face_vector[handle_id].push_back (face_idx); + m_element_face_vector[neigh_id].push_back (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 + } face f; - f.m_type = neigh_is_remote ? face_type::MPI_CONFORMAL : face_type::CONFORMAL; - f.m_sides.push_back ({ handle_id, iface, /*orientation=*/0, LOCAL_RANK }); - f.m_sides.push_back ({ neigh_handle, dual_faces[0], /*orientation=*/0, neigh_rank }); + f.m_type = face_type::MPI_MORTAR; + const int orientation_neigh + = t8_forest_leaf_face_orientation (forest, neighs[0]->get_tree_id (), t8_forest_get_scheme (forest), + neighs[0]->get_element (), dual_faces[0]); + f.m_sides.push_back ({ neigh_id, dual_faces[0], orientation_neigh, neigh_rank }); + const int orientation = t8_forest_leaf_face_orientation ( + forest, elem.get_tree_id (), t8_forest_get_scheme (forest), elem.get_element (), iface); + f.m_sides.push_back ({ handle_id, iface, orientation, LOCAL_RANK }); const int face_idx = static_cast (m_faces.size ()); m_faces.push_back (std::move (f)); m_element_face_vector[handle_id].push_back (face_idx); - m_element_face_vector[neigh_handle].push_back (face_idx); + m_element_face_vector[neigh_id].push_back (face_idx); } else { - /* The neighbour (local) will insert this face; just record - * that we will participate so the index can be filled in when - * the neighbour processes its face. Because we iterate in - * handle_id order the neighbour has already done so. */ - /* Find the face index already inserted by the neighbour. */ - for (int fi : m_element_face_vector[neigh_handle]) { - const face& candidate = m_faces[fi]; - if (candidate.m_sides.size () >= 2 - && ((candidate.m_sides[0].m_element_id == neigh_handle - && candidate.m_sides[1].m_element_id == handle_id) - || (candidate.m_sides[1].m_element_id == neigh_handle - && candidate.m_sides[0].m_element_id == handle_id))) { - m_element_face_vector[handle_id].push_back (fi); - break; - } - } + SC_ABORT ("Not possible.\n"); } } - 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 are on a small side; do not insert, but record the - * face index once the large side has inserted it. - * Because we may reach the small side before the large side - * has been processed (it is a ghost), we defer: the large - * side's insertion loop will also update our entry. */ - /* Nothing to insert here; the large side handles it. */ - } else { - /* ---- Large side of a MORTAR or MPI_MORTAR (num_neighs==1) ---- - * Our level < neighbour level: we are the coarse side. - * This branch is reached when a coarse local element borders a - * single finer element (possible in 2D with triangles). */ - face f; - f.m_type = neigh_is_remote ? face_type::MPI_MORTAR : face_type::MORTAR; - f.m_sides.push_back ({ handle_id, iface, /*orientation=*/0, LOCAL_RANK }); - f.m_sides.push_back ({ neigh_handle, dual_faces[0], /*orientation=*/0, neigh_rank }); - + /* --- Large side of a MORTAR or MPI_MORTAR (num_neighs > 1) --- */ + // We are the coarse element; neighs are all finer small sides. + bool any_remote = false; + for (int ineigh = 0; ineigh < num_neighs; ++in) { + if (mesh.get_rank (neighs[ineigh]->get_element_handle_id ()) != LOCAL_RANK) { + any_remote = true; + break; + } + } + // Face index for the face we are about to insert. const int face_idx = static_cast (m_faces.size ()); - m_faces.push_back (std::move (f)); m_element_face_vector[handle_id].push_back (face_idx); - m_element_face_vector[neigh_handle].push_back (face_idx); - } - } - else { - /* ---- Large side of a MORTAR or MPI_MORTAR (num_neighs > 1) ---- - * We are the coarse element; neighs are all finer small sides. */ - bool any_remote = false; - for (int in = 0; in < num_neighs; ++in) { - if (neighs[in]->get_rank () != LOCAL_RANK) { - any_remote = true; - break; - } - } - - face f; - f.m_type = any_remote ? face_type::MPI_MORTAR : face_type::MORTAR; - /* m_sides[0] = large side (us) */ - f.m_sides.push_back ({ handle_id, iface, /*orientation=*/0, LOCAL_RANK }); - /* m_sides[1..N] = small sides in face-corner order */ - for (int in = 0; in < num_neighs; ++in) { - const int nh = neighs[in]->get_element_handle_id (); - const int nr = neighs[in]->get_rank (); - f.m_sides.push_back ({ nh, dual_faces[in], /*orientation=*/0, nr }); - } - - const int face_idx = static_cast (m_faces.size ()); - m_faces.push_back (std::move (f)); - /* Record this face index for the large side and all small sides. */ - m_element_face_vector[handle_id].push_back (face_idx); - for (int in = 0; in < num_neighs; ++in) { - m_element_face_vector[neighs[in]->get_element_handle_id ()].push_back (face_idx); + face f; + f.m_type = any_remote ? face_type::MPI_MORTAR : face_type::MORTAR; + const int orientation = t8_forest_leaf_face_orientation ( + forest, elem.get_tree_id (), t8_forest_get_scheme (forest), elem.get_element (), iface); + f.m_sides.push_back ({ handle_id, iface, orientation, LOCAL_RANK }); + 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); + const int orientation_neigh + = t8_forest_leaf_face_orientation (forest, neighs[ineigh]->get_tree_id (), t8_forest_get_scheme (forest), + neighs[0]->get_element (), dual_faces[ineigh]); + f.m_sides.push_back ({ neigh_id, dual_faces[ineigh], orientation_neigh, neigh_rank }); + // Record face index for the small side element. + m_element_face_vector[neighs[in]->get_element_handle_id ()].push_back (face_idx); + } + // Insert face f. + m_faces.push_back (std::move (f)); } } } } -} -protected: - 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; + 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/mesh.hxx b/mesh_handle/mesh.hxx index c063d965ed..0a984dbf62 100644 --- a/mesh_handle/mesh.hxx +++ b/mesh_handle/mesh.hxx @@ -301,6 +301,7 @@ class mesh: public TMeshCompetencePack::template applyincomplete_trees, "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); @@ -418,6 +419,24 @@ class mesh: public TMeshCompetencePack::template apply Date: Mon, 11 May 2026 16:13:56 +0200 Subject: [PATCH 08/12] add competences folder to test --- mesh_handle/competence_pack.hxx | 4 +++ mesh_handle/competences/dg_competences.hxx | 33 ++++++++++++++++--- mesh_handle/mesh.hxx | 2 +- test/CMakeLists.txt | 7 ++-- test/mesh_handle/README.md | 8 +++-- .../t8_gtest_cache_element_competence.cxx} | 0 .../t8_gtest_custom_competence.cxx | 0 .../competences/t8_gtest_dg_competences.cxx | 0 .../t8_gtest_handle_data.cxx | 0 9 files changed, 42 insertions(+), 12 deletions(-) rename test/mesh_handle/{t8_gtest_cache_competence.cxx => competences/t8_gtest_cache_element_competence.cxx} (100%) rename test/mesh_handle/{ => competences}/t8_gtest_custom_competence.cxx (100%) create mode 100644 test/mesh_handle/competences/t8_gtest_dg_competences.cxx rename test/mesh_handle/{ => competences}/t8_gtest_handle_data.cxx (100%) diff --git a/mesh_handle/competence_pack.hxx b/mesh_handle/competence_pack.hxx index 310f555a97..d983d6819c 100644 --- a/mesh_handle/competence_pack.hxx +++ b/mesh_handle/competence_pack.hxx @@ -32,6 +32,7 @@ #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/dg_competences.hxx b/mesh_handle/competences/dg_competences.hxx index 0906f7422a..553862579f 100644 --- a/mesh_handle/competences/dg_competences.hxx +++ b/mesh_handle/competences/dg_competences.hxx @@ -21,6 +21,8 @@ along with t8code; if not, write to the Free Software Foundation, Inc., */ /** \file dg_element_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. * TODO: IDEA for test: write rank on element data and then check if this is correct for the ghosts */ @@ -59,7 +61,6 @@ struct remote_ranks_mesh_competence: public t8_crtp_operatorunderlying (); const t8_locidx_t num_local = mesh.get_num_local_elements (); const t8_locidx_t num_ghosts = mesh.get_num_ghosts (); @@ -125,7 +126,7 @@ enum class face_type { MPI_MORTAR, ///< Mortar where at least one side (large or small) is remote. }; -/** Class for the face side of an element. One face can have multiple face sides of different elements. */ +/** Class for the face side of an element. One \ref face can have multiple face sides of different elements. */ struct face_side { int m_element_id; ///< Mesh element handle id of the side's element. @@ -190,7 +191,7 @@ struct face_vector_mesh_competence: public t8_crtp_operatorincomplete_trees, "This functionality is not specified for incomplete trees.\n"); + 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 (), {}); @@ -270,7 +271,7 @@ struct face_vector_mesh_competence: public t8_crtp_operator 1) --- */ // We are the coarse element; neighs are all finer small sides. bool any_remote = false; - for (int ineigh = 0; ineigh < num_neighs; ++in) { + for (int ineigh = 0; ineigh < num_neighs; ++ineigh) { if (mesh.get_rank (neighs[ineigh]->get_element_handle_id ()) != LOCAL_RANK) { any_remote = true; break; @@ -293,7 +294,7 @@ struct face_vector_mesh_competence: public t8_crtp_operatorget_element (), dual_faces[ineigh]); f.m_sides.push_back ({ neigh_id, dual_faces[ineigh], orientation_neigh, neigh_rank }); // Record face index for the small side element. - m_element_face_vector[neighs[in]->get_element_handle_id ()].push_back (face_idx); + m_element_face_vector[neighs[ineigh]->get_element_handle_id ()].push_back (face_idx); } // Insert face f. m_faces.push_back (std::move (f)); @@ -302,6 +303,28 @@ struct face_vector_mesh_competence: public t8_crtp_operator& + 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; diff --git a/mesh_handle/mesh.hxx b/mesh_handle/mesh.hxx index 0a984dbf62..4c84076acc 100644 --- a/mesh_handle/mesh.hxx +++ b/mesh_handle/mesh.hxx @@ -301,7 +301,7 @@ class mesh: public TMeshCompetencePack::template applyincomplete_trees, "The mesh handle can't adapt forests with incomplete trees.\n"); + SC_CHECK_ABORT (m_forest->incomplete_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); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 7d4ad406e6..686347342c 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -223,8 +223,9 @@ if( T8CODE_BUILD_MESH_HANDLE ) add_t8_cpp_test( NAME t8_gtest_mesh_handle_parallel SOURCES mesh_handle/t8_gtest_mesh_handle.cxx ) add_t8_cpp_test( NAME t8_gtest_compare_handle_to_forest_serial SOURCES mesh_handle/t8_gtest_compare_handle_to_forest.cxx ) add_t8_cpp_test( NAME t8_gtest_handle_ghost_parallel SOURCES mesh_handle/t8_gtest_ghost.cxx ) - add_t8_cpp_test( NAME t8_gtest_custom_competence_serial SOURCES mesh_handle/t8_gtest_custom_competence.cxx ) - add_t8_cpp_test( NAME t8_gtest_cache_competence_serial SOURCES mesh_handle/t8_gtest_cache_competence.cxx ) - add_t8_cpp_test( NAME t8_gtest_handle_data_parallel SOURCES mesh_handle/t8_gtest_handle_data.cxx ) add_t8_cpp_test( NAME t8_gtest_adapt_partition_balance_parallel SOURCES mesh_handle/t8_gtest_adapt_partition_balance.cxx ) + + add_t8_cpp_test( NAME t8_gtest_custom_competence_serial SOURCES mesh_handle/competences/t8_gtest_custom_competence.cxx ) + add_t8_cpp_test( NAME t8_gtest_cache_competence_serial SOURCES mesh_handle/competences/t8_gtest_cache_element_competence.cxx ) + add_t8_cpp_test( NAME t8_gtest_handle_data_parallel SOURCES mesh_handle/competences/t8_gtest_handle_data.cxx ) endif() diff --git a/test/mesh_handle/README.md b/test/mesh_handle/README.md index e76d668e68..e48ca10a2a 100644 --- a/test/mesh_handle/README.md +++ b/test/mesh_handle/README.md @@ -6,9 +6,11 @@ The testing strategy is written down here to ensure consistency. Therefore, the tests are structured as follows: - [t8_gtest_mesh_handle.cxx](t8_gtest_mesh_handle.cxx) tests some basic functionality, e.g., if the mesh handle generally works with default competences. The handle's functionality is not tested in detail, but instead using some *exemplary functions*. -- [t8_gtest_custom_competence.cxx](t8_gtest_custom_competence.cxx) checks that user defined competences can be used for the mesh handle. - [t8_gtest_compare_handle_to_forest.cxx](t8_gtest_compare_handle_to_forest.cxx) tests that the functionality of the handle gives the same results as if worked with the forest directly. This test checks *all functions*. -- [t8_gtest_cache_competence.cxx](t8_gtest_cache_competence.cxx) tests that *all predefined caching competences* work as intended. - [t8_gtest_ghost.cxx](t8_gtest_ghost.cxx) checks ghosts and the neighbor algorithm. Furthermore tests if *all functions work also for ghost cells* if applicable. -- [t8_gtest_handle_data.cxx](t8_gtest_handle_data.cxx) tests that user and element data functionality works as intended. - [t8_gtest_adapt_partition_balance.cxx](t8_gtest_adapt_partition_balance.cxx) implements tests for the adapt, partition and balance routines of mesh handle. + +The folder [competences/](competences/) checks element and mesh competences of the mesh handle. We have the following tests: +- [t8_gtest_custom_competence.cxx](t8_gtest_custom_competence.cxx) checks that user defined competences can be used for the mesh handle. +- [t8_gtest_cache_competence.cxx](t8_gtest_cache_competence.cxx) tests that *all predefined caching competences* work as intended. +- [t8_gtest_handle_data.cxx](t8_gtest_handle_data.cxx) tests that user and element data functionality works as intended. The element data functionality is provided via a competence. \ No newline at end of file diff --git a/test/mesh_handle/t8_gtest_cache_competence.cxx b/test/mesh_handle/competences/t8_gtest_cache_element_competence.cxx similarity index 100% rename from test/mesh_handle/t8_gtest_cache_competence.cxx rename to test/mesh_handle/competences/t8_gtest_cache_element_competence.cxx diff --git a/test/mesh_handle/t8_gtest_custom_competence.cxx b/test/mesh_handle/competences/t8_gtest_custom_competence.cxx similarity index 100% rename from test/mesh_handle/t8_gtest_custom_competence.cxx rename to test/mesh_handle/competences/t8_gtest_custom_competence.cxx 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..e69de29bb2 diff --git a/test/mesh_handle/t8_gtest_handle_data.cxx b/test/mesh_handle/competences/t8_gtest_handle_data.cxx similarity index 100% rename from test/mesh_handle/t8_gtest_handle_data.cxx rename to test/mesh_handle/competences/t8_gtest_handle_data.cxx From 1deb0a854458e6eac06fac5e7594d59f3454bef0 Mon Sep 17 00:00:00 2001 From: Lena Ploetzke Date: Tue, 12 May 2026 14:43:47 +0200 Subject: [PATCH 09/12] First working version of a test --- mesh_handle/competences/dg_competences.hxx | 105 ++++----- mesh_handle/element.hxx | 13 +- test/CMakeLists.txt | 1 + .../competences/t8_gtest_dg_competences.cxx | 216 ++++++++++++++++++ 4 files changed, 281 insertions(+), 54 deletions(-) diff --git a/mesh_handle/competences/dg_competences.hxx b/mesh_handle/competences/dg_competences.hxx index 553862579f..c4ab11086e 100644 --- a/mesh_handle/competences/dg_competences.hxx +++ b/mesh_handle/competences/dg_competences.hxx @@ -52,10 +52,10 @@ struct remote_ranks_mesh_competence: public t8_crtp_operator (element_handle_id) < ranks.size ()); + return ranks[element_handle_id]; } protected: - mutable std::vector m_ranks; ///< The rank of the owner for each element. + mutable std::vector ranks; ///< The rank of the owner for each element. }; /** Type of the face. */ @@ -129,29 +129,26 @@ enum class face_type { /** Class for the face side of an element. One \ref face can have multiple face sides of different elements. */ struct face_side { - int m_element_id; ///< Mesh element handle id of the side's element. - int m_local_face_id; ///< The face of that element pointing to this interface. - int m_orientation; ///< Face orientation code for coordinate permutation. - int m_rank = LOCAL_RANK; ///< LOCAL_RANK if owned locally, else MPI rank of owner. + 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 orientation = 0; ///< Face orientation code for coordinate permutation. + 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 m_type; ///< The face type of the face, see \ref face_type for the possible types and their meaning. + 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: m_sides[0] = the single local side - * - MORTAR / MPI_MORTAR: m_sides[0] = large side; - * m_sides[1..N] = small sides (in face-corner order of the large element) - * - CONFORMAL / MPI_CONFORMAL: m_sides[0] = primary (smaller handle id); m_sides[1] = secondary + * - 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 m_sides; - /** Boundary tag; meaningful only for BOUNDARY. - * \note Currently unused but may be filled with t8_cmesh_set_attribute. - */ - int m_boundary_tag = 0; + std::vector sides; }; /** @@ -165,7 +162,7 @@ struct face * 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 m_sides. + * 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. @@ -182,6 +179,10 @@ struct face_vector_mesh_competence: public t8_crtp_operatorunderlying (); // Ensure that rank vector is set. if constexpr (mesh.has_remote_ranks_mesh_competence ()) { @@ -208,8 +209,8 @@ struct face_vector_mesh_competence: public t8_crtp_operator (m_faces.size ()); m_faces.push_back (std::move (f)); @@ -224,15 +225,15 @@ struct face_vector_mesh_competence: public t8_crtp_operatorget_tree_id (), t8_forest_get_scheme (forest), - neighs[0]->get_element (), dual_faces[0]); - f.m_sides.push_back ({ neigh_id, dual_faces[0], orientation_neigh, neigh_rank }); + f.sides.push_back ({ handle_id, iface, orientation, LOCAL_RANK }); + const int orientation_neigh = t8_forest_leaf_face_orientation ( + forest, neighs[0]->get_local_tree_id (), t8_forest_get_scheme (forest), + neighs[0]->get_forest_element (), dual_faces[0]); + f.sides.push_back ({ neigh_id, dual_faces[0], orientation_neigh, neigh_rank }); const int face_idx = static_cast (m_faces.size ()); m_faces.push_back (std::move (f)); @@ -240,7 +241,7 @@ struct face_vector_mesh_competence: public t8_crtp_operatorget_level ()) { + 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. @@ -249,14 +250,14 @@ struct face_vector_mesh_competence: public t8_crtp_operatorget_tree_id (), t8_forest_get_scheme (forest), - neighs[0]->get_element (), dual_faces[0]); - f.m_sides.push_back ({ neigh_id, dual_faces[0], orientation_neigh, neigh_rank }); + = t8_forest_leaf_face_orientation (forest, neighs[0]->get_local_tree_id (), t8_forest_get_scheme (forest), + neighs[0]->get_forest_element (), dual_faces[0]); + f.sides.push_back ({ neigh_id, dual_faces[0], orientation_neigh, neigh_rank }); const int orientation = t8_forest_leaf_face_orientation ( - forest, elem.get_tree_id (), t8_forest_get_scheme (forest), elem.get_element (), iface); - f.m_sides.push_back ({ handle_id, iface, orientation, LOCAL_RANK }); + forest, elem.get_local_tree_id (), t8_forest_get_scheme (forest), elem.get_forest_element (), iface); + f.sides.push_back ({ handle_id, iface, orientation, LOCAL_RANK }); const int face_idx = static_cast (m_faces.size ()); m_faces.push_back (std::move (f)); @@ -264,7 +265,7 @@ struct face_vector_mesh_competence: public t8_crtp_operatorget_element_handle_id (); const int neigh_rank = mesh.get_rank (neigh_id); - const int orientation_neigh - = t8_forest_leaf_face_orientation (forest, neighs[ineigh]->get_tree_id (), t8_forest_get_scheme (forest), - neighs[0]->get_element (), dual_faces[ineigh]); - f.m_sides.push_back ({ neigh_id, dual_faces[ineigh], orientation_neigh, neigh_rank }); + const int orientation_neigh = t8_forest_leaf_face_orientation ( + forest, neighs[ineigh]->get_local_tree_id (), t8_forest_get_scheme (forest), + neighs[ineigh]->get_forest_element (), dual_faces[ineigh]); + f.sides.push_back ({ neigh_id, dual_faces[ineigh], orientation_neigh, neigh_rank }); // Record face index for the small side element. m_element_face_vector[neighs[ineigh]->get_element_handle_id ()].push_back (face_idx); } @@ -309,7 +310,7 @@ struct face_vector_mesh_competence: public t8_crtp_operator& get_unique_face_vector () const { - T8_ASSERTF (m_faces.empty (), "m_faces has not been set. Call set_unique_face_vector() first."); + T8_ASSERTF (!m_faces.empty (), "m_faces has not been set. Call set_unique_face_vector() first."); return m_faces; } @@ -320,7 +321,7 @@ struct face_vector_mesh_competence: public t8_crtp_operator>& get_element_face_vector () const { - T8_ASSERTF (m_element_face_vector.empty (), + 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; } diff --git a/mesh_handle/element.hxx b/mesh_handle/element.hxx index fce2936187..a04b24e840 100644 --- a/mesh_handle/element.hxx +++ b/mesh_handle/element.hxx @@ -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) { @@ -382,7 +381,7 @@ class element: public TCompetences>... { } if (num_neighbors > 0) { // Free allocated memory. - t8_forest_get_scheme (m_mesh->m_forest)->element_destroy (get_tree_class (), num_neighbors, neighbors); + t8_forest_get_scheme (m_mesh->m_forest)->element_destroy (neigh_class, num_neighbors, neighbors); T8_FREE (neighbors); T8_FREE (dual_faces_internal); T8_FREE (neighids); @@ -521,6 +520,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/test/CMakeLists.txt b/test/CMakeLists.txt index 686347342c..c54437a2d4 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -228,4 +228,5 @@ if( T8CODE_BUILD_MESH_HANDLE ) add_t8_cpp_test( NAME t8_gtest_custom_competence_serial SOURCES mesh_handle/competences/t8_gtest_custom_competence.cxx ) add_t8_cpp_test( NAME t8_gtest_cache_competence_serial SOURCES mesh_handle/competences/t8_gtest_cache_element_competence.cxx ) add_t8_cpp_test( NAME t8_gtest_handle_data_parallel SOURCES mesh_handle/competences/t8_gtest_handle_data.cxx ) + add_t8_cpp_test( NAME t8_gtest_dg_competences_parallel SOURCES mesh_handle/competences/t8_gtest_dg_competences.cxx ) endif() diff --git a/test/mesh_handle/competences/t8_gtest_dg_competences.cxx b/test/mesh_handle/competences/t8_gtest_dg_competences.cxx index e69de29bb2..b02c367376 100644 --- a/test/mesh_handle/competences/t8_gtest_dg_competences.cxx +++ b/test/mesh_handle/competences/t8_gtest_dg_competences.cxx @@ -0,0 +1,216 @@ +/* +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 work as expected. + */ +#include +#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_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. + * \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; +} + +/** TODO + */ +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 (); + + if ((mesh->get_dimension () > 1) && (mesh->get_num_local_elements () > 1)) { + // Ensure that we actually test with ghost elements. + ASSERT_GT (mesh->get_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 (mesh->get_num_local_elements (), { mpirank }); + mesh->set_element_data (std::move (element_data)); + // Get element data and check that the remote ranks functionality 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 = mesh->get_num_local_elements (); + ighost < mesh->get_num_local_elements () + mesh->get_num_ghosts (); ighost++) { + EXPECT_EQ ((*mesh)[ighost].get_element_data ().rank, mesh->get_rank ((*mesh)[ighost].get_element_handle_id ())); + } +} + +TEST (t8_gtest_dg_competences, face_vector_mesh_competence) +{ + const int level = 1; + 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 (); + + if ((mesh->get_dimension () > 1) && (mesh->get_num_local_elements () > 1)) { + ASSERT_GT (mesh->get_num_ghosts (), 0); + } + // 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 (const auto& elem : *mesh) { + const int handle_id = elem.get_element_handle_id (); + EXPECT_EQ (static_cast (element_face_vector[handle_id].size ()), elem.get_num_faces ()) + << "Element " << handle_id << " has wrong number of face indices."; + // Check that the element_face_vector points to valid face indices in the faces vector. + for (int iface = 0; iface < elem.get_num_faces (); ++iface) { + EXPECT_GE (element_face_vector[handle_id][iface], 0); + EXPECT_LE (element_face_vector[handle_id][iface], static_cast (faces.size ())); + } + } + for (t8_locidx_t ighost = num_local; ighost < num_local + num_ghosts; ighost++) { + const int element_face_vec_size = element_face_vector[ighost].size (); + EXPECT_GE (element_face_vec_size, 1) << "Ghost element " << ighost << " must have at least one face index."; + for (int iface = 0; iface < element_face_vec_size; ++iface) { + EXPECT_GE (element_face_vector[ighost][iface], 0); + EXPECT_LE (element_face_vector[ighost][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) << "Face " << iface << " must have at least one side."; + const auto elem_first = (*mesh)[face.sides[0].element_id]; + if (elem_first.is_ghost_element ()) { + continue; + } + std::vector dual_faces; + auto neighs = elem_first.get_face_neighbors (face.sides[0].local_face_id, dual_faces); + // Check for the case that this face is a 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 must be local."; + EXPECT_EQ (neighs.size (), 0) << "BOUNDARY side must not have neighbors."; + } + 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.sides[0].orientation, face.sides[1].orientation) + << "ATTENTION if this is correct we can pull orientation information to the face instead of face sides."; + } + else if (face.type == face_type::MORTAR || face.type == face_type::MPI_MORTAR) { + EXPECT_GT (face.sides.size (), 2) << "MORTAR face must have more than 2 sides."; + // Check that the first side is the large side (num_neighs>1). + EXPECT_EQ (neighs.size (), face.sides.size () - 1) + << "MORTAR side must have as many neighbors as small sides of the face."; + 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); + 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."; + if (!neighs[ineigh]->is_ghost_element ()) { + 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.sides[0].orientation, (*neigh_face_side).orientation) + << "ATTENTION if this is correct we can pull orientation information to the face instead of face sides."; + } + } + + // This ensures that the element_face_vector is filled correctly and together with the check that + for (const auto& side : face.sides) { + bool has_face_in_element_face_vec + = std::any_of (element_face_vector[side.element_id].begin (), element_face_vector[side.element_id].end (), + [iface] (const int faceidx) { return iface == faceidx; }); + EXPECT_TRUE (has_face_in_element_face_vec) + << "Face index in element_face_vector does not match face index in faces vector."; + } + } +} From 1e1e5d5eebb9eec7428e943f97833eb266650b0a Mon Sep 17 00:00:00 2001 From: Lena Ploetzke Date: Wed, 13 May 2026 16:33:50 +0200 Subject: [PATCH 10/12] current status --- mesh_handle/competences/dg_competences.hxx | 64 ++++-- .../competences/element_data_competences.hxx | 3 +- .../competences/t8_gtest_dg_competences.cxx | 190 +++++++++--------- 3 files changed, 139 insertions(+), 118 deletions(-) diff --git a/mesh_handle/competences/dg_competences.hxx b/mesh_handle/competences/dg_competences.hxx index c4ab11086e..94c7a3ea9c 100644 --- a/mesh_handle/competences/dg_competences.hxx +++ b/mesh_handle/competences/dg_competences.hxx @@ -34,7 +34,7 @@ along with t8code; if not, write to the Free Software Foundation, Inc., #include #include #include -#include +#include namespace t8_mesh_handle { @@ -185,7 +185,7 @@ struct face_vector_mesh_competence: public t8_crtp_operatorunderlying (); // Ensure that rank vector is set. - if constexpr (mesh.has_remote_ranks_mesh_competence ()) { + if constexpr (TUnderlying::has_remote_ranks_mesh_competence ()) { mesh.set_rank_vector (); } else { @@ -200,7 +200,15 @@ struct face_vector_mesh_competence: public t8_crtp_operator dual_faces; auto neighs = elem.get_face_neighbors (iface, dual_faces); @@ -214,7 +222,7 @@ struct face_vector_mesh_competence: public t8_crtp_operator (m_faces.size ()); m_faces.push_back (std::move (f)); - m_element_face_vector[handle_id].push_back (face_idx); + m_element_face_vector[handle_id][iface] = face_idx; } else if (num_neighs == 1) { const int neigh_id = neighs[0]->get_element_handle_id (); @@ -237,8 +245,11 @@ struct face_vector_mesh_competence: public t8_crtp_operator (m_faces.size ()); m_faces.push_back (std::move (f)); - m_element_face_vector[handle_id].push_back (face_idx); - m_element_face_vector[neigh_id].push_back (face_idx); + 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 ()) { @@ -249,20 +260,32 @@ struct face_vector_mesh_competence: public t8_crtp_operatorget_local_tree_id (), t8_forest_get_scheme (forest), neighs[0]->get_forest_element (), dual_faces[0]); - f.sides.push_back ({ neigh_id, dual_faces[0], orientation_neigh, neigh_rank }); const int 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, orientation, LOCAL_RANK }); + face f { face_type::MPI_MORTAR, + { { neigh_id, dual_faces[0], orientation_neigh, neigh_rank }, + { handle_id, iface, orientation, LOCAL_RANK } } }; const int face_idx = static_cast (m_faces.size ()); m_faces.push_back (std::move (f)); - m_element_face_vector[handle_id].push_back (face_idx); - m_element_face_vector[neigh_id].push_back (face_idx); + 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 { SC_ABORT ("Not possible without incomplete trees.\n"); @@ -271,16 +294,12 @@ struct face_vector_mesh_competence: public t8_crtp_operator 1) --- */ // We are the coarse element; neighs are all finer small sides. - bool any_remote = false; - for (int ineigh = 0; ineigh < num_neighs; ++ineigh) { - if (mesh.get_rank (neighs[ineigh]->get_element_handle_id ()) != LOCAL_RANK) { - any_remote = true; - break; - } - } + 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].push_back (face_idx); + m_element_face_vector[handle_id][iface] = face_idx; face f; f.type = any_remote ? face_type::MPI_MORTAR : face_type::MORTAR; @@ -295,7 +314,10 @@ struct face_vector_mesh_competence: public t8_crtp_operatorget_forest_element (), dual_faces[ineigh]); f.sides.push_back ({ neigh_id, dual_faces[ineigh], orientation_neigh, neigh_rank }); // Record face index for the small side element. - m_element_face_vector[neighs[ineigh]->get_element_handle_id ()].push_back (face_idx); + 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)); diff --git a/mesh_handle/competences/element_data_competences.hxx b/mesh_handle/competences/element_data_competences.hxx index 445146ceb6..3d52fda837 100644 --- a/mesh_handle/competences/element_data_competences.hxx +++ b/mesh_handle/competences/element_data_competences.hxx @@ -63,8 +63,7 @@ class element_data_mesh_competence_impl: public t8_crtp_basic { void set_element_data (std::vector element_data) { - const auto num_local_elements = this->underlying ().get_num_local_elements (); - T8_ASSERT (element_data.size () == static_cast (num_local_elements)); + T8_ASSERT (element_data.size () == static_cast (this->underlying ().get_num_local_elements ())); // element_data is moved, not copied. m_element_data = std::move (element_data); } diff --git a/test/mesh_handle/competences/t8_gtest_dg_competences.cxx b/test/mesh_handle/competences/t8_gtest_dg_competences.cxx index b02c367376..b5a9193bbf 100644 --- a/test/mesh_handle/competences/t8_gtest_dg_competences.cxx +++ b/test/mesh_handle/competences/t8_gtest_dg_competences.cxx @@ -116,101 +116,101 @@ TEST (t8_gtest_dg_competences, face_vector_mesh_competence) // 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 (const auto& elem : *mesh) { - const int handle_id = elem.get_element_handle_id (); - EXPECT_EQ (static_cast (element_face_vector[handle_id].size ()), elem.get_num_faces ()) - << "Element " << handle_id << " has wrong number of face indices."; - // Check that the element_face_vector points to valid face indices in the faces vector. - for (int iface = 0; iface < elem.get_num_faces (); ++iface) { - EXPECT_GE (element_face_vector[handle_id][iface], 0); - EXPECT_LE (element_face_vector[handle_id][iface], static_cast (faces.size ())); - } - } - for (t8_locidx_t ighost = num_local; ighost < num_local + num_ghosts; ighost++) { - const int element_face_vec_size = element_face_vector[ighost].size (); - EXPECT_GE (element_face_vec_size, 1) << "Ghost element " << ighost << " must have at least one face index."; - for (int iface = 0; iface < element_face_vec_size; ++iface) { - EXPECT_GE (element_face_vector[ighost][iface], 0); - EXPECT_LE (element_face_vector[ighost][iface], static_cast (faces.size ())); - } - } + // 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 (const auto& elem : *mesh) { + // const int handle_id = elem.get_element_handle_id (); + // EXPECT_EQ (static_cast (element_face_vector[handle_id].size ()), elem.get_num_faces ()) + // << "Element " << handle_id << " has wrong number of face indices."; + // // Check that the element_face_vector points to valid face indices in the faces vector. + // for (int iface = 0; iface < elem.get_num_faces (); ++iface) { + // EXPECT_GE (element_face_vector[handle_id][iface], 0); + // EXPECT_LE (element_face_vector[handle_id][iface], static_cast (faces.size ())); + // } + // } + // for (t8_locidx_t ighost = num_local; ighost < num_local + num_ghosts; ighost++) { + // const int element_face_vec_size = element_face_vector[ighost].size (); + // EXPECT_GE (element_face_vec_size, 1) << "Ghost element " << ighost << " must have at least one face index."; + // for (int iface = 0; iface < element_face_vec_size; ++iface) { + // EXPECT_GE (element_face_vector[ighost][iface], 0); + // EXPECT_LE (element_face_vector[ighost][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) << "Face " << iface << " must have at least one side."; - const auto elem_first = (*mesh)[face.sides[0].element_id]; - if (elem_first.is_ghost_element ()) { - continue; - } - std::vector dual_faces; - auto neighs = elem_first.get_face_neighbors (face.sides[0].local_face_id, dual_faces); - // Check for the case that this face is a 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 must be local."; - EXPECT_EQ (neighs.size (), 0) << "BOUNDARY side must not have neighbors."; - } - 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.sides[0].orientation, face.sides[1].orientation) - << "ATTENTION if this is correct we can pull orientation information to the face instead of face sides."; - } - else if (face.type == face_type::MORTAR || face.type == face_type::MPI_MORTAR) { - EXPECT_GT (face.sides.size (), 2) << "MORTAR face must have more than 2 sides."; - // Check that the first side is the large side (num_neighs>1). - EXPECT_EQ (neighs.size (), face.sides.size () - 1) - << "MORTAR side must have as many neighbors as small sides of the face."; - 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); - 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."; - if (!neighs[ineigh]->is_ghost_element ()) { - 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.sides[0].orientation, (*neigh_face_side).orientation) - << "ATTENTION if this is correct we can pull orientation information to the face instead of face sides."; - } - } - - // This ensures that the element_face_vector is filled correctly and together with the check that - for (const auto& side : face.sides) { - bool has_face_in_element_face_vec - = std::any_of (element_face_vector[side.element_id].begin (), element_face_vector[side.element_id].end (), - [iface] (const int faceidx) { return iface == faceidx; }); - EXPECT_TRUE (has_face_in_element_face_vec) - << "Face index in element_face_vector does not match face index in 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) << "Face " << iface << " must have at least one side."; + // const auto elem_first = (*mesh)[face.sides[0].element_id]; + // if (elem_first.is_ghost_element ()) { + // continue; + // } + // std::vector dual_faces; + // auto neighs = elem_first.get_face_neighbors (face.sides[0].local_face_id, dual_faces); + // // Check for the case that this face is a 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 must be local."; + // EXPECT_EQ (neighs.size (), 0) << "BOUNDARY side must not have neighbors."; + // } + // 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.sides[0].orientation, face.sides[1].orientation) + // << "ATTENTION if this is correct we can pull orientation information to the face instead of face sides."; + // } + // else if (face.type == face_type::MORTAR || face.type == face_type::MPI_MORTAR) { + // //EXPECT_GT (face.sides.size (), 2) << "MORTAR face must have more than 2 sides."; + // // Check that the first side is the large side (num_neighs>1). + // EXPECT_EQ (neighs.size (), face.sides.size () - 1) + // << "MORTAR side must have as many neighbors as small sides of the face."; + // 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); + // // 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."; + // // if (!neighs[ineigh]->is_ghost_element ()) { + // // 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.sides[0].orientation, (*neigh_face_side).orientation) + // // << "ATTENTION if this is correct we can pull orientation information to the face instead of face sides."; + // // } + // } + + // // This ensures that the element_face_vector is filled correctly and together with the check that + // for (const auto& side : face.sides) { + // bool has_face_in_element_face_vec + // = std::any_of (element_face_vector[side.element_id].begin (), element_face_vector[side.element_id].end (), + // [iface] (const int faceidx) { return iface == faceidx; }); + // EXPECT_TRUE (has_face_in_element_face_vec) + // << "Face index in element_face_vector does not match face index in faces vector."; + // } + // } } From 604349cc75e4d380550e27980ebd8939aed02cfa Mon Sep 17 00:00:00 2001 From: Lena Ploetzke Date: Wed, 20 May 2026 11:23:19 +0200 Subject: [PATCH 11/12] Working test --- mesh_handle/competences/dg_competences.hxx | 31 +-- mesh_handle/mesh.hxx | 2 +- .../competences/t8_gtest_dg_competences.cxx | 187 +++++++++--------- 3 files changed, 106 insertions(+), 114 deletions(-) diff --git a/mesh_handle/competences/dg_competences.hxx b/mesh_handle/competences/dg_competences.hxx index 94c7a3ea9c..b5c1c4c04f 100644 --- a/mesh_handle/competences/dg_competences.hxx +++ b/mesh_handle/competences/dg_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 dg_element_competences.hxx +/** \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. * TODO: IDEA for test: write rank on element data and then check if this is correct for the ghosts @@ -156,7 +156,7 @@ struct face * 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: + * \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), @@ -184,6 +184,7 @@ struct face_vector_mesh_competence: public t8_crtp_operatorunderlying (); + T8_ASSERT (mesh.is_balanced ()); // Ensure that rank vector is set. if constexpr (TUnderlying::has_remote_ranks_mesh_competence ()) { mesh.set_rank_vector (); @@ -191,6 +192,7 @@ struct face_vector_mesh_competence: public t8_crtp_operatorincomplete_trees == 0, "This functionality is not specified for incomplete trees.\n"); @@ -205,10 +207,10 @@ struct face_vector_mesh_competence: public t8_crtp_operator dual_faces; auto neighs = elem.get_face_neighbors (iface, dual_faces); @@ -216,10 +218,7 @@ struct face_vector_mesh_competence: public t8_crtp_operator (m_faces.size ()); m_faces.push_back (std::move (f)); m_element_face_vector[handle_id][iface] = face_idx; @@ -234,10 +233,11 @@ struct face_vector_mesh_competence: public t8_crtp_operatorget_local_tree_id (), t8_forest_get_scheme (forest), neighs[0]->get_forest_element (), dual_faces[0]); @@ -258,9 +258,10 @@ struct face_vector_mesh_competence: public t8_crtp_operator (m_faces.size ()); @@ -293,7 +294,7 @@ struct face_vector_mesh_competence: public t8_crtp_operator 1) --- */ - // We are the coarse element; neighs are all finer small sides. + // 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; }); @@ -313,7 +314,7 @@ struct face_vector_mesh_competence: public t8_crtp_operatorget_local_tree_id (), t8_forest_get_scheme (forest), neighs[ineigh]->get_forest_element (), dual_faces[ineigh]); f.sides.push_back ({ neigh_id, dual_faces[ineigh], orientation_neigh, neigh_rank }); - // Record face index for the small side element. + // Record face index for the small mortar. if (m_element_face_vector[neigh_id].empty ()) { m_element_face_vector[neigh_id].assign (neighs[ineigh]->get_num_faces (), -1); } diff --git a/mesh_handle/mesh.hxx b/mesh_handle/mesh.hxx index 37edd8642d..4c88c0a4bf 100644 --- a/mesh_handle/mesh.hxx +++ b/mesh_handle/mesh.hxx @@ -177,7 +177,7 @@ class mesh: public TMeshCompetencePack::template applyset_ghost (); mesh->commit (); - if ((mesh->get_dimension () > 1) && (mesh->get_num_local_elements () > 1)) { - ASSERT_GT (mesh->get_num_ghosts (), 0); - } // 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 (const auto& elem : *mesh) { - // const int handle_id = elem.get_element_handle_id (); - // EXPECT_EQ (static_cast (element_face_vector[handle_id].size ()), elem.get_num_faces ()) - // << "Element " << handle_id << " has wrong number of face indices."; - // // Check that the element_face_vector points to valid face indices in the faces vector. - // for (int iface = 0; iface < elem.get_num_faces (); ++iface) { - // EXPECT_GE (element_face_vector[handle_id][iface], 0); - // EXPECT_LE (element_face_vector[handle_id][iface], static_cast (faces.size ())); - // } - // } - // for (t8_locidx_t ighost = num_local; ighost < num_local + num_ghosts; ighost++) { - // const int element_face_vec_size = element_face_vector[ighost].size (); - // EXPECT_GE (element_face_vec_size, 1) << "Ghost element " << ighost << " must have at least one face index."; - // for (int iface = 0; iface < element_face_vec_size; ++iface) { - // EXPECT_GE (element_face_vector[ighost][iface], 0); - // EXPECT_LE (element_face_vector[ighost][iface], static_cast (faces.size ())); - // } - // } + 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) << "Face " << iface << " must have at least one side."; - // const auto elem_first = (*mesh)[face.sides[0].element_id]; - // if (elem_first.is_ghost_element ()) { - // continue; - // } - // std::vector dual_faces; - // auto neighs = elem_first.get_face_neighbors (face.sides[0].local_face_id, dual_faces); - // // Check for the case that this face is a 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 must be local."; - // EXPECT_EQ (neighs.size (), 0) << "BOUNDARY side must not have neighbors."; - // } - // 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.sides[0].orientation, face.sides[1].orientation) - // << "ATTENTION if this is correct we can pull orientation information to the face instead of face sides."; - // } - // else if (face.type == face_type::MORTAR || face.type == face_type::MPI_MORTAR) { - // //EXPECT_GT (face.sides.size (), 2) << "MORTAR face must have more than 2 sides."; - // // Check that the first side is the large side (num_neighs>1). - // EXPECT_EQ (neighs.size (), face.sides.size () - 1) - // << "MORTAR side must have as many neighbors as small sides of the face."; - // 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); - // // 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."; - // // if (!neighs[ineigh]->is_ghost_element ()) { - // // 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.sides[0].orientation, (*neigh_face_side).orientation) - // // << "ATTENTION if this is correct we can pull orientation information to the face instead of face sides."; - // // } - // } - - // // This ensures that the element_face_vector is filled correctly and together with the check that - // for (const auto& side : face.sides) { - // bool has_face_in_element_face_vec - // = std::any_of (element_face_vector[side.element_id].begin (), element_face_vector[side.element_id].end (), - // [iface] (const int faceidx) { return iface == faceidx; }); - // EXPECT_TRUE (has_face_in_element_face_vec) - // << "Face index in element_face_vector does not match face index in 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."; + } + // --- 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.sides[0].orientation, face.sides[1].orientation) + << "ATTENTION if this is correct we can pull orientation information to the face instead of face sides."; + } + 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 and has more than one neighbors. + // 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_GT (neighs[0]->get_level (), elem_first.get_level ()) << "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.sides[0].orientation, (*neigh_face_side).orientation) + << "ATTENTION if this is correct we can pull orientation information to the face instead of face sides."; + } + } + + // 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); + } + } } From 2e1a160f03bee82ddeb6f7450fa5d9a58bc87d34 Mon Sep 17 00:00:00 2001 From: Lena Ploetzke Date: Wed, 20 May 2026 13:44:10 +0200 Subject: [PATCH 12/12] documentation and orientation --- mesh_handle/competences/README.md | 3 +- mesh_handle/competences/dg_competences.hxx | 62 ++++++++----------- test/mesh_handle/README.md | 3 +- .../competences/t8_gtest_dg_competences.cxx | 57 +++++++++++------ .../t8_gtest_adapt_partition_balance.cxx | 2 +- 5 files changed, 68 insertions(+), 59 deletions(-) diff --git a/mesh_handle/competences/README.md b/mesh_handle/competences/README.md index 39c668a5ad..c556568907 100644 --- a/mesh_handle/competences/README.md +++ b/mesh_handle/competences/README.md @@ -4,6 +4,7 @@ Definition of additional competences/functionalities that can be used for the t8 - [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: @@ -11,7 +12,7 @@ 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, this is required. +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. diff --git a/mesh_handle/competences/dg_competences.hxx b/mesh_handle/competences/dg_competences.hxx index b5c1c4c04f..2cf20858f9 100644 --- a/mesh_handle/competences/dg_competences.hxx +++ b/mesh_handle/competences/dg_competences.hxx @@ -23,7 +23,6 @@ along with t8code; if not, write to the Free Software Foundation, Inc., /** \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. - * TODO: IDEA for test: write rank on element data and then check if this is correct for the ghosts */ #pragma once @@ -64,16 +63,17 @@ struct remote_ranks_mesh_competence: public t8_crtp_operatorunderlying (); 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 --- */ + /* --- 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. --- */ + /* --- 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; @@ -101,6 +101,7 @@ struct remote_ranks_mesh_competence: public t8_crtp_operator sides; + int orientation = 0; ///< Face orientation code for coordinate permutation. }; /** @@ -173,8 +174,8 @@ struct face_vector_mesh_competence: public t8_crtp_operator (m_faces.size ()); m_faces.push_back (std::move (f)); m_element_face_vector[handle_id][iface] = face_idx; @@ -233,15 +234,10 @@ struct face_vector_mesh_competence: public t8_crtp_operatorget_local_tree_id (), t8_forest_get_scheme (forest), - neighs[0]->get_forest_element (), dual_faces[0]); - f.sides.push_back ({ neigh_id, dual_faces[0], orientation_neigh, neigh_rank }); const int face_idx = static_cast (m_faces.size ()); m_faces.push_back (std::move (f)); @@ -263,22 +259,18 @@ struct face_vector_mesh_competence: public t8_crtp_operatorget_local_tree_id (), t8_forest_get_scheme (forest), - neighs[0]->get_forest_element (), dual_faces[0]); + 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], orientation_neigh, neigh_rank }, // Large mortar first. - { handle_id, iface, orientation, LOCAL_RANK } } }; + { { 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)); @@ -288,9 +280,7 @@ struct face_vector_mesh_competence: public t8_crtp_operator 1) --- */ @@ -304,17 +294,15 @@ struct face_vector_mesh_competence: public t8_crtp_operatorget_element_handle_id (); const int neigh_rank = mesh.get_rank (neigh_id); - const int orientation_neigh = t8_forest_leaf_face_orientation ( - forest, neighs[ineigh]->get_local_tree_id (), t8_forest_get_scheme (forest), - neighs[ineigh]->get_forest_element (), dual_faces[ineigh]); - f.sides.push_back ({ neigh_id, dual_faces[ineigh], orientation_neigh, neigh_rank }); - // Record face index for the small mortar. + 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); } @@ -338,9 +326,9 @@ struct face_vector_mesh_competence: public t8_crtp_operator>& get_element_face_vector () const { diff --git a/test/mesh_handle/README.md b/test/mesh_handle/README.md index e48ca10a2a..b48729350a 100644 --- a/test/mesh_handle/README.md +++ b/test/mesh_handle/README.md @@ -13,4 +13,5 @@ Therefore, the tests are structured as follows: The folder [competences/](competences/) checks element and mesh competences of the mesh handle. We have the following tests: - [t8_gtest_custom_competence.cxx](t8_gtest_custom_competence.cxx) checks that user defined competences can be used for the mesh handle. - [t8_gtest_cache_competence.cxx](t8_gtest_cache_competence.cxx) tests that *all predefined caching competences* work as intended. -- [t8_gtest_handle_data.cxx](t8_gtest_handle_data.cxx) tests that user and element data functionality works as intended. The element data functionality is provided via a competence. \ No newline at end of file +- [t8_gtest_handle_data.cxx](t8_gtest_handle_data.cxx) tests that user and element data functionality works as intended. The element data functionality is provided via a competence. +- [t8_gtest_dg_competences.cxx](t8_gtest_dg_competences.cxx) checks the competences related to the discontinuous Galerkin methods. \ No newline at end of file diff --git a/test/mesh_handle/competences/t8_gtest_dg_competences.cxx b/test/mesh_handle/competences/t8_gtest_dg_competences.cxx index 039638ae23..55d7c38f69 100644 --- a/test/mesh_handle/competences/t8_gtest_dg_competences.cxx +++ b/test/mesh_handle/competences/t8_gtest_dg_competences.cxx @@ -22,10 +22,9 @@ along with t8code; if not, write to the Free Software Foundation, Inc., /** * \file t8_gtest_dg_competences.cxx - * Checks that the competences for discontinuous Galerkin methods work as expected. + * Checks that the competences for discontinuous Galerkin methods defined in \ref dg_competences.hxx work as expected. */ #include -#include #include #include @@ -38,9 +37,10 @@ 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_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. @@ -59,7 +59,8 @@ mesh_adapt_callback_test_refine_second ([[maybe_unused]] const TMeshClass& mesh, return 0; } -/** TODO +/** 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) { @@ -74,9 +75,11 @@ TEST (t8_gtest_dg_competences, remote_ranks) mesh->set_ghost (); mesh->commit (); - if ((mesh->get_dimension () > 1) && (mesh->get_num_local_elements () > 1)) { + 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 (mesh->get_num_ghosts (), 0); + ASSERT_GT (num_ghosts, 0); } int mpirank; @@ -84,24 +87,27 @@ TEST (t8_gtest_dg_competences, remote_ranks) SC_CHECK_MPI (mpiret); // Set local rank for all local mesh elements. - std::vector element_data (mesh->get_num_local_elements (), { 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 functionality works as expected. + // 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 = mesh->get_num_local_elements (); - ighost < mesh->get_num_local_elements () + mesh->get_num_ghosts (); ighost++) { + 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 = 1; + 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); @@ -144,14 +150,19 @@ TEST (t8_gtest_dg_competences, face_vector_mesh_competence) std::vector dual_faces; auto neighs = elem_first.get_face_neighbors (face.sides[0].local_face_id, dual_faces); - // --- BOUNDARY --- + /* --- 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 --- + /* --- 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."; @@ -166,20 +177,25 @@ TEST (t8_gtest_dg_competences, face_vector_mesh_competence) 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.sides[0].orientation, face.sides[1].orientation) - << "ATTENTION if this is correct we can pull orientation information to the face instead of face sides."; + 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 and has more than one neighbors. + // 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_GT (neighs[0]->get_level (), elem_first.get_level ()) << "MORTAR first element must be the large mortar."; + 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."; @@ -194,8 +210,11 @@ TEST (t8_gtest_dg_competences, face_vector_mesh_competence) << "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.sides[0].orientation, (*neigh_face_side).orientation) - << "ATTENTION if this is correct we can pull orientation information to the face instead of face sides."; + 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."; } } 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.