diff --git a/.typos.toml b/.typos.toml index b956ef5a00..eb06f3ed0d 100644 --- a/.typos.toml +++ b/.typos.toml @@ -1,6 +1,7 @@ [default.extend-words] eles = "eles" packageid = "packageid" +countr = "countr" [files] extend-exclude = ["scripts/indent.sh", "thirdparty/", "t8code_logo.png", "cmake/FindOpenCASCADE.cmake", "src/t8_misc/t8_with_macro_error.h", "doc/Doxyfile.in"] diff --git a/example/CMakeLists.txt b/example/CMakeLists.txt index d2d4d918be..e7b7355d26 100644 --- a/example/CMakeLists.txt +++ b/example/CMakeLists.txt @@ -95,6 +95,9 @@ add_t8_example( NAME t8_example_spheres SOURCES remove/t8_exampl add_t8_example( NAME t8_example_gauss_blob SOURCES remove/t8_example_gauss_blob.cxx ) add_t8_example( NAME t8_example_empty_trees SOURCES remove/t8_example_empty_trees.cxx ) +add_t8_example( NAME t8_example_hanging_nodes_shock SOURCES subelements/t8_hanging_nodes_shock.cxx ) +add_t8_example( NAME t8_example_hanging_nodes SOURCES subelements/t8_quads_hanging_nodes.cxx ) + add_t8_example( NAME t8_version SOURCES version/t8_version.cxx ) # NOTE: The following examples are (currently) deprecated and no longer compiled. diff --git a/example/subelements/t8_hanging_nodes_shock.cxx b/example/subelements/t8_hanging_nodes_shock.cxx new file mode 100644 index 0000000000..8574994f30 --- /dev/null +++ b/example/subelements/t8_hanging_nodes_shock.cxx @@ -0,0 +1,252 @@ +/* + 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_hanging_nodes_shock.cxx + * Example demonstrating hanging-node resolution on a hybrid 2D mesh. + * + * The program builds a uniform forest on a hybrid (quad + triangle) 2D hypercube + * using the subelement scheme, then repeatedly grades the mesh around a circle and + * resolves the resulting hanging nodes. Each stage is written to VTK so the whole + * process can be inspected: + * 1. uniform forest, + * 2. adapt (refine near a circle) -> hanging nodes appear, + * 3. balance (enforce a 2:1 level difference between neighbors), + * 4. remove hanging nodes -> transition cells are split into subelements, + * 5. discard subelements -> back to a plain (recursively refined) forest, + * 6. a second adapt / balance / remove cycle to show the process is repeatable. + * + * Subelements are the mechanism that keeps the mesh conformal: where balancing + * leaves a coarse element adjacent to finer ones (a hanging node), that element is + * transitioned into a fan of smaller subelements so no hanging nodes remain. + */ + +#include /* General t8code header, always include this. */ +#include /* cmesh definition and basic interface. */ +#include /* A collection of exemplary cmeshes */ +#include /* forest definition and basic interface. */ +#include /* save forest */ +#include /* geometrical information of the forest */ +#include /* Function for adding subelements. */ +#include /* Subelement refinement scheme. */ +#include /* Basic operations on 3D vectors. */ +#include /* Element class (eclass) definitions. */ + +/** User data passed to the adaptation callback \ref t8_adapt_callback. + * Defines the circle the mesh is refined around and the level bounds. */ +struct t8_adapt_data +{ + double midpoint[3]; /* Center of the circle the mesh is refined around. */ + double radius; /* Radius of that circle; the mesh is refined near its boundary. */ + double delta; /* Width of the transition band around the circle over which the level ranges. */ + int minlevel; /* Coarsest level, reached at distance >= delta from the circle. */ + int maxlevel; /* Finest level, reached on the circle itself. */ +}; + +/** The adaptation callback function. + * Adapts the mesh around a circle of radius \a radius centered at \a midpoint: + * Elements on the circle are refined to \a maxlevel, relaxing linearly to \a minlevel over a band of width \a delta. + * The closer an element is to the circle, the finer it gets: elements right on the circle are refined to maxlevel, + * elements delta or more away stay at minlevel, and in between the level scales linearly with the distance to the + * circle. + * \param [in] forest The current forest that is in construction. + * \param [in] forest_from The forest from which we adapt (here, the uniform forest). + * \param [in] which_tree The process local id of the current tree. + * \param [in] tree_class The eclass of \a which_tree. + * \param [in] lelement_id The tree local index of the current element (or first of the family). + * \param [in] scheme The refinement scheme for this tree's element class. + * \param [in] is_family If 1, the first \a num_elements entries in \a elements form a family. + * \param [in] num_elements The number of entries in \a elements that are defined. + * \param [in] elements The element or family to consider for refinement/coarsening. + * \return 1 to refine, -1 to coarsen the family, 0 to leave unchanged. + */ +int +t8_adapt_callback (t8_forest_t forest, t8_forest_t forest_from, t8_locidx_t which_tree, t8_eclass_t tree_class, + [[maybe_unused]] t8_locidx_t lelement_id, const t8_scheme *scheme, const int is_family, + [[maybe_unused]] const int num_elements, t8_element_t *elements[]) +{ + const auto *adapt_data = (const struct t8_adapt_data *) t8_forest_get_user_data (forest); + + T8_ASSERT (adapt_data != NULL); + + /* Compute the element's centroid coordinates. */ + double centroid[3]; + t8_forest_element_centroid (forest_from, which_tree, elements[0], centroid); + + /* Distance to the circle center, then to the circle boundary. */ + double radius = t8_dist (centroid, adapt_data->midpoint); + double abs_to_radius = fabs (radius - adapt_data->radius); + const int level = scheme->element_get_level (tree_class, elements[0]); + + /* Normalized shell distance in [0, 1]; target level is between maxlevel and minlevel. */ + double alpha = std::min (abs_to_radius / adapt_data->delta, 1.0); + int target_level + = adapt_data->maxlevel - static_cast (std::round (alpha * (adapt_data->maxlevel - adapt_data->minlevel))); + + if ((level < target_level) && (level < adapt_data->maxlevel)) { + /* Coarser than target: refine. */ + return 1; + } + else if ((is_family && level > target_level) && (level > adapt_data->minlevel)) { + /* Finer than target: coarsen the family. Check is_family first. */ + return -1; + } + /* At target level: leave unchanged. */ + return 0; +} + +/** Adapt a forest around the circle of radius 0.45 (first adaptation cycle). + * \param[in] forest Forest to be adapted. */ +t8_forest_t +t8_adapt_forest (t8_forest_t forest) +{ + struct t8_adapt_data adapt_data = { + { 0, 1, 0 }, /* Center of the circle. */ + 0.45, /* Radius */ + 0.1, /* Delta (transition band width) */ + 2, /* Minlevel */ + 6 /* Maxlevel */ + }; + + t8_forest_t forest_adapt; + forest_adapt = t8_forest_new_adapt (forest, t8_adapt_callback, 1, 0, &adapt_data); + return forest_adapt; +} + +/** Adapt a forest around the circle of radius 0.6 (second adaptation cycle). + * Same criterion as \ref t8_adapt_forest but with a larger radius. + * \param[in] forest Forest to be adapted. + */ +t8_forest_t +t8_adapt_forest_2and (t8_forest_t forest) +{ + struct t8_adapt_data adapt_data = { + { 0, 1, 0 }, /* Center of the circle. */ + 0.6, /* Radius */ + 0.1, /* Delta (transition band width) */ + 2, /* Minlevel */ + 6 /* Maxlevel */ + }; + + t8_forest_t forest_adapt; + forest_adapt = t8_forest_new_adapt (forest, t8_adapt_callback, 1, 0, &adapt_data); + return forest_adapt; +} + +/** Balance a forest, i.e. enforce that neighboring elements differ by at most one + * refinement level (a 2:1 balance). + * \param[in] forest Forest to be balanced. + */ +t8_forest_t +t8_forest_balance (t8_forest_t forest) +{ + t8_forest_t forest_new; + t8_forest_init (&forest_new); + t8_forest_set_balance (forest_new, forest, 0); + t8_forest_commit (forest_new); + return forest_new; +} + +/** Entry point of the program. + * + * Runs the full demonstration pipeline (uniform -> adapt -> balance -> + * remove hanging nodes -> discard -> adapt -> balance -> remove), writing the + * forest to VTK after each stage. + */ +int +main (int argc, char **argv) +{ + /* Initialize MPI, libsc, and t8code. */ + int mpiret = sc_MPI_Init (&argc, &argv); + SC_CHECK_MPI (mpiret); + sc_init (sc_MPI_COMM_WORLD, 1, 1, NULL, SC_LP_ESSENTIAL); + t8_init (SC_LP_PRODUCTION); + /* We will use MPI_COMM_WORLD as a communicator. */ + sc_MPI_Comm comm = sc_MPI_COMM_WORLD; + + /* --- Setup: build the cmesh and a uniform forest. --- */ + /* Hybrid 2D hypercube: a mesh containing both quad and triangle trees. */ + t8_cmesh_t cmesh; + t8_cmesh_init (&cmesh); + t8_cmesh_new_2D_hypercube_hybrid (cmesh, comm); + /* Uniform forest using the subelement scheme (required for hanging-node resolution). */ + const int level = 0; + t8_forest_t forest = t8_forest_new_uniform (cmesh, t8_scheme_new_subelement (), level, 0, comm); + const char *prefix = "t8_uniform"; + t8_forest_write_vtk (forest, prefix); + t8_global_productionf (" [subelements] Uniform forest wrote to file: %s*\n", prefix); + + /* --- Adapt the forest: refine near the first circle, creating hanging nodes. --- */ + forest = t8_adapt_forest (forest); + std::cout << "Subelements before removing: " << t8_forest_has_global_subelements (forest) << std::endl; + prefix = "t8_adapted1"; + t8_forest_write_vtk (forest, prefix); + t8_global_productionf (" [subelements] Wrote adapted forest with hanging nodes to vtu files: %s*\n", prefix); + + /* --- Balance the forest (2:1 balance between neighboring elements). --- */ + forest = t8_forest_balance (forest); + prefix = "t8_balanced1"; + t8_forest_write_vtk (forest, prefix); + t8_global_productionf (" [subelements] Balanced and wrote to file: %s*\n", prefix); + + /* --- Resolve hanging nodes by transitioning elements into subelements. --- */ + forest = t8_forest_remove_hanging_nodes (forest); + std::cout << "Subelements after removing: " << t8_forest_has_global_subelements (forest) << std::endl; + const char *prefix_without_hanging_nodes = "t8_resolved_hanging_nodes1"; + t8_forest_write_vtk (forest, prefix_without_hanging_nodes); + t8_global_productionf (" [subelements] Wrote adapted forest with resolved hanging nodes to vtu files: %s*\n", + prefix_without_hanging_nodes); + + /* --- Discard the subelements to recover a plain, recursively refined forest. --- */ + /* This is the inverse of the previous step and is required before adapting again. */ + forest = t8_forest_discard_subelements (forest); + std::cout << "Subelements removed: " << t8_forest_has_global_subelements (forest) << std::endl; + const char *prefix_removed_sub = "t8_discarded_subelements1"; + t8_forest_write_vtk (forest, prefix_removed_sub); + t8_global_productionf (" [subelements] Wrote adapted forest with discarded subelements to vtu files: %s*\n", + prefix_removed_sub); + + /* --- Second cycle: adapt around the larger circle. --- */ + forest = t8_adapt_forest_2and (forest); + prefix = "t8_adapted2"; + t8_forest_write_vtk (forest, prefix); + t8_global_productionf (" [subelements] Adapted again and wrote to file: %s*\n", prefix); + + /* --- Balance again. --- */ + forest = t8_forest_balance (forest); + prefix = "t8_balanced2"; + t8_forest_write_vtk (forest, prefix); + t8_global_productionf (" [subelements] Balanced again and wrote to file: %s*\n", prefix); + + /* --- Resolve hanging nodes again. --- */ + forest = t8_forest_remove_hanging_nodes (forest); + prefix = "t8_resolved_hanging_nodes2"; + t8_forest_write_vtk (forest, prefix); + t8_global_productionf (" [subelements] Removed hanging nodes after second adaptation and wrote to : %s*\n", prefix); + + /* --- Cleanup: free the forest and finalize t8code / libsc / MPI. --- */ + t8_forest_unref (&forest); + sc_finalize (); + mpiret = sc_MPI_Finalize (); + SC_CHECK_MPI (mpiret); + + return 0; +} diff --git a/example/subelements/t8_quads_hanging_nodes.cxx b/example/subelements/t8_quads_hanging_nodes.cxx new file mode 100644 index 0000000000..d95c15b929 --- /dev/null +++ b/example/subelements/t8_quads_hanging_nodes.cxx @@ -0,0 +1,124 @@ +/* + 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_quads_hanging_nodes.cxx + * This is an example to demonstrate hanging node resolution for quads. + */ + +#include /* General t8code header, always include this. */ +#include /* cmesh definition and basic interface. */ +#include /* A collection of exemplary cmeshes */ +#include /* forest definition and basic interface. */ +#include /* save forest */ +#include /* geometrical information of the forest */ +#include /* Function for adding subelements. */ +#include /* Subelement refinement scheme. */ +#include /* Basic operations on 3D vectors. */ +#include + +/** The adaptation callback function. This refines every second element (with even global id). + * \param [in] forest The current forest that is in construction. + * \param [in] forest_from The forest from which we adapt the current forest (in our case, the uniform forest) + * \param [in] which_tree The process local id of the current tree. + * \param [in] tree_class The eclass of \a which_tree. + * \param [in] lelement_id The tree local index of the current element (or the first of the family). + * \param [in] scheme The refinement scheme for this tree's element class. + * \param [in] is_family If 1, the first \a num_elements entries in \a elements form a family. If 0, they do not. + * \param [in] num_elements The number of entries in \a elements elements that are defined. + * \param [in] elements The element or family of elements to consider for refinement/coarsening. + */ +int +t8_adapt_callback ([[maybe_unused]] t8_forest_t forest, [[maybe_unused]] t8_forest_t forest_from, + t8_locidx_t which_tree, [[maybe_unused]] t8_eclass_t tree_class, t8_locidx_t lelement_id, + [[maybe_unused]] const t8_scheme *scheme, [[maybe_unused]] const int is_family, + [[maybe_unused]] const int num_elements, [[maybe_unused]] t8_element_t *elements[]) +{ + if ((t8_forest_get_tree_element_offset (forest_from, which_tree) + lelement_id) % 2 == 0) { + return 1; + } + return 0; +} + +/** Adapt forest according to callback. */ +t8_forest_t +t8_adapt_forest (t8_forest_t forest) +{ + t8_forest_t forest_adapt; + forest_adapt = t8_forest_new_adapt (forest, t8_adapt_callback, 0, 0, NULL); + return forest_adapt; +} + +/** Entry point of the program. */ +int +main (int argc, char **argv) +{ + /* The uniform refinement level of the forest. */ + const int level = 3; + + int mpiret = sc_MPI_Init (&argc, &argv); + SC_CHECK_MPI (mpiret); + sc_init (sc_MPI_COMM_WORLD, 1, 1, NULL, SC_LP_ESSENTIAL); + t8_init (SC_LP_PRODUCTION); + + /* We will use MPI_COMM_WORLD as a communicator. */ + sc_MPI_Comm comm = sc_MPI_COMM_WORLD; + + /* --- Setup. Build cmesh and uniform forest.--- */ + t8_cmesh_t cmesh; + t8_cmesh_init (&cmesh); + t8_cmesh_new_hypercube (&cmesh, T8_ECLASS_QUAD, comm, 0, 0, 0); + t8_forest_t forest = t8_forest_new_uniform (cmesh, t8_scheme_new_subelement (), level, 0, comm); + + /* --- Adapt the forest. --- */ + forest = t8_adapt_forest (forest); + std::cout << "Subelements before removing: " << t8_forest_has_local_subelements (forest) << std::endl; + const char *prefix_with_hanging_nodes = "t8_with_hanging_nodes"; + t8_forest_write_vtk (forest, prefix_with_hanging_nodes); + t8_global_productionf (" [subelements] Wrote adapted forest with hanging nodes to vtu files: %s*\n", + prefix_with_hanging_nodes); + + /* --- Remove hanging nodes. --- */ + forest = t8_forest_remove_hanging_nodes (forest); + std::cout << "Subelements after removing: " << t8_forest_has_local_subelements (forest) << std::endl; + // Output to vtk. + const char *prefix_without_hanging_nodes = "t8_without_hanging_nodes"; + t8_forest_write_vtk (forest, prefix_without_hanging_nodes); + t8_global_productionf (" [subelements] Wrote adapted forest without hanging nodes to vtu files: %s*\n", + prefix_without_hanging_nodes); + + /* ---Discard subelements. --- */ + forest = t8_forest_discard_subelements (forest); + std::cout << "Subelements removed: " << t8_forest_has_local_subelements (forest) << std::endl; + // Now output to vtk. + const char *prefix_removed_sub = "t8_removed_sub"; + t8_forest_write_vtk (forest, prefix_removed_sub); + t8_global_productionf (" [subelements] Wrote adapted forest with discarded subelements to vtu files: %s*\n", + prefix_removed_sub); + // --- Cleanup. --- + t8_forest_unref (&forest); + + sc_finalize (); + mpiret = sc_MPI_Finalize (); + SC_CHECK_MPI (mpiret); + + return 0; +} diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 6d5197d07e..92252725d4 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -171,6 +171,7 @@ target_sources( T8 PRIVATE t8_forest/t8_forest_ghost.cxx t8_forest/t8_forest_iterate.cxx t8_forest/t8_forest_balance.cxx + t8_forest/t8_forest_subelement.cxx t8_forest/t8_forest_search/t8_forest_search.cxx t8_geometry/t8_geometry.cxx t8_geometry/t8_geometry_helpers.c @@ -207,6 +208,7 @@ target_sources( T8 PRIVATE t8_schemes/t8_default/t8_default_vertex/t8_default_vertex.cxx t8_types/t8_vec.cxx t8_schemes/t8_standalone/t8_standalone.cxx + t8_schemes/t8_subelement/t8_subelement.cxx t8_vtk/t8_vtk.c t8_vtk/t8_vtk_writer.cxx t8_vtk/t8_vtk_write_ASCII.cxx diff --git a/src/t8_cmesh/t8_cmesh_examples.cxx b/src/t8_cmesh/t8_cmesh_examples.cxx index afe23d4d55..ef4d49c92e 100644 --- a/src/t8_cmesh/t8_cmesh_examples.cxx +++ b/src/t8_cmesh/t8_cmesh_examples.cxx @@ -681,6 +681,81 @@ t8_cmesh_new_hypercube_hybrid (t8_cmesh_t cmesh, sc_MPI_Comm comm, int periodic) t8_cmesh_commit (cmesh, comm); } +void +t8_cmesh_new_2D_hypercube_hybrid (t8_cmesh_t cmesh, sc_MPI_Comm comm) +{ + + T8_ASSERT (cmesh != NULL); + T8_ASSERT (t8_cmesh_is_initialized (cmesh)); + T8_ASSERT (!t8_cmesh_is_committed (cmesh, 0)); + T8_ASSERT (t8_cmesh_stash_is_empty (cmesh)); + { + /* clang-format off */ + double vertices[60] = { /* All vertices of all trees. Partly duplicated */ + 0, 0, 0, /* tree 0, triangle */ + 0.5, 0, 0, + 0.5, 0.5, 0, + 0, 0, 0, /* tree 1, triangle */ + 0.5, 0.5, 0, + 0, 0.5, 0, + 0.5, 0, 0, /* tree 2, quad */ + 1, 0, 0, 0.5, + 0.5, 0, 1, 0.5, + 0, 0, 0.5, 0, /* tree 3, quad */ + 0.5, 0.5, 0, + 0, 1, 0, + 0.5, 1, 0, + 0.5, 0.5, 0, /* tree 4, triangle */ + 1, 0.5, 0, + 1, 1, 0, + 0.5, 0.5, 0, /* tree 5, triangle */ + 1, 1, 0, + 0.5, 1, 0 + }; + /* clang-format on */ + + /* + * This is how the cmesh looks like. The numbers are the tree numbers: + * + * +---+---+ + * | |5 /| + * | 3 | / | + * | |/ 4| + * +---+---+ + * |1 /| | + * | / | 2 | + * |/0 | | + * +---+---+ + */ + + /* Use linear geometry */ + t8_cmesh_register_geometry (cmesh); + + t8_cmesh_set_tree_class (cmesh, 0, T8_ECLASS_TRIANGLE); + t8_cmesh_set_tree_class (cmesh, 1, T8_ECLASS_TRIANGLE); + t8_cmesh_set_tree_class (cmesh, 2, T8_ECLASS_QUAD); + t8_cmesh_set_tree_class (cmesh, 3, T8_ECLASS_QUAD); + t8_cmesh_set_tree_class (cmesh, 4, T8_ECLASS_TRIANGLE); + t8_cmesh_set_tree_class (cmesh, 5, T8_ECLASS_TRIANGLE); + + t8_cmesh_set_tree_vertices (cmesh, 0, vertices, 3); + t8_cmesh_set_tree_vertices (cmesh, 1, vertices + 9, 3); + t8_cmesh_set_tree_vertices (cmesh, 2, vertices + 18, 4); + t8_cmesh_set_tree_vertices (cmesh, 3, vertices + 30, 4); + t8_cmesh_set_tree_vertices (cmesh, 4, vertices + 42, 3); + t8_cmesh_set_tree_vertices (cmesh, 5, vertices + 51, 3); + + t8_cmesh_set_join (cmesh, 0, 1, 1, 2, 0); + t8_cmesh_set_join (cmesh, 0, 2, 0, 0, 0); + t8_cmesh_set_join (cmesh, 1, 3, 0, 2, 1); + t8_cmesh_set_join (cmesh, 2, 4, 3, 2, 0); + t8_cmesh_set_join (cmesh, 3, 5, 1, 1, 0); + t8_cmesh_set_join (cmesh, 4, 5, 1, 2, 0); + + t8_cmesh_commit (cmesh, comm); + } +} + /* The unit cube is constructed from trees of the same eclass. * For triangles the square is divided along the (0,0) -- (1,1) diagonal. * For prisms the front (y=0) and back (y=1) face are divided into triangles diff --git a/src/t8_cmesh/t8_cmesh_examples.h b/src/t8_cmesh/t8_cmesh_examples.h index 24b1dcd267..a50362732d 100644 --- a/src/t8_cmesh/t8_cmesh_examples.h +++ b/src/t8_cmesh/t8_cmesh_examples.h @@ -214,6 +214,14 @@ t8_cmesh_new_hypercube_pad_ext (t8_cmesh_t cmesh, const t8_eclass_t eclass, sc_M void t8_cmesh_new_hypercube_hybrid (t8_cmesh_t cmesh, sc_MPI_Comm comm, int periodic); +/** Construct a unit square of two quads and four triangles. + * \param [in,out] cmesh An initialized, but not committed cmesh, as created by \ref t8_cmesh_init. + * Filled and committed in place. + * \param [in] comm The mpi communicator to use. + */ +void +t8_cmesh_new_2D_hypercube_hybrid (t8_cmesh_t cmesh, sc_MPI_Comm comm); + /** Construct a unit interval/square/cube coarse mesh that is periodic in each direction. * Element class? * Hypercube? diff --git a/src/t8_forest/t8_forest_adapt.cxx b/src/t8_forest/t8_forest_adapt.cxx index 1c8e00d382..5ad5dfb1d9 100644 --- a/src/t8_forest/t8_forest_adapt.cxx +++ b/src/t8_forest/t8_forest_adapt.cxx @@ -27,7 +27,9 @@ #include #include #include +#include #include +#include #include /* We want to export the whole implementation to be callable from "C" */ @@ -431,6 +433,10 @@ t8_forest_adapt (t8_forest_t forest) T8_ASSERT (forest->trees->elem_count == forest_from->trees->elem_count); if (forest->set_adapt_recursive) { + if (t8_scheme_has_subelement_scheme (t8_forest_get_scheme (forest_from))) { + SC_CHECK_ABORT (!t8_forest_has_local_subelements (forest_from), + "Recursive adaptation is currently not implemented for subelement schemes."); + } refine_list = sc_list_new (nullptr); } forest->local_num_leaf_elements = 0; @@ -619,6 +625,23 @@ t8_forest_adapt (t8_forest_t forest) } el_considered++; } + else if (refine > 1) { // Subelement case. + T8_ASSERT (t8_eclass_scheme_is_subelement (t8_forest_get_scheme (forest_from), T8_ECLASS_QUAD)); + /* The subelement-callback function returns refine = subelement_type + 1 to avoid subelement_type = 1. + * We undo this (e.g. to use the subelement_type-values that match the binary encoding of the neighbour + * structure for hanging node resolution). + */ + int subelement_type = refine - 1; + + int num_subelements = t8_element_get_number_of_subelements (scheme, tree->eclass, subelement_type); + (void) t8_element_array_push_count (telements, num_subelements); + for (int zz = 0; zz < num_subelements; zz++) { + elements[zz] = t8_element_array_index_locidx_mutable (telements, el_inserted + zz); + } + t8_refine_element_in_subelements (scheme, tree->eclass, elements_from[0], subelement_type, elements); + el_inserted += (t8_locidx_t) num_subelements; + el_considered++; + } else { /* Remove the element */ T8_ASSERT (refine == -2); diff --git a/src/t8_forest/t8_forest_subelement.cxx b/src/t8_forest/t8_forest_subelement.cxx new file mode 100644 index 0000000000..5f50d3e785 --- /dev/null +++ b/src/t8_forest/t8_forest_subelement.cxx @@ -0,0 +1,186 @@ +/* + 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_forest_subelement.cxx + * Implementation of functionality in \ref t8_forest_subelement.hxx. + */ + +#include "t8_forest_subelement.hxx" +#include +#include "t8_forest_general.h" +#include "t8_forest_types.h" +#include "t8_forest_private.h" +#include +#include +#include +#include +#include "t8_forest_adapt.h" + +/** Namespace to hide implementation details.*/ +namespace detail +{ + +/** Adapt callback for \ref t8_forest_discard_subelements. All subelements are coarsened such that the mesh using only + * recursive refinement is restored and the subelements are discarded. This is necessary for another adaption cycle. # + * \param [in] forest The forest to which the new elements belong. + * \param [in] forest_from The forest that is adapted. + * \param [in] which_tree The local tree containing \a elements. + * \param [in] tree_class The eclass of \a which_tree. + * \param [in] lelement_id The local element id in \a forest_from in the tree of the current element. + * \param [in] scheme The scheme of the forest. + * \param [in] is_family If 1, the first \a num_elements entries in \a elements form a family. If 0, they do not. + * \param [in] num_elements The number of entries in \a elements that are defined + * \param [in] elements Pointers to a family or, if \a is_family is zero, pointer to one element. + * \return -1 for subelements, 0 else. + */ +int +discard_subelements_callback ([[maybe_unused]] t8_forest_t forest, [[maybe_unused]] t8_forest_t forest_from, + [[maybe_unused]] t8_locidx_t which_tree, t8_eclass_t tree_class, + [[maybe_unused]] t8_locidx_t lelement_id, const t8_scheme *scheme, + [[maybe_unused]] const int is_family, [[maybe_unused]] const int num_elements, + t8_element_t *elements[]) +{ + // Coarsen if the element is a subelement. + if (t8_element_is_subelement (scheme, tree_class, elements[0]) && is_family) { + return -1; + } + return 0; +} + +/** Adapt callback for hanging node resolution. + * We use the face enumeration to determine which subelement type to use for the transition cell. + * Every face has a flag parameter, which is set to 1, if there is a neighbour with a higher level + * and to 0, if the level of the neighbour is at most the level of the element. + * If all faces are hanging, we use the normal 1:8 refinement and return 1. + * Otherwise, we use subelements and add 1 to every type, to avoid refine = 1. + * \param [in] forest The forest to which the new elements belong. + * \param [in] forest_from The forest that is adapted. + * \param [in] which_tree The local tree containing \a elements. + * \param [in] tree_class The eclass of \a which_tree. + * \param [in] lelement_id The local element id in \a forest_from in the tree of the current element. + * \param [in] scheme The scheme of the forest. + * \param [in] is_family If 1, the first \a num_elements entries in \a elements form a family. If 0, they do not. + * \param [in] num_elements The number of entries in \a elements that are defined + * \param [in] elements Pointers to a family or, if \a is_family is zero, pointer to one element. + * \return The subelement type + 1 to be used for the transition cell, which is a binary encoding of the hanging faces. + */ +int +t8_remove_hanging_nodes_callback ([[maybe_unused]] t8_forest_t forest, t8_forest_t forest_from, t8_locidx_t which_tree, + [[maybe_unused]] t8_eclass_t tree_class, [[maybe_unused]] t8_locidx_t lelement_id, + const t8_scheme *scheme, [[maybe_unused]] const int is_family, + [[maybe_unused]] const int num_elements, t8_element_t *elements[]) +{ + // Determine the hanging faces of the element. This is stored in the subelement type. + int subelement_type = 0; + const int num_faces = scheme->element_get_num_faces (tree_class, elements[0]); + for (int iface = 0; iface < num_faces; iface++) { + const t8_element_t **neighbors; /**< Neighboring elements. */ + int *dual_faces_internal; /**< Face indices of the neighbor elements. */ + int num_neighbors; /**< Number of neighboring elements. */ + t8_locidx_t *neighids; /**< Neighboring elements ids. */ + t8_eclass_t neigh_class; /**< Neighboring elements tree class. */ + + t8_forest_leaf_face_neighbors (forest_from, which_tree, elements[0], &neighbors, iface, &dual_faces_internal, + &num_neighbors, &neighids, &neigh_class); + if (num_neighbors > 1) { + /* Store in correct cell of the binary format. We encode it as f0 -> bit (num_faces-1), ..., f_{n-1} -> bit 0. + * This means (f0 f1 ... f_{n-1}). */ + subelement_type += 1 << ((num_faces - 1) - iface); + } + + // Free allocated memory. + if (num_neighbors > 0) { + T8_FREE (neighbors); + T8_FREE (dual_faces_internal); + T8_FREE (neighids); + } + } + + /* Returning the correct subelement type. */ + if (subelement_type == 0) { /* In this case, there are no hanging faces and we do nothing. */ + return 0; + } + else if (subelement_type == 15) { /* Normal 1:8 refinement. */ + return 1; + } + else { /* Use subelements and add 1 to every type, to avoid refine = 1. */ + return subelement_type + 1; + } +} + +} // namespace detail + +t8_forest_t +t8_forest_remove_hanging_nodes (t8_forest_t forest) +{ + t8_global_productionf ("Into t8_forest_remove_hanging_nodes.\n"); + forest = t8_forest_new_adapt (forest, detail::t8_remove_hanging_nodes_callback, 0, 0, NULL); + t8_global_productionf ("Done t8_forest_remove_hanging_nodes.\n"); + return forest; +} + +t8_forest_t +t8_forest_discard_subelements (t8_forest_t forest) +{ + if (!t8_forest_has_local_subelements (forest)) { + return forest; + } + return t8_forest_new_adapt (forest, detail::discard_subelements_callback, 0, 0, NULL); +} + +bool +t8_forest_has_local_subelements (const t8_forest_t forest) +{ + auto scheme = t8_forest_get_scheme (forest); + if (!t8_scheme_has_subelement_scheme (scheme)) { + return false; + } + for (t8_locidx_t itree = 0; itree < t8_forest_get_num_local_trees (forest); ++itree) { + auto eclass = t8_forest_get_eclass (forest, itree); + if (!t8_eclass_scheme_is_subelement (scheme, eclass)) { + continue; + } + for (t8_locidx_t ielem = 0; ielem < t8_forest_get_tree_num_leaf_elements (forest, itree); ++ielem) { + const t8_element_t *elem = t8_forest_get_leaf_element_in_tree (forest, itree, ielem); + if (t8_element_is_subelement (scheme, eclass, elem)) { + return true; + } + } + } + return false; +} + +bool +t8_forest_has_global_subelements (const t8_forest_t forest) +{ + /* Extract the MPI communicator from the forest */ + sc_MPI_Comm comm = t8_forest_get_mpicomm (forest); + + /* Convert boolean condition to MPI-compatible integer */ + int local = t8_forest_has_local_subelements (forest) ? 1 : 0; + int global = 0; + + const int mpiret = sc_MPI_Allreduce (&local, &global, 1, sc_MPI_INT, sc_MPI_LOR, comm); + SC_CHECK_MPI (mpiret); + + return global != 0; +} diff --git a/src/t8_forest/t8_forest_subelement.hxx b/src/t8_forest/t8_forest_subelement.hxx new file mode 100644 index 0000000000..1222a290c1 --- /dev/null +++ b/src/t8_forest/t8_forest_subelement.hxx @@ -0,0 +1,63 @@ +/* + This file is part of t8code. + t8code is a C library to manage a collection (a forest) of multiple + connected adaptive space-trees of general element classes in parallel. + + Copyright (C) 2025 the developers + + t8code is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + t8code is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with t8code; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +*/ + +/** \file t8_forest_subelement.hxx + * Functionality to handle subelements in a forest. + */ +#pragma once + +#include +#include "t8_forest_general.h" +#include +#include +#include + +/** Remove hanging nodes from the forest by transitioning elements with hanging nodes into subelements. + * \param [in] forest The input forest, which may contain hanging nodes. + * \a forest must be committed before calling this function. Please note that the scheme provided with the + * forest has to be a fitting subelement scheme. + * \return A new forest with the same number of trees and the same connectivity, but conformal without hanging nodes. + */ +t8_forest_t +t8_forest_remove_hanging_nodes (t8_forest_t forest); + +/** Remove all subelements from a forest. This is required to restore the original mesh using only recursive refinement + * and to be able to adapt again. + * \param [in] forest The input forest which may contain subelements. + * \return A new forest with the same number of trees and the same connectivity, but without subelements. + */ +t8_forest_t +t8_forest_discard_subelements (t8_forest_t forest); + +/** Check if a forest contains subelements locally. + * \param [in] forest The forest to be checked. + * \return true if there are subelements in the forest, false otherwise. + */ +bool +t8_forest_has_local_subelements (const t8_forest_t forest); + +/** Check if a forest contains subelements globally. + * \param [in] forest The forest to be checked. + * \return true if there are subelements in the forest, false otherwise. + */ +bool +t8_forest_has_global_subelements (const t8_forest_t forest); diff --git a/src/t8_schemes/t8_scheme.cxx b/src/t8_schemes/t8_scheme.cxx index 8bdbaeae6f..6457bc91d0 100644 --- a/src/t8_schemes/t8_scheme.cxx +++ b/src/t8_schemes/t8_scheme.cxx @@ -24,10 +24,12 @@ * Implements functions declared in \ref t8_scheme.h. */ +#include "t8.h" #include #include #include #include +#include void t8_scheme_ref (t8_scheme_c *scheme) @@ -435,3 +437,30 @@ t8_element_MPI_Unpack (const t8_scheme_c *scheme, const t8_eclass_t tree_class, { return scheme->element_MPI_Unpack (tree_class, recvbuf, buffer_size, position, elements, count, comm); } + +int +t8_element_is_subelement (const t8_scheme_c *scheme, const t8_eclass_t tree_class, const t8_element_t *elem) +{ + SC_CHECK_ABORT (t8_eclass_scheme_is_subelement (scheme, tree_class), + "t8_element_is_subelement was called for a scheme or eclass that does not support subelements.\n"); + return scheme->element_is_subelement (tree_class, elem); +} + +int +t8_element_get_number_of_subelements (const t8_scheme_c *scheme, const t8_eclass_t tree_class, int subelement_type) +{ + SC_CHECK_ABORT ( + t8_eclass_scheme_is_subelement (scheme, tree_class), + "t8_element_get_number_of_subelements was called for a scheme or eclass that does not support subelements.\n"); + return scheme->element_get_number_of_subelements (tree_class, subelement_type); +} + +void +t8_refine_element_in_subelements (const t8_scheme_c *scheme, const t8_eclass_t tree_class, const t8_element_t *elem, + int type, t8_element_t *c[]) +{ + SC_CHECK_ABORT ( + t8_eclass_scheme_is_subelement (scheme, tree_class), + "t8_refine_element_in_subelements was called for a scheme or eclass that does not support subelements.\n"); + scheme->refine_element_in_subelements (tree_class, elem, type, c); +} diff --git a/src/t8_schemes/t8_scheme.h b/src/t8_schemes/t8_scheme.h index 0c3d16ef49..1c2a4feb6c 100644 --- a/src/t8_schemes/t8_scheme.h +++ b/src/t8_schemes/t8_scheme.h @@ -860,6 +860,34 @@ void t8_element_MPI_Unpack (const t8_scheme_c *scheme, const t8_eclass_t tree_class, void *recvbuf, const int buffer_size, int *position, t8_element_t **elements, const unsigned int count, sc_MPI_Comm comm); +/** Check if \a elem is a subelement. + * \param [in] scheme The scheme of the forest. + * \param [in] tree_class The eclass of the current tree. + * \param [in] elem The elem to be checked. + */ +int +t8_element_is_subelement (const t8_scheme_c *scheme, const t8_eclass_t tree_class, const t8_element_t *elem); + +/** Get the number of subelements an element is refined into for a specific type. + * \param [in] scheme The scheme of the forest. + * \param [in] tree_class The eclass of the current tree. + * \param [in] subelement_type The subelement type used for refinement. + */ +int +t8_element_get_number_of_subelements (const t8_scheme_c *scheme, const t8_eclass_t tree_class, int subelement_type); + +/** This defines how an element is refined in subelements using a specified subelement type. + * \param [in] scheme The scheme of the forest. + * \param [in] tree_class The eclass of the current tree. + * \param [in] elem The element to be refined. + * \param [in] type The subelement type to be used for refinement. + * \param [in, out] c An array of allocated elements that will be filled with the subelements of \a elem. + * The number of subelements is determined by \ref t8_element_get_number_of_subelements. + */ +void +t8_refine_element_in_subelements (const t8_scheme_c *scheme, const t8_eclass_t tree_class, const t8_element_t *elem, + int type, t8_element_t *c[]); + T8_EXTERN_C_END (); #endif /* !T8_SCHEME_H */ diff --git a/src/t8_schemes/t8_scheme.hxx b/src/t8_schemes/t8_scheme.hxx index fdaf9510c4..593d5c5960 100644 --- a/src/t8_schemes/t8_scheme.hxx +++ b/src/t8_schemes/t8_scheme.hxx @@ -43,6 +43,9 @@ #include #include #include +#include +#include +#include #include #if T8_ENABLE_DEBUG // Only needed for t8_debug_print_type @@ -96,10 +99,14 @@ struct t8_scheme t8_default_scheme_tet, t8_default_scheme_prism, t8_default_scheme_pyramid, + /* Standalone schemes */ t8_standalone_scheme, t8_standalone_scheme, t8_standalone_scheme, - t8_standalone_scheme + t8_standalone_scheme, + /* Subelement schemes */ + t8_subelementquad_scheme, + t8_subelementtri_scheme >; /* clang-format on */ @@ -1188,6 +1195,67 @@ struct t8_scheme [&] (auto &&scheme) { return scheme.element_MPI_Unpack (recvbuf, buffer_size, position, elements, count, comm); }, eclass_schemes[tree_class]); }; + + /** Check if \a elem is a subelement. + * \param [in] tree_class The eclass of the current tree. + * \param [in] elem The elem to be checked. + */ + inline int + element_is_subelement (const t8_eclass_t tree_class, const t8_element_t *elem) const + { + return std::visit ( + [&] (auto &&scheme) -> int { + if constexpr (requires { scheme.element_is_subelement (elem); }) { + return scheme.element_is_subelement (elem); + } + else { + SC_ABORT ("element_is_subelement not supported by this scheme"); + } + }, + eclass_schemes[tree_class]); + }; + + /** Get the number of subelements an element is refined into for a specific type. + * \param [in] tree_class The eclass of the current tree. + * \param [in] subelement_type The subelement type used for refinement. + */ + inline int + element_get_number_of_subelements (const t8_eclass_t tree_class, int subelement_type) const + { + return std::visit ( + [&] (auto &&scheme) -> int { + if constexpr (requires { scheme.element_get_number_of_subelements (subelement_type); }) { + return scheme.element_get_number_of_subelements (subelement_type); + } + else { + SC_ABORT ("element_get_number_of_subelements not supported by this scheme"); + } + }, + eclass_schemes[tree_class]); + } + + /** This defines how an element is refined in subelements using a specified subelement type. + * \param [in] tree_class The eclass of the current tree. + * \param [in] elem The element to be refined. + * \param [in] type The subelement type to be used for refinement. + * \param [in, out] c An array of allocated elements that will be filled with the subelements of \a elem. + * The number of subelements is determined by \ref element_get_number_of_subelements. + */ + inline void + refine_element_in_subelements (const t8_eclass_t tree_class, const t8_element_t *elem, int type, + t8_element_t *c[]) const + { + std::visit ( + [&] (auto &&scheme) -> void { + if constexpr (requires { scheme.refine_element_in_subelements (elem, type, c); }) { + scheme.refine_element_in_subelements (elem, type, c); + } + else { + SC_ABORT ("refine_element_in_subelements not supported by this scheme"); + } + }, + eclass_schemes[tree_class]); + } }; #endif /* !T8_SCHEME_HXX */ diff --git a/src/t8_schemes/t8_subelement/README.md b/src/t8_schemes/t8_subelement/README.md new file mode 100644 index 0000000000..5d83b92f70 --- /dev/null +++ b/src/t8_schemes/t8_subelement/README.md @@ -0,0 +1,25 @@ +# t8_subelement + +This folder provides **subelement schemes**. Subelements are inserted *after* the standard recursive refinement and enable one additional refinement level that uses a different scheme. + +This is useful, for example, to resolve hanging nodes left behind by the recursive refinement, or to add a uniform subgrid to each mesh element after adaptation, which is beneficial on GPUs. + +Subelements should **never be refined any further**: they always sit at the very bottom of the refinement tree. This keeps the efficient forest-of-trees strategy intact, so that we only store the leaf elements and can recreate the whole forest from them. Before the next adaptation cycle, the subelements are removed again. This also ensures that the parent mesh bounds the mesh quality. + +![](https://github.com/user-attachments/assets/999bc7d9-5617-4dde-bb2a-abb22b921fc7) + +## Implementation details and file structure + +The scheme is built from a **common base** that provides the logic shared by all subelement schemes, plus one **specialization per element class** that supplies the parts that differ between element types. + +- [t8_subelement.hxx](./t8_subelement.hxx) / [t8_subelement.cxx](./t8_subelement.cxx): Main access point to the subelement schemes. Provides the constructor `t8_scheme_new_subelement()`, which assembles the full scheme for all element classes of t8code using the subelement schemes for element classes that are already implemented and for all other standalone/default implementations. + +- [t8_subelement_type.hxx](./t8_subelement_type.hxx): Defines the element class of a subelement. A subelement always consists of an underlying element plus a subelement **type** and **id** that define how the underlying element is transitioned into a subelement. + +- [t8_subelement_scheme.hxx](./t8_subelement_scheme.hxx): The common scheme (`t8_subelement_scheme_common`) implementing the functionality shared by all subelement schemes: construction and destruction, the element memory pool, element sizing, and the general element interface. It is templated on the underlying element class and on a specialization scheme; whenever logic is needed that is *not* identical for all subelements, it delegates to that specialization. + +- [t8_subelement_traits.hxx](./t8_subelement_traits.hxx): Trait definitions that map each concrete subelement scheme to its underlying scheme and subelement type. For example, quadrilateral subelements build on the standalone quad scheme, while triangular subelements build on the default triangle scheme. This is needed for the common subelement scheme implementation. + +- [specializations/](./specializations): Per–element-class specializations providing the subelement logic that is *not* shared by the common scheme: + - [t8_scheme_quads.hxx](./specializations/t8_scheme_quads.hxx): `t8_subelementquad_scheme`, the subelement scheme to resolve hanging nodes for quadrilateral elements. A quad is transitioned into triangular subelements; the subelement type is a binary code over the four faces indicating which of them are hanging. + - [t8_scheme_tri.hxx](./specializations/t8_scheme_tri.hxx): `t8_subelementtri_scheme`, the subelement scheme for triangular elements. diff --git a/src/t8_schemes/t8_subelement/specializations/t8_scheme_quads.hxx b/src/t8_schemes/t8_subelement/specializations/t8_scheme_quads.hxx new file mode 100644 index 0000000000..8cb3872fd7 --- /dev/null +++ b/src/t8_schemes/t8_subelement/specializations/t8_scheme_quads.hxx @@ -0,0 +1,375 @@ +/* + 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_scheme_quads.hxx + * Subelement scheme specialization for quadrilateral elements. A quad is transitioned into + * triangular subelements (e.g. to resolve hanging nodes). The subelement type is a binary + * code over the quad's four faces indicating which of them are hanging; type 0 means the + * element is not a subelement (it is just the underlying standalone quad). This file only + * implements the quad-specific logic. The functionality shared by all subelement schemes + * lives in \ref t8_subelement_scheme.hxx. + */ + +#pragma once +#include +#include "../t8_subelement_scheme.hxx" +#include + +/** Maximum subelement type. The subelement type ranges from 0 (=no subelement, normal standalone quad) to 14. +* The type 15 would mean in binary representation that all faces are hanging but in this case, the element just get +* refined by the standard recursive refinement in 4 standalone quads.*/ +#define T8_SUB_QUAD_MAX_SUBELEMENT_TYPE 14 + +/** Subelement scheme for quadrilateral elements. + * A quad is transitioned into triangular subelements. The subelement type encodes which of + * the quad's four faces are hanging as a binary code (one bit per face). Type 0 means the + * element is not a subelement and is just the underlying standalone quad. For hanging-node + * resolution the type determines the number of subelements the quad is split into, and the + * subelement id runs from 0 to num_subelement - 1. Valid types are 1..14; the all-faces- + * hanging code 15 is excluded. + * + * \verbatim + f3 1 + x - - x - - x x - - x - - x + | | | \ | / | + | | | \ | / | | f0 | f1 | f2 | f3 | + f0 x | f1 --> 1 x - - x | 0 | 1 | 0 | 0 | 1 | = 9 + | | | / \ | + | elem | | / \ | binary code following the face + x - - - - - x x - - - - - x enumeration (here faces f0 and f3 + f2 0 are hanging) + * \endverbatim + * Also have a look at \a vertex_coords_of_subelement for the definition of the subelement ids for quads and the + * order of vertices. + */ +struct t8_subelementquad_scheme: public t8_subelement_scheme_common +{ + public: + /** The recursive scheme used for the underlying (standalone) quad elements. Whenever the + * subelement logic is not needed, the scheme forwards to this underlying scheme. */ + using TUnderlyingScheme = typename t8_subelement_traits::UnderlyingScheme; + /** The subelement element type (an underlying element plus a subelement type and id). */ + using TSubelementType = typename t8_subelement_traits::SubelementType; + + TUnderlyingScheme underlying_scheme {}; /**< Instance of the underlying standalone scheme. */ + + /** Compute the number of corners of an element. + * \param [in] elem The subelement. + * \return The number of corners of \a elem. + */ + static int + subelement_get_num_corners ([[maybe_unused]] const TSubelementType *elem) noexcept + { + return T8_ELEMENT_NUM_CORNERS[T8_ECLASS_TRIANGLE]; + } + + /** Compute the number of faces of a given element. + * \param [in] elem The element. + * \return The number of faces of \a elem. + */ + static int + subelement_get_num_faces ([[maybe_unused]] const TSubelementType *elem) noexcept + { + return T8_ELEMENT_NUM_FACES[T8_ECLASS_TRIANGLE]; + } + + /** Compute the maximum number of faces of a given element and all of its descendants. + * \param [in] elem The element. + * \return The maximum number of faces of \a elem and its descendants. + */ + static int + subelement_get_max_num_faces (const TSubelementType *elem) noexcept + { + return subelement_get_num_faces (elem); + } + + /** Return the shape of an allocated element. + * \param [in] elem The element to be considered + * \return The shape of the element as an eclass + */ + static t8_element_shape_t + subelement_get_shape ([[maybe_unused]] const TSubelementType *elem) noexcept + { + return T8_ECLASS_TRIANGLE; + } + + /** Compute the shape of the face of an element. + * \param [in] elem The element. + * \param [in] face A face of \a elem. + * \return The element shape of the face. As we are in 2D, here always LINE. + */ + static t8_element_shape_t + subelement_get_face_shape ([[maybe_unused]] const TSubelementType *elem, [[maybe_unused]] const int face) noexcept + { + return T8_ECLASS_LINE; + } + + /** Return the max number of children if an element is refined into subelements. + * \return The maximum number of subelements for the quad. + */ + static int + subelement_get_max_num_children () noexcept + { + return 7; + } + + /** Return the number of valid subelement types for a quad. + * \return The maximum valid subelement type. Subelement types run from 0 to this value. + */ + static int + subelement_get_number_of_valid_types () noexcept + { + return T8_SUB_QUAD_MAX_SUBELEMENT_TYPE; + } + + /** Get the number of subelements an element is refined into for a specific type. + * \param [in] subelement_type The subelement type used for refinement. + * \return The number of subelements the quad is split into for \a subelement_type. + */ + static int + element_get_number_of_subelements (int subelement_type) + { + int num_hanging_faces = 0; + /* Count the number of ones of the binary subelement type. This number equals the number of hanging faces. */ + for (int i = 0; i < T8_ELEMENT_NUM_FACES[T8_ECLASS_QUAD]; ++i) { + num_hanging_faces += (subelement_type & (1 << i)) >> i; + } + return T8_ELEMENT_NUM_FACES[T8_ECLASS_QUAD] + num_hanging_faces; + } + + /** This defines how an element is refined into subelements using a specified subelement type. + * \param [in] elem The element to be refined. + * \param [in] type The subelement type to be used for refinement. This is a binary encoding of the hanging faces. + * \param [in, out] c An array of allocated elements that will be filled with the subelements of \a elem. + * The number of subelements is determined by \ref element_get_number_of_subelements. + * \note The different subelement types (up to rotation) are: + * \verbatim + x - - - - - - x x - - - - - x x - - - - - x x - - - - - x x - - x - - x + | | | \ 2 / | | \ / | | \ / | | \ | / | + | | | 1 \ / | | \ / | | \ / | | \ | / | + | | --> x - - X 3 | or x - - x | or x - - x - - x or x - - x - - x + | | | 0 / \ | | / | \ | | / \ | | / \ | + | elem | | / 4 \ | | / | \ | | / \ | | / \ | + + - - - - - - x x - - - - - x x - - x - - x x - - - - - x x - - - - - x + * \endverbatim + * Subelement ids are counted clockwise, starting with the (lower) left subelement with id 0. + * Note that we do not change the underlying quadrant. + */ + void + refine_element_in_subelements (const t8_element_t *elem, int type, t8_element_t *c[]) const noexcept + { + const TSubelementType *element = this->as_subelement (elem); + TSubelementType **subelements = reinterpret_cast (c); + const int num_subelements = this->element_get_number_of_subelements (type); + + T8_ASSERT (type >= 1 && type <= T8_SUB_QUAD_MAX_SUBELEMENT_TYPE); + T8_ASSERT (!this->element_is_subelement (elem)); + T8_ASSERT (this->element_is_valid (elem)); +#if T8_ENABLE_DEBUG + { + for (int j = 0; j < num_subelements; j++) { + T8_ASSERT (this->element_is_valid (c[j])); + } + } +#endif + + /* Setting the parameter values for different subelements. */ + for (int sub_id_counter = 0; sub_id_counter < num_subelements; sub_id_counter++) { + TUnderlyingScheme::element_copy (this->subelement_to_standalone (element), + this->subelement_to_standalone (subelements[sub_id_counter])); + subelements[sub_id_counter]->subelement_type = type; + subelements[sub_id_counter]->subelement_id = sub_id_counter; + T8_ASSERT (this->element_is_valid (c[sub_id_counter])); + } + } + + /** Convert a point in the reference space of a (triangular) subelement to a point in the + * reference space of the tree. + * \param [in] elem The subelement. + * \param [in] ref_coords The coordinates in \f$ [0,1]^2 \f$ of the points in the subelement's reference space. + * \param [in] num_coords The number of points to convert. + * \param [out] out_coords The coordinates of the points in the reference space of the tree. + */ + void + subelement_get_reference_coords (const t8_element_t *elem, const double *ref_coords, const size_t num_coords, + double *out_coords) const noexcept + { + + /* Get the 3 integer vertex coords of the subelement triangle. */ + std::array, 3> vertex_coords; + vertex_coords_of_subelement (elem, vertex_coords); + + /* Normalize to [0,1] by dividing by root length. */ + const double root_len = (1 << T8_ELEMENT_MAXLEVEL[T8_ECLASS_QUAD]); + const double n0[2] = { vertex_coords[0][0] / root_len, vertex_coords[0][1] / root_len }; + const double n1[2] = { vertex_coords[1][0] / root_len, vertex_coords[1][1] / root_len }; + const double n2[2] = { vertex_coords[2][0] / root_len, vertex_coords[2][1] / root_len }; + + for (size_t coord = 0; coord < num_coords; ++coord) { + const double u = ref_coords[coord * 2 + 0]; + const double v = ref_coords[coord * 2 + 1]; + + /* Mapping: (0,0) -> n0, (1,0) -> n1, (1,1) -> n2. */ + out_coords[coord * 2 + 0] = (1.0 - u) * n0[0] + (u - v) * n1[0] + v * n2[0]; + out_coords[coord * 2 + 1] = (1.0 - u) * n0[1] + (u - v) * n1[1] + v * n2[1]; + } + } + + private: + /** Check whether a given face of the parent quad is hanging (and therefore split in half). + * \param [in] type The subelement type (binary code over the faces, order is (f0 ,..., f_{numfaces-1})). + * \param [in] iface The face to check. + * \return True if \a iface is hanging for \a type. + */ + static bool + face_is_split (const unsigned type, const int iface) noexcept + { + return ((type >> ((T8_ELEMENT_NUM_FACES[T8_ECLASS_QUAD] - 1) - iface)) & 1u) != 0u; + } + + /** For each parent face, its two vertices in clockwise order. */ + static constexpr int face_to_clockwise_vertex[4][2] = { { 0, 2 }, { 3, 1 }, { 1, 0 }, { 2, 3 } }; + + /** Compute the integer coordinates of all three vertices of a triangular subelement. + * We use the following order of subelements in a quad: + * Subelement ids are counted clockwise, starting with the (lower) left subelement with id 0. + * The vertices are enumerated clockwise, starting at the center of the transition cell. + * Therefore vertex 0 is always the center of the transition cell. + * For example: + * \verbatim + * f3 V1 + * x - - - - - x x + * | \ 2 / | / | + * | 1 \ / 3 | / 3 | + * f0 x - - + - - x f1 --> + - - x + * | 0 / | \ 4 | V0 V2 + * | / 6 | 5 \ | + * x - - x - - x + * f2 + * \endverbatim + * \param [in] elem The subelement. + * \param [out] vertex_coords The three (x, y) integer vertex coordinates of the subelement. + */ + void + vertex_coords_of_subelement (const t8_element_t *elem, + std::array, 3> &vertex_coords) const noexcept + { + T8_ASSERT (this->element_is_valid (elem)); + T8_ASSERT (this->element_is_subelement (elem)); + const auto *subelement = this->as_subelement (elem); + + /* The length of the parent quadrant and its lower left corner. */ + const int len = this->parent_element_get_len (subelement); + const int origin[2] = { subelement->element.coords[0], subelement->element.coords[1] }; + + // Fill location information. + const std::array location = element_get_location_of_subelement (elem); + const int face_number = location[0]; + const int split = location[1]; + const int sub_face_id = location[2]; + /* Check, whether the get_location function provides meaningful location data. */ + T8_ASSERT (face_number == 0 || face_number == 1 || face_number == 2 || face_number == 3); + T8_ASSERT ((split == 0 && sub_face_id == 0) || (split == 1 && (sub_face_id == 0 || sub_face_id == 1))); + + /** The vertex offsets of a quad (as multiples of its edge length). */ + static constexpr int vertex_offset[4][2] = { { 0, 0 }, { 1, 0 }, { 0, 1 }, { 1, 1 } }; + /** Function lambda to get the vertex coordinates of the parent element. */ + const auto vertex_coords_parent = [&] (const int vertex) { + return std::array { origin[0] + len * vertex_offset[vertex][0], + origin[1] + len * vertex_offset[vertex][1] }; + }; + /** Function lambda to get the midpoint of a face of the parent element. */ + const auto vertex_midpoint_coords_parent = [&] (const int face) { + const int face_vertex1 = face_to_clockwise_vertex[face][0]; + const int face_vertex2 = face_to_clockwise_vertex[face][1]; + return std::array { + origin[0] + (len * (vertex_offset[face_vertex1][0] + vertex_offset[face_vertex2][0])) / 2, + origin[1] + (len * (vertex_offset[face_vertex1][1] + vertex_offset[face_vertex2][1])) / 2 + }; + }; + + /* Vertex 0 is always the centre of the transition cell. */ + vertex_coords[0] = { origin[0] + len / 2, origin[1] + len / 2 }; + /* Vertices 1 and 2 are the face's two clockwise vertices, unless the face is split: then one of + * them is replaced by the face midpoint. */ + const std::array vertex_start = vertex_coords_parent (face_to_clockwise_vertex[face_number][0]); + const std::array vertex_end = vertex_coords_parent (face_to_clockwise_vertex[face_number][1]); + const std::array face_midpoint = vertex_midpoint_coords_parent (face_number); + + vertex_coords[1] = (split && sub_face_id) ? face_midpoint : vertex_start; + vertex_coords[2] = (split && !sub_face_id) ? face_midpoint : vertex_end; + } + + /** Determine the location of a subelement within its transition cell. + * \param [in] elem The subelement. + * \return Three values: + * - the face of the parent quad the subelement is adjacent to ({0,1,2,3}) + * - whether that face is split in half ({0,1}) + * - and whether it is the first or second subelement at the face ({0,1}). + * + * For a subelement of type 14 the location array is {1,1,0} for id 3 and {2,1,1} for id 6. + * \verbatim + * f3 V1 + * x - - - - - x x + * | \ 2 / | / | + * | 1 \ / 3 | / 3 | + * f0 x - - + - - x f1 --> + - - x + * | 0 / | \ 4 | V0 V2 + * | / 6 | 5 \ | + * x - - x - - x + * f2 + * \endverbatim + */ + std::array + element_get_location_of_subelement (const t8_element_t *elem) const + { + T8_ASSERT (this->element_is_subelement (elem)); + T8_ASSERT (this->element_is_valid (elem)); + const auto *subelement = this->as_subelement (elem); + const unsigned type = static_cast (subelement->subelement_type); + const int sub_id = subelement->subelement_id; + T8_ASSERT (sub_id < element_get_number_of_subelements (static_cast (type))); + /** The parent face at each clockwise position, starting at the left face: left (f0), top (f3), + * right (f1), bottom (f2). Subelement ids are assigned in this order. */ + const int clockwise_ordering_to_parent_face[4] = { 0, 3, 1, 2 }; + + /* Walk the faces in clockwise order (the order in which subelement ids are assigned). Each face + * contributes one subelement, or two if it is split. The subelement lies at the first face whose + * running count exceeds sub_id. */ + int clockwise_face = 0; + int split = 0; + int subelements_up_to = 0; // The current clockwise face iface contains subelements with ids < this number. + for (clockwise_face = 0; clockwise_face < T8_ELEMENT_NUM_FACES[T8_ECLASS_QUAD]; ++clockwise_face) { + split = face_is_split (type, clockwise_ordering_to_parent_face[clockwise_face]); + subelements_up_to += split + 1; + if (sub_id < subelements_up_to) { + break; + } + } // Now split and the clockwise face are set correctly. + /* On a split face the two subelements take the last two ids of its range. Determine which one. + * It is the second subelement if sub_id + 1 == subelements_up_to (as this is the last sub id that is contained in + * this face), otherwise it is the first. (and 0=false if not split). */ + const int sub_face_id = split && (sub_id + 1 == subelements_up_to); + + return { clockwise_ordering_to_parent_face[clockwise_face], split, sub_face_id }; + } +}; diff --git a/src/t8_schemes/t8_subelement/specializations/t8_scheme_tri.hxx b/src/t8_schemes/t8_subelement/specializations/t8_scheme_tri.hxx new file mode 100644 index 0000000000..717256de51 --- /dev/null +++ b/src/t8_schemes/t8_subelement/specializations/t8_scheme_tri.hxx @@ -0,0 +1,335 @@ +/* + 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_scheme_tri.hxx + * Subelement scheme specialization for triangular elements. A triangle is transitioned into + * triangular subelements (e.g. to resolve hanging nodes). The subelement type is a binary code + * over the triangle's three faces indicating which of them are hanging; type 0 means the element + * is not a subelement (it is just the underlying standalone triangle). This file implements only + * the triangle-specific logic. The functionality shared by all subelement schemes lives in + * t8_subelement_scheme.hxx. + */ + +#pragma once + +#include +#include +#include "../t8_subelement_scheme.hxx" +#include +#include + +/** Maximum subelement type. The subelement type ranges from 0 (= no subelement, normal standalone + * triangle) to 6. Type 7 would mean in binary representation that all faces are hanging, but in + * that case the element is refined by the standard recursive refinement instead. */ +#define T8_TRI_MAX_SUBELEMENT_TYPE 6 + +/** Subelement scheme for triangular elements. + * A triangle is transitioned into triangular subelements. The subelement type encodes which of the + * triangle's three faces are hanging as a binary code (one bit per face). Type 0 means the element + * is not a subelement and is just the underlying standalone triangle. Valid types are 1..6; the + * all-faces-hanging code 7 is excluded, since that is a normal recursive refinement. + * + * Please have a look at \a vertex_coords_of_subelement for the definition of the subelement ids for triangles. + */ +struct t8_subelementtri_scheme: public t8_subelement_scheme_common +{ + public: + /** The recursive scheme used for the underlying (standalone) triangle elements. Whenever the + * subelement logic is not needed, the scheme forwards to this underlying scheme. */ + using TUnderlyingScheme = typename t8_subelement_traits::UnderlyingScheme; + /** The subelement element type (an underlying element plus a subelement type and id). */ + using TSubelementType = typename t8_subelement_traits::SubelementType; + + TUnderlyingScheme underlying_scheme {}; /**< Instance of the underlying standalone scheme. */ + + /** Compute the number of corners of an element. + * \param [in] elem The subelement. + * \return The number of corners of \a elem. + */ + static int + subelement_get_num_corners ([[maybe_unused]] const TSubelementType *elem) noexcept + { + return T8_ELEMENT_NUM_CORNERS[T8_ECLASS_TRIANGLE]; + } + + /** Compute the number of faces of a given element. + * \param [in] elem The element. + * \return The number of faces of \a elem. + */ + static int + subelement_get_num_faces ([[maybe_unused]] const TSubelementType *elem) noexcept + { + return T8_ELEMENT_NUM_FACES[T8_ECLASS_TRIANGLE]; + } + + /** Compute the maximum number of faces of a given element and all of its descendants. + * \param [in] elem The element. + * \return The maximum number of faces of \a elem and its descendants. + */ + static int + subelement_get_max_num_faces (const TSubelementType *elem) noexcept + { + return subelement_get_num_faces (elem); + } + + /** Return the shape of an allocated element. + * \param [in] elem The element to be considered + * \return The shape of the element as an eclass + */ + static t8_element_shape_t + subelement_get_shape ([[maybe_unused]] const TSubelementType *elem) noexcept + { + return T8_ECLASS_TRIANGLE; + } + + /** Compute the shape of the face of an element. + * \param [in] elem The element. + * \param [in] face A face of \a elem. + * \return The element shape of the face. As we are in 2D, here always LINE. + */ + static t8_element_shape_t + subelement_get_face_shape ([[maybe_unused]] const TSubelementType *elem, [[maybe_unused]] const int face) noexcept + { + return T8_ECLASS_LINE; + } + + /** Return the max number of children if an element is refined into subelements. + * \return The maximum number of subelements for the triangle. + */ + static int + subelement_get_max_num_children () noexcept + { + return 3; + } + + /** Return the number of valid subelement types for a triangle. + * \return The maximum valid subelement type. Subelement types run from 0 to this value. + */ + static int + subelement_get_number_of_valid_types () noexcept + { + return T8_TRI_MAX_SUBELEMENT_TYPE; + } + + /** Get the number of subelements an element is refined into for a specific type. + * \param [in] subelement_type The subelement type used for refinement. + * \return The number of subelements the triangle is split into (hanging faces + 1). + */ + static int + element_get_number_of_subelements (int subelement_type) + { + int num_hanging_faces = 0; + /* Count the number of ones of the binary subelement type. This number equals the number of hanging faces. */ + for (int i = 0; i < T8_ELEMENT_NUM_FACES[T8_ECLASS_TRIANGLE]; ++i) { + num_hanging_faces += (subelement_type & (1 << i)) >> i; + } + return num_hanging_faces + 1; + } + + /** This defines how an element is refined into subelements using a specified subelement type. + * \param [in] elem The element to be refined. + * \param [in] type The subelement type to be used for refinement. This is a binary encoding of the hanging faces. + * \param [in, out] c An array of allocated elements that will be filled with the subelements of \a elem. + * The number of subelements is determined by \ref element_get_number_of_subelements. + */ + void + refine_element_in_subelements (const t8_element_t *elem, int type, t8_element_t *c[]) const noexcept + { + const TSubelementType *element = this->as_subelement (elem); + TSubelementType **subelements = reinterpret_cast (c); + const int num_subelements = this->element_get_number_of_subelements (type); + + T8_ASSERT (type >= 1 && type <= T8_TRI_MAX_SUBELEMENT_TYPE); + T8_ASSERT (!this->element_is_subelement (elem)); + T8_ASSERT (this->element_is_valid (elem)); +#if T8_ENABLE_DEBUG + { + for (int j = 0; j < num_subelements; j++) { + T8_ASSERT (this->element_is_valid (c[j])); + } + } +#endif + + /* Setting the parameter values for different subelements. */ + for (int sub_id_counter = 0; sub_id_counter < num_subelements; sub_id_counter++) { + underlying_scheme.element_copy (this->subelement_to_standalone (element), + this->subelement_to_standalone (subelements[sub_id_counter])); + subelements[sub_id_counter]->subelement_type = type; + subelements[sub_id_counter]->subelement_id = sub_id_counter; + T8_ASSERT (this->element_is_valid (c[sub_id_counter])); + } + } + + /** Convert points in the reference space of a (triangular) subelement to points in the reference + * space of the tree. + * \param [in] elem The subelement. + * \param [in] ref_coords The coordinates in \f$ [0,1]^2 \f$ of the points in the subelement's reference space. + * \param [in] num_coords The number of points to convert. + * \param [out] out_coords The coordinates of the points in the reference space of the tree. + */ + void + subelement_get_reference_coords (const t8_element_t *elem, const double *ref_coords, const size_t num_coords, + double *out_coords) const noexcept + { + + /* Get the 3 integer vertex coords of the subelement triangle. */ + std::array, 3> vertex_coords; + vertex_coords_of_subelement (elem, vertex_coords); + + /* Normalize to [0,1] by dividing by root length. */ + const double root_len = (1 << T8_ELEMENT_MAXLEVEL[T8_ECLASS_TRIANGLE]); + double n0[2] = { vertex_coords[0][0] / root_len, vertex_coords[0][1] / root_len }; + double n1[2] = { vertex_coords[1][0] / root_len, vertex_coords[1][1] / root_len }; + double n2[2] = { vertex_coords[2][0] / root_len, vertex_coords[2][1] / root_len }; + + for (size_t coord = 0; coord < num_coords; ++coord) { + const double u = ref_coords[coord * 2 + 0]; + const double v = ref_coords[coord * 2 + 1]; + + /* Mapping: (0,0) -> n0, (1,0) -> n1, (1,1) -> n2. */ + out_coords[coord * 2 + 0] = (1.0 - u) * n0[0] + (u - v) * n1[0] + v * n2[0]; + out_coords[coord * 2 + 1] = (1.0 - u) * n0[1] + (u - v) * n1[1] + v * n2[1]; + } + } + + private: + /** Check whether a given face of the parent triangle is hanging. + * \param [in] type The subelement type (binary code over the faces, f0 is the most significant bit). + * \param [in] iface The face to check. + * \return True if \a iface is hanging for \a type. + */ + static bool + face_is_hanging (const unsigned type, const int iface) noexcept + { + // Get the bit corresponding to iface. + // If that bit is 1, the face is hanging. + // 1u is for lowest bit extraction. + return ((type >> ((T8_ELEMENT_NUM_FACES[T8_ECLASS_TRIANGLE] - 1) - iface)) & 1u) != 0u; + } + + /** Compute the integer coordinates of the three vertices of a triangular subelement. + * + * For this, we first define the order of the subelements and the subelement vertices: + * All subelements of a transition cell share one common point, which we define as \a m_c, and each subelement is + * spanned by \a m_c together with two consecutive points of a \b path along the boundary of the + * parent triangle. This defines the numbering completely: + * - The \a main \a face is the lowest-indexed hanging face \a fA. + * - The common point \a m_c is the midpoint of the main face. Every subelement contains it. + * - Let \a v_a < \a v_b be the two end vertices of the main face fA. The \b path walks the parent + * boundary from \a v_a to \a v_b the way that does not traverse the main face fA (so the other way around + * such that we go over all other faces). That walk passes through exactly one other vertex, \a v_c, and + * traverses exactly two faces: the edge (\a v_a, \a v_c) and the edge (\a v_c, \a v_b). The midpoint of + * each traversed face that is hanging is inserted at its position on the walk. + * - The subelement with id \a i is then defined through the vertices: + * vertex 0 = \a m_c, vertex 1 = path[ \a i ], vertex 2 = path[ \a i+1 ]. + * + * Since the path has (number of hanging faces + 2) points, there are (number of hanging faces + 1) + * subelements, which matches \ref element_get_number_of_subelements. No case distinction is needed: + * one hanging face yields a path of three points, two hanging faces a path of four. + * + * \verbatim + f2 hanging f2 hanging + one hanging face (here f2) two hanging faces (here f1 and f2) (not nicely displayed) + + v2 v2 + /| \ / \ + / | \ / 2 \ + f1 / | \ f0 M1 x – \ f0 + / | \ / | \ \ + / 0 | 1 \ / 0 \ –– \ + / | \ / | 1 \\ + v0 ---- M2----- v1 v0 ---- M2------ v1 + f2 f2 + + main face = f2, \a m_c = M2 main face = f1 (lowest), \a m_c = M1 + path: v0 -> v2 -> v1 path: v0 -> M2 -> v1 -> v2 + * \endverbatim + * + * \param [in] elem The subelement. + * \param [out] vertex_coords The three (x, y) integer vertex coordinates of the subelement. + */ + void + vertex_coords_of_subelement (const t8_element_t *elem, + std::array, 3> &vertex_coords) const noexcept + { + T8_ASSERT (this->element_is_valid (elem)); + T8_ASSERT (this->element_is_subelement (elem)); + const auto *subelement = this->as_subelement (elem); + const unsigned type = static_cast (subelement->subelement_type); + const unsigned id = static_cast (subelement->subelement_id); + [[maybe_unused]] const int num_hanging_faces = std::popcount (type); + T8_ASSERT (num_hanging_faces == 1 || num_hanging_faces == 2); + + /* The corners of the parent triangle. */ + std::array, 3> parent_coords; + for (int icorner = 0; icorner < T8_ELEMENT_NUM_CORNERS[T8_ECLASS_TRIANGLE]; ++icorner) { + underlying_scheme.element_get_vertex_integer_coords (this->subelement_to_standalone (subelement), icorner, + parent_coords[icorner].data ()); + } + + /* Lambda for the midpoints of a face of the parent triangle. */ + const auto face_midpoint = [&parent_coords] (const int iface) { + const std::array &first = parent_coords[t8_face_vertex_to_tree_vertex[T8_ECLASS_TRIANGLE][iface][0]]; + const std::array &second = parent_coords[t8_face_vertex_to_tree_vertex[T8_ECLASS_TRIANGLE][iface][1]]; + return std::array { (first[0] + second[0]) / 2, (first[1] + second[1]) / 2 }; + }; + + /* The main face is the lowest-indexed hanging face; m_c is its midpoint. */ + int main_face = 0; + while (!face_is_hanging (type, main_face)) { + ++main_face; + } + T8_ASSERT (main_face < T8_ELEMENT_NUM_FACES[T8_ECLASS_TRIANGLE]); + const std::array m_c = face_midpoint (main_face); + + /* Build the path: Walk the parent edges from the first to the second end vertex of the main face, the way + * that does not traverse the main face itself. + */ + const int start_vertex = t8_face_vertex_to_tree_vertex[T8_ECLASS_TRIANGLE][main_face][0]; + const int end_vertex = t8_face_vertex_to_tree_vertex[T8_ECLASS_TRIANGLE][main_face][1]; + + // The path has maximal length 4 for 2 hanging faces. + // For the path we use the property of the triangle enumeration that the face has always the id of the opposite + // vertex (so the only vertex it is not adjacent to). Therefore we can use the face ids to get the midpoints of the hanging faces. + std::array, 4> path; + int path_length = 0; + path[path_length++] = parent_coords[start_vertex]; + // The next face to traverse is the face opposite to the end vertex. Therefore it has the id "end_vertex". + if (face_is_hanging (type, end_vertex)) { + path[path_length++] = face_midpoint (end_vertex); + } + // Next vertex has the id of the main face. + path[path_length++] = parent_coords[main_face]; + if (face_is_hanging (type, start_vertex)) { + path[path_length++] = face_midpoint (start_vertex); + } + path[path_length++] = parent_coords[end_vertex]; + + /* Path length should be 4 for 2 hanging faces and 3 for 1. */ + T8_ASSERT (path_length == num_hanging_faces + 2); + T8_ASSERT (static_cast (id) + 1 < path_length); + + vertex_coords[0] = m_c; + vertex_coords[1] = path[id]; + vertex_coords[2] = path[id + 1]; + } +}; diff --git a/src/t8_schemes/t8_subelement/t8_subelement.cxx b/src/t8_schemes/t8_subelement/t8_subelement.cxx new file mode 100644 index 0000000000..f3a09c3c51 --- /dev/null +++ b/src/t8_schemes/t8_subelement/t8_subelement.cxx @@ -0,0 +1,72 @@ +/* + This file is part of t8code. + t8code is a C library to manage a collection (a forest) of multiple + connected adaptive space-trees of general element classes in parallel. + + Copyright (C) 2025 the developers + + t8code is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + t8code is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with t8code; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +*/ + +/** \file t8_subelement.cxx + * Implements functions declared in \ref t8_subelement.hxx. + */ + +#include +#include +#include +#include "t8_subelement_scheme.hxx" +#include "specializations/t8_scheme_quads.hxx" +#include "specializations/t8_scheme_tri.hxx" + +const t8_scheme * +t8_scheme_new_subelement (void) +{ + t8_scheme_builder builder; + + builder.add_eclass_scheme> (); + builder.add_eclass_scheme> (); + builder.add_eclass_scheme (); + builder.add_eclass_scheme (); + builder.add_eclass_scheme> (); + builder.add_eclass_scheme (); + builder.add_eclass_scheme (); + builder.add_eclass_scheme (); + return builder.build_scheme (); +} + +int +t8_eclass_scheme_is_subelement (const t8_scheme *scheme, const t8_eclass_t eclass) +{ + switch (eclass) { + case T8_ECLASS_QUAD: + return scheme->check_eclass_scheme_type (T8_ECLASS_QUAD); + case T8_ECLASS_TRIANGLE: + return scheme->check_eclass_scheme_type (T8_ECLASS_TRIANGLE); + default: + return 0; /* Default return value false. */ + } +} + +bool +t8_scheme_has_subelement_scheme (const t8_scheme *scheme) +{ + for (int ieclass = T8_ECLASS_ZERO; ieclass < T8_ECLASS_COUNT; ++ieclass) { + if (t8_eclass_scheme_is_subelement (scheme, static_cast (ieclass))) { + return true; + } + } + return false; +} diff --git a/src/t8_schemes/t8_subelement/t8_subelement.hxx b/src/t8_schemes/t8_subelement/t8_subelement.hxx new file mode 100644 index 0000000000..50ef75ed1f --- /dev/null +++ b/src/t8_schemes/t8_subelement/t8_subelement.hxx @@ -0,0 +1,49 @@ +/* + 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_subelement.hxx + * Define the subelement scheme interface. + */ + +#pragma once + +#include + +/** Return the subelement scheme implementation of t8code. */ +const t8_scheme * +t8_scheme_new_subelement (void); + +/** Check whether a given eclass_scheme is one of the subelement schemes. + * \param [in] scheme A (pointer to a) scheme. + * \param [in] eclass The eclass to check. + * \return True if \a scheme is one of the subelement schemes for the element class, false otherwise. + */ +bool +t8_eclass_scheme_is_subelement (const t8_scheme *scheme, const t8_eclass_t eclass); + +/** Check if \a scheme uses a subelement scheme for any eclass. + * This means that it checks if \ref t8_eclass_scheme_is_subelement is true for any eclass. + * \param [in] scheme A (pointer to a) scheme. + * \return True if \a scheme uses a subelement scheme for any eclass, false otherwise. + */ +bool +t8_scheme_has_subelement_scheme (const t8_scheme *scheme); diff --git a/src/t8_schemes/t8_subelement/t8_subelement_scheme.hxx b/src/t8_schemes/t8_subelement/t8_subelement_scheme.hxx new file mode 100644 index 0000000000..9aaa72af3d --- /dev/null +++ b/src/t8_schemes/t8_subelement/t8_subelement_scheme.hxx @@ -0,0 +1,1252 @@ +/* + 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_subelement_scheme.hxx + * Common functionality for schemes that support subelements. + * Subelements are temporary elements used after common recursive adaptation and are discarded before the next + * adaptation cycle. + * Scheme-specific subelement functionality is provided by the corresponding specialization in the specializations + * folder and this file only implements functionality that is equal for all subelements. + */ +#pragma once + +#include +#include +#include +#include +#include "t8_subelement_traits.hxx" +#include +#include + +/** Scheme for the common functionality of all subelements. + * Subelements are discarded before the next adaptation cycle and do not have children. + * \tparam TEclass The element class of the underlying elements which we want to define subelements for. + * The subelements themselves could have another eclass. + * \tparam TSubelementSchemeSpecialization Specialization scheme for the subelements. Every time we need the subelement + * logic which is not equal for all subelements, the scheme calls the functionality of this subelement scheme. + */ +template +struct t8_subelement_scheme_common: + public t8_scheme_helpers> +{ + public: + /** The subelement type used by this subelement scheme defined by a trait. */ + using TSubelementType = typename t8_subelement_traits::SubelementType; + + /** Constructor. */ + t8_subelement_scheme_common () noexcept + : element_size (sizeof (TSubelementType)), scheme_context (sc_mempool_new (element_size)) {}; + + protected: + size_t element_size; /**< The size in bytes of an element. */ + void *scheme_context; /**< Anonymous implementation context. */ + + public: + // #################################____Constructor & Destructor...____############################################### + /** Destructor. */ + ~t8_subelement_scheme_common () + { + T8_ASSERT (scheme_context != NULL); + SC_ASSERT (((sc_mempool_t *) scheme_context)->elem_count == 0); + sc_mempool_destroy ((sc_mempool_t *) scheme_context); + } + + /** Move constructor */ + t8_subelement_scheme_common (t8_subelement_scheme_common &&other) noexcept + : element_size (other.element_size), scheme_context (std::exchange (other.scheme_context, nullptr)) + { + } + + /** Move assignment operator */ + t8_subelement_scheme_common & + operator= (t8_subelement_scheme_common &&other) noexcept + { + if (this != &other) { + // Free existing resources of moved-to object + if (scheme_context) { + sc_mempool_destroy ((sc_mempool_t *) scheme_context); + } + // Transfer ownership of resources + element_size = other.element_size; + scheme_context = other.scheme_context; + // Leave the source object in a valid state + other.scheme_context = nullptr; + } + return *this; + } + + /** Copy constructor */ + t8_subelement_scheme_common (const t8_subelement_scheme_common &other) + : element_size (other.element_size), scheme_context (sc_mempool_new (other.element_size)) {}; + + /** Copy assignment operator */ + t8_subelement_scheme_common & + operator= (const t8_subelement_scheme_common &other) + { + if (this != &other) { + // Free existing resources of assigned-to object + if (scheme_context) { + sc_mempool_destroy ((sc_mempool_t *) scheme_context); + } + // Copy the values from the source object + element_size = other.element_size; + scheme_context = sc_mempool_new (other.element_size); + } + return *this; + } + + // ################################################____GENERAL INFO____############################################### + + /** Return the size of any element. + * \return The size of an element. + */ + static constexpr size_t + get_element_size (void) noexcept + { + return sizeof (TSubelementType); + } + + /** Returns true, if there is one element in the tree, that does not refine into 2^dim children, false otherwise. + * \return Always true as subelements may occur. + */ + static constexpr int + refines_irregular (void) noexcept + { + return true; // Potentially there are subelements. + } + + /** Return the maximum allowed level for any element of a given class. + * \return The maximum allowed level for elements of class \b ts. + */ + constexpr int + get_maxlevel (void) const noexcept + { + return derived ().underlying_scheme.get_maxlevel () - 1; // We need to reserve one level for the subelements. + } + + // ################################################____SHAPE INFORMATION____########################################## + + /** Compute the number of corners of an element. + * \param [in] elem The element. + * \return The number of corners of \a elem. + */ + int + element_get_num_corners (const t8_element_t *elem) const noexcept + { + T8_ASSERT (element_is_valid (elem)); + if (!element_is_subelement (elem)) { + return derived ().underlying_scheme.element_get_num_corners (element_to_standalone (elem)); + } + return TSubelementSchemeSpecialization::subelement_get_num_corners (as_subelement (elem)); + } + + /** Compute the number of faces of a given element. + * \param [in] elem The element. + * \return The number of faces of \a elem. + */ + int + element_get_num_faces (const t8_element_t *elem) const noexcept + { + T8_ASSERT (element_is_valid (elem)); + if (!element_is_subelement (elem)) { + return derived ().underlying_scheme.element_get_num_faces (element_to_standalone (elem)); + } + return TSubelementSchemeSpecialization::subelement_get_num_faces (as_subelement (elem)); + } + + /** Compute the maximum number of faces of a given element and all of its descendants. + * \param [in] elem The element. + * \return The maximum number of faces of \a elem and its descendants. + */ + int + element_get_max_num_faces (const t8_element_t *elem) const noexcept + { + T8_ASSERT (element_is_valid (elem)); + return std::max (derived ().underlying_scheme.element_get_max_num_faces (element_to_standalone (elem)), + TSubelementSchemeSpecialization::subelement_get_max_num_faces (as_subelement (elem))); + } + + /** Return the shape of an allocated element. + * \param [in] elem The element to be considered + * \return The shape of the element as an eclass + */ + t8_element_shape_t + element_get_shape (const t8_element_t *elem) const noexcept + { + T8_ASSERT (element_is_valid (elem)); + if (!element_is_subelement (elem)) { + return derived ().underlying_scheme.element_get_shape (element_to_standalone (elem)); + } + return TSubelementSchemeSpecialization::subelement_get_shape (as_subelement (elem)); + } + + /** Not implemented for this scheme. + * \param [in] element The element. + * \param [in] face A face index for \a element. + * \param [in] corner A corner index for the face 0 <= \a corner < num_face_corners. + * \return The corner number of the \a corner-th vertex of \a face. + */ + static int + element_get_face_corner ([[maybe_unused]] const t8_element_t *element, [[maybe_unused]] const int face, + [[maybe_unused]] const int corner) noexcept + { + SC_ABORT ("element_get_face_corner is not implemented for subelements yet.\n"); + } + + /** Not implemented for this scheme. + * \param [in] element The element. + * \param [in] corner A corner index for the face. + * \param [in] face A face index for \a corner. + * \return The face number of the \a face-th face at \a corner. + */ + static int + element_get_corner_face ([[maybe_unused]] const t8_element_t *element, [[maybe_unused]] const int corner, + [[maybe_unused]] const int face) noexcept + { + SC_ABORT ("element_get_corner_face is not implemented for subelements yet.\n"); + } + + /** Compute the shape of the face of an element. + * \param [in] elem The element. + * \param [in] face A face of \a elem. + * \return The element shape of the face. + */ + t8_element_shape_t + element_get_face_shape (const t8_element_t *elem, const int face) const noexcept + { + T8_ASSERT (element_is_valid (elem)); + T8_ASSERT (0 <= face && face < element_get_num_faces (elem)); + if (!element_is_subelement (elem)) { + return derived ().underlying_scheme.element_get_face_shape (element_to_standalone (elem), face); + } + return TSubelementSchemeSpecialization::subelement_get_face_shape (as_subelement (elem), face); + } + + /** Return the level of a particular element. For subelements, the level is the same as the level of the parent. + * \param [in] elem The element whose level should be returned. + * \return The level of \b elem. + */ + int + element_get_level (const t8_element_t *elem) const noexcept + { + T8_ASSERT (element_is_valid (elem)); + return derived ().underlying_scheme.element_get_level (element_to_standalone (elem)); + } + + // ################################################____GENERAL HELPER____############################################# + + /** Copy all entries of \b source to \b dest. \b dest must be an existing + * element. No memory is allocated by this function. + * \param [in] source The element whose entries will be copied to \b dest. + * \param [in,out] dest This element's entries will be overwrite with the entries of \b source. + * \note \a source and \a dest may point to the same element. + */ + void + element_copy (const t8_element_t *source, t8_element_t *dest) const noexcept + { + T8_ASSERT (element_is_valid (source)); + if (source == dest) + return; + memcpy (as_subelement (dest), as_subelement (source), sizeof (TSubelementType)); + T8_ASSERT (element_is_valid (dest)); + } + + /** Check if two elements are equal. + * \note For subelements, it is only checked that the type is equal and not the id!! + * \param [in] elem1 The first element. + * \param [in] elem2 The second element. + * \return true if the elements are equal, false if they are not equal + */ + int + element_is_equal (const t8_element_t *elem1, const t8_element_t *elem2) const noexcept + { + T8_ASSERT (element_is_valid (elem1) && element_is_valid (elem2)); + const auto *el1 = as_subelement (elem1); + const auto *el2 = as_subelement (elem2); + if (el1->subelement_type != el2->subelement_type) { + return 0; + } + return derived ().underlying_scheme.element_is_equal (subelement_to_standalone (el1), + subelement_to_standalone (el2)); + } + + // ################################################____REFINEMENT____################################################ + /** Create the root element. + * \param [in,out] elem The element that is filled with the root. + */ + void + set_to_root (t8_element_t *elem) const noexcept + { + auto *subelement = as_subelement (elem); + reset_subelement_values (subelement); + derived ().underlying_scheme.set_to_root (subelement_to_standalone (subelement)); + } + + /** Compute the parent of a given element \b elem and store it in \b parent. + * \b parent needs to be an existing element. No memory is allocated by this function. \b elem and \b parent can + * point to the same element, then the entries of \b elem are overwritten by the ones of its parent. + * \param [in] elem The element whose parent will be computed. + * \param [in,out] parent This element's entries will be overwritten by those of \b elem's parent. + * The storage for this element must exist and match the element class of the parent. + */ + void + element_get_parent (const t8_element_t *elem, t8_element_t *parent) const noexcept + { + T8_ASSERT (element_is_valid (elem)); + const auto *subelement = as_subelement (elem); + auto *parent_subelement = as_subelement (parent); + reset_subelement_values (parent_subelement); + if (element_is_subelement (elem)) { + // For subelements, the parent is the element from which they are refined. + derived ().underlying_scheme.element_copy (subelement_to_standalone (subelement), + subelement_to_standalone (parent_subelement)); + return; + } + derived ().underlying_scheme.element_get_parent (subelement_to_standalone (subelement), + subelement_to_standalone (parent_subelement)); + } + + /** Compute the number of siblings of an element. That is the number of elements with the same parent (if available). + * \param [in] elem The element. + * \return The number of siblings of \a element. + * Note that this number is >= 1, since we count the element itself as a sibling.. + */ + int + element_get_num_siblings (const t8_element_t *elem) const noexcept + { + T8_ASSERT (element_is_valid (elem)); + if (!element_is_subelement (elem)) { + return derived ().underlying_scheme.element_get_num_siblings (element_to_standalone (elem)); + } + return element_get_number_of_subelements (as_subelement (elem)->subelement_type); + } + + /** Not implemented for this scheme + * \param [in] elem The element whose sibling will be computed. + * \param [in] sibid The id of the sibling computed. + * \param [in,out] sibling This element's entries will be overwritten by those of \b elem's sibid-th sibling. + */ + static void + element_get_sibling ([[maybe_unused]] const t8_element_t *elem, [[maybe_unused]] const int sibid, + [[maybe_unused]] t8_element_t *sibling) noexcept + { + SC_ABORT ("element_get_sibling not implemented yet.\n"); + } + + /** As subelements are discarded before the next adaptation cycle, they do not have children. + * \param [in] elem This must be a valid element, bigger than maxlevel. + * \param [in] childid The number of the child to construct. + * \param [in,out] child The storage for this element must exist and match the element class of the child. + */ + void + element_get_child (const t8_element_t *elem, const int childid, t8_element_t *child) const noexcept + { + SC_CHECK_ABORT (!element_is_subelement (elem), "element_get_child: Cannot construct child of a subelement.\n"); + derived ().underlying_scheme.element_get_child (element_to_standalone (elem), childid, + element_to_standalone (child)); + } + + /** Return the number of children of an element when it is refined. Not for subelements as they do not have children. + * \param [in] elem The element whose number of children is returned. + * \return The number of children of \a elem if it is to be refined. + */ + int + element_get_num_children (const t8_element_t *elem) const noexcept + { + SC_CHECK_ABORT (!element_is_subelement (elem), + "element_get_num_children: Cannot construct child of a subelement.\n"); + return derived ().underlying_scheme.element_get_num_children (element_to_standalone (elem)); + } + + /** Return the max number of children of an eclass. + * \return Maximum number of possible children (maximum of normal refinement and subelement refinement). + */ + int + get_max_num_children () const noexcept + { + return std::max (derived ().underlying_scheme.get_max_num_children (), + TSubelementSchemeSpecialization::subelement_get_max_num_children ()); + } + + /** + * Indicates if an element is refinable. Possible reasons for being not refinable could be + * that the element has reached its max level or is a subelement. + * \param [in] elem The element to check. + * \return True if the element is refinable. + */ + bool + element_is_refinable (const t8_element_t *elem) const noexcept + { + T8_ASSERT (element_is_valid (elem)); + if (element_is_subelement (elem)) { + // Subelements are not refinable, as they are discarded for the next adaptation cycle. + return false; + } + return derived ().underlying_scheme.element_get_level (element_to_standalone (elem)) < get_maxlevel (); + } + + /** Construct all children of a given element. Not possible for subelements as they have no children. + * \param [in] elem This must be a valid element, bigger than maxlevel. + * \param [in] length The length of the output array \a c must match the number of children. + * See \ref element_get_num_children. + * \param [in,out] c The storage for these \a length elements must exist and match the element class in + * the children's ordering. On output, all children are valid. + * It is valid to call this function with elem = c[0]. + */ + void + element_get_children (const t8_element_t *elem, const int length, t8_element_t *c[]) const noexcept + { + SC_CHECK_ABORT (!element_is_subelement (elem), "element_get_children: Cannot construct child of a subelement.\n"); + T8_ASSERT (element_is_valid (elem)); + const auto *subelement = as_subelement (elem); + + t8_element_t **standalone_children_ptrs = T8_ALLOC (t8_element_t *, length); + for (int ichild = 0; ichild < length; ++ichild) { + auto *child = as_subelement (c[ichild]); + standalone_children_ptrs[ichild] = subelement_to_standalone (child); + reset_subelement_values (child); + } + derived ().underlying_scheme.element_get_children (subelement_to_standalone (subelement), length, + standalone_children_ptrs); + T8_FREE (standalone_children_ptrs); + } + + /** Compute the child id of an element. + * \param [in] elem This must be a valid element. + * \return The child id of elem. + */ + int + element_get_child_id (const t8_element_t *elem) const noexcept + { + T8_ASSERT (element_is_valid (elem)); + const auto *subelement = as_subelement (elem); + if (element_is_subelement (elem)) { + // For subelements, the child id is the subelement id. + return subelement->subelement_id; + } + return derived ().underlying_scheme.element_get_child_id (subelement_to_standalone (subelement)); + } + + /** Compute the ancestor id of an element, that is the child id at a given level. + * \param [in] elem This must be a valid element. + * \param [in] level A refinement level. Must satisfy \a level < elem.level + * \return The child_id of \a elem in regard to its \a level ancestor. + */ + int + element_get_ancestor_id (const t8_element_t *elem, const t8_element_level level) const noexcept + { + SC_CHECK_ABORT (!element_is_subelement (elem), "element_get_ancestor_id is not implemented for subelements yet.\n"); + return derived ().underlying_scheme.element_get_ancestor_id (element_to_standalone (elem), level); + } + + /** Query whether a given set of elements is a family or not. + * \param [in] fam An array of as many elements as an element of class + * \b ts has siblings. + * \return Zero if \b fam is not a family, nonzero if it is. + * \note level 0 elements do not form a family. + */ + int + elements_are_family (t8_element_t *const *fam) const noexcept + { + const int num_siblings = element_get_num_siblings (fam[0]); +#if T8_ENABLE_DEBUG + for (int isib = 0; isib < num_siblings; isib++) { + T8_ASSERT (element_is_valid (fam[isib])); + } +#endif + /* If the first element is a subelement, the remaining elements also have to be subelements and the standalone + * elements must be equal. */ + if (element_is_subelement (fam[0])) { + auto element_0 = element_to_standalone (fam[0]); + for (int isib = 1; isib < num_siblings; ++isib) { + if (!element_is_subelement (fam[isib]) + || !derived ().underlying_scheme.element_is_equal (element_0, element_to_standalone (fam[isib]))) { + return 0; + } + } + return 1; + } + /* If the first element is no subelement, the remaining elements should not be subelements and + * they must form a family. */ + t8_element_t **standalone_children_ptrs = T8_ALLOC (t8_element_t *, num_siblings); + for (int isib = 0; isib < num_siblings; ++isib) { + if (element_is_subelement (fam[isib])) { + T8_FREE (standalone_children_ptrs); + return 0; + } + standalone_children_ptrs[isib] = element_to_standalone (fam[isib]); + } + + bool are_family = derived ().underlying_scheme.elements_are_family (standalone_children_ptrs); + T8_FREE (standalone_children_ptrs); + return are_family; + } + + /** Query whether element A is an ancestor of the element B. + * An element A is ancestor of an element B if A == B or if B can + * be obtained from A via successive refinement. + * \param [in] element_A An element of class \a eclass in scheme \a scheme. + * \param [in] element_B An element of class \a eclass in scheme \a scheme. + * \return True if and only if \a element_A is an ancestor of \a element_B. + */ + bool + element_is_ancestor (const t8_element_t *element_A, const t8_element_t *element_B) const noexcept + { + T8_ASSERT (element_is_valid (element_A)); + T8_ASSERT (element_is_valid (element_B)); + if (element_is_equal (element_A, element_B)) { + return true; + } + if (element_is_subelement (element_A)) { + // Subelements are not ancestors of any element, as they are discarded for the next adaptation cycle. + // B could be a subelement if the underlying element is an ancestor of A. + return false; + } + return derived ().underlying_scheme.element_is_ancestor (element_to_standalone (element_A), + element_to_standalone (element_B)); + } + + /** Compute the nearest common ancestor of two elements. Not implemented yet. + * \param [in] elem1 The first of the two input elements. + * \param [in] elem2 The second of the two input elements. + * \param [in,out] nca The storage for this element must exist and match the element class of the child. + * On output the unique nearest common ancestor of \b elem1 and \b elem2. + */ + void + element_get_nca ([[maybe_unused]] const t8_element_t *elem1, [[maybe_unused]] const t8_element_t *elem2, + [[maybe_unused]] t8_element_t *nca) const noexcept + { + SC_CHECK_ABORT ((!element_is_subelement (elem1)) && (!element_is_subelement (elem2)), + "element_get_nca is not implemented for subelements yet.\n"); + derived ().underlying_scheme.element_get_nca (element_to_standalone (elem1), element_to_standalone (elem2), + element_to_standalone (nca)); + } + + /** Compute the first descendant of a given element. + * The first descendant of a subelement is the descendant of the parent element, + * as they are discarded for the next adaptation cycle. + * \param [in] elem The element whose descendant is computed. + * \param [out] desc The first element in a uniform refinement of \a elem of the given level. + * \param [in] level The level, at which the descendant is computed. + */ + void + element_get_first_descendant (const t8_element_t *elem, t8_element_t *desc, + const t8_element_level level) const noexcept + { + derived ().underlying_scheme.element_get_first_descendant (element_to_standalone (elem), + element_to_standalone (desc), level); + reset_subelement_values (as_subelement (desc)); + } + + /** Compute the last descendant of a given element. + * The last descendant of a subelement is the descendant of the parent element, + * as they are discarded for the next adaptation cycle. + * \param [in] elem The element whose descendant is computed. + * \param [out] desc The last element in a uniform refinement of \a elem of the given level. + * \param [in] level The level, at which the descendant is computed. + */ + void + element_get_last_descendant (const t8_element_t *elem, t8_element_t *desc, + const t8_element_level level) const noexcept + { + derived ().underlying_scheme.element_get_last_descendant (element_to_standalone (elem), + element_to_standalone (desc), level); + reset_subelement_values ((TSubelementType *) desc); + } + + // ################################################____FACE REFINEMENT____############################################ + + /** Return the number of children of an element's face when the element is refined. + * \note This is not implemented for subelements. + * \param [in] elem The element whose face is considered. + * \param [in] face A face of \a elem. + * \return The number of children of \a face if \a elem is to be refined. + */ + int + element_get_num_face_children ([[maybe_unused]] const t8_element_t *elem, + [[maybe_unused]] const int face) const noexcept + { + SC_CHECK_ABORT (!element_is_subelement (elem), + "element_get_num_face_children is not implemented for subelements yet.\n"); + return derived ().underlying_scheme.element_get_num_face_children (element_to_standalone (elem), face); + } + + /** Given an element and a face of the element, compute all children of the element that touch the face. + * \note This is not implemented for subelements. + * \param [in] elem The element. + * \param [in] face A face of \a elem. + * \param [in,out] children Allocated elements, in which the children of \a elem that share a face with + * \a face are stored. They will be stored in order of their linear id. + * \param [in] num_children The number of elements in \a children. Must match the number of children + * that touch \a face. See \ref element_get_num_face_children. + * \param [in,out] child_indices If not NULL, an array of num_children integers must be given, + * on output its i-th entry is the child_id of the i-th face_child. + * It is valid to call this function with elem = children[0]. + */ + void + element_get_children_at_face ([[maybe_unused]] const t8_element_t *elem, const int face, t8_element_t *children[], + const int num_children, [[maybe_unused]] int *child_indices) const noexcept + { + SC_CHECK_ABORT (!element_is_subelement (elem), + "element_get_children_at_face is not implemented for subelements yet.\n"); + t8_element_t **standalone_children_ptrs = T8_ALLOC (t8_element_t *, num_children); + for (int ichild = 0; ichild < num_children; ++ichild) { + auto *child = as_subelement (children[ichild]); + standalone_children_ptrs[ichild] = subelement_to_standalone (child); + reset_subelement_values (child); + } + derived ().underlying_scheme.element_get_children_at_face (element_to_standalone (elem), face, + standalone_children_ptrs, num_children, child_indices); + T8_FREE (standalone_children_ptrs); + } + + /** Given a face of an element and a child number of a child of that face, return the face number + * of the child of the element that matches the child face. + * \note This is not implemented for subelements. + * \param [in] elem The element. + * \param [in] face Then number of the face. + * \param [in] face_child A number 0 <= \a face_child < num_face_children, specifying a child of \a elem that + * shares a face with \a face. These children are counted in linear order. + * This coincides with the order of children from a call to \ref element_get_children_at_face. + * \return The face number of the face of a child of \a elem that coincides with \a face_child. + */ + int + element_face_get_child_face ([[maybe_unused]] const t8_element_t *elem, const int face, + [[maybe_unused]] const int face_child) const noexcept + { + SC_CHECK_ABORT (!element_is_subelement (elem), + "element_face_get_child_face is not implemented for subelements yet.\n"); + return derived ().underlying_scheme.element_face_get_child_face (element_to_standalone (elem), face, face_child); + } + + /** Given a face of an element return the face number of the parent of the element that matches the element's face. + * \note This is not implemented for subelements. + * \param [in] elem The element. + * \param [in] face Then number of the face. + * \return If \a face of \a elem is also a face of \a elem's parent, the face number of this face. + * Otherwise -1. + * \note For the root element this function always returns \a face. + */ + int + element_face_get_parent_face ([[maybe_unused]] const t8_element_t *elem, const int face) const noexcept + { + SC_CHECK_ABORT (!element_is_subelement (elem), + "element_face_get_parent_face is not implemented for subelements yet.\n"); + return derived ().underlying_scheme.element_face_get_parent_face (element_to_standalone (elem), face); + } + + /** Construct the first descendant of an element at a given level that touches a given face. + * \note This is not implemented for subelements. + * \param [in] elem The input element. + * \param [in] face A face of \a elem. + * \param [in, out] first_desc An allocated element. This element's data will be filled with the data of the first + * descendant of \a elem that shares a face with \a face. + * \param [in] level The level, at which the first descendant is constructed + */ + void + element_get_first_descendant_face (const t8_element_t *elem, const int face, t8_element_t *first_desc, + const t8_element_level level) const noexcept + { + SC_CHECK_ABORT (!element_is_subelement (elem), + "element_get_first_descendant_face is not implemented for subelements yet.\n"); + derived ().underlying_scheme.element_get_first_descendant_face (element_to_standalone (elem), face, + element_to_standalone (first_desc), level); + } + + /** Construct the last descendant of an element at a given level that touches a given face. + * \note This is not implemented for subelements. + * \param [in] elem The input element. + * \param [in] face A face of \a elem. + * \param [in, out] last_desc An allocated element. This element's data will be filled with the data of the last + * descendant of \a elem that shares a face with \a face. + * \param [in] level The level, at which the last descendant is constructed + */ + void + element_get_last_descendant_face ([[maybe_unused]] const t8_element_t *elem, const int face, t8_element_t *last_desc, + const t8_element_level level) const noexcept + { + SC_CHECK_ABORT (!element_is_subelement (elem), + "element_get_last_descendant_face is not implemented for subelements yet.\n"); + derived ().underlying_scheme.element_get_last_descendant_face (element_to_standalone (elem), face, + element_to_standalone (last_desc), level); + } + + // ################################################____FACE NEIGHBOR____############################################## + + /** Compute whether a given element shares a given face with its root tree. + * \note This is not implemented for subelements. + * \param [in] elem The input element. + * \param [in] face A face of \a elem. + * \return True if \a face is a subface of the element's root element. + * \note You can compute the corresponding face number of the tree via \ref element_get_tree_face. + */ + int + element_is_root_boundary (const t8_element_t *elem, const int face) const noexcept + { + SC_CHECK_ABORT (!element_is_subelement (elem), + "element_is_root_boundary is not implemented for subelements yet.\n"); + return derived ().underlying_scheme.element_is_root_boundary (element_to_standalone (elem), face); + } + + /** Given an element and a face of this element. If the face lies on the tree boundary, return the face number + * of the tree face. If not the return value is arbitrary. + * You can call \ref t8_element_is_root_boundary to query whether the face is at the tree boundary. + * \note This is not implemented for subelements. + * \param [in] elem The element. + * \param [in] face The index of a face of \a elem. + * \return The index of the tree face that \a face is a subface of, if \a face is on a tree boundary. + * Any arbitrary integer if \a is not at a tree boundary. + * \warning The return value may look like a valid face of the tree even if the element does not lie on the root boundary. + */ + int + element_get_tree_face (const t8_element_t *elem, const int face) const noexcept + { + SC_CHECK_ABORT (!element_is_subelement (elem), "element_get_tree_face is not implemented for subelements yet.\n"); + return derived ().underlying_scheme.element_get_tree_face (element_to_standalone (elem), face); + } + + /** Construct the face neighbor of a given element if this face neighbor is inside the root tree. Return 0 otherwise. + * \note This is not implemented for subelements. + * \param [in] elem The element to be considered. + * \param [in,out] neigh If the face neighbor of \a elem along \a face is inside the root tree, this element's data + * is filled with the data of the face neighbor. Otherwise the data can be modified arbitrarily. + * \param [in] face The number of the face along which the neighbor should be constructed. + * \param [out] neigh_face The number of \a face as viewed from \a neigh. + * An arbitrary value, if the neighbor is not inside the root tree. + * \return True if \a neigh is inside the root tree. + * False if not. In this case \a neigh's data can be arbitrary on output. + */ + int + element_get_face_neighbor_inside (const t8_element_t *elem, t8_element_t *neigh, const int face, + int *neigh_face) const noexcept + { + SC_CHECK_ABORT (!element_is_subelement (elem), + "element_get_face_neighbor_inside is not implemented for subelements yet.\n"); + return derived ().underlying_scheme.element_get_face_neighbor_inside ( + element_to_standalone (elem), element_to_standalone (neigh), face, neigh_face); + } + + // ######################################____TREE FACE TRANSFORMATION____############################################# + + /** \note This is not implemented for this scheme. + * \param [in] elem1 + * \param [in,out] elem2 + * \param [in] orientation + * \param [in] sign + * \param [in] is_smaller_face + */ + static void + element_transform_face ([[maybe_unused]] const t8_element_t *elem1, [[maybe_unused]] t8_element_t *elem2, + [[maybe_unused]] const int orientation, [[maybe_unused]] const int sign, + [[maybe_unused]] const int is_smaller_face) noexcept + { + SC_ABORT ("Not implemented.\n"); + } + + /** \note This is not implemented for this scheme. + * \param [in] face + * \param [in,out] elem + * \param [in] root_face + * \param [in] scheme + */ + int + element_extrude_face ([[maybe_unused]] const t8_element_t *face, [[maybe_unused]] t8_element_t *elem, + [[maybe_unused]] const int root_face, [[maybe_unused]] const t8_scheme *scheme) const noexcept + { + SC_CHECK_ABORT (!element_is_subelement (elem), "element_extrude_face is not implemented for subelements yet.\n"); + return derived ().underlying_scheme.element_extrude_face (element_to_standalone (face), + element_to_standalone (elem), root_face, scheme); + } + + /** \note This is not implemented for this scheme. + * \param [in] elem + * \param [in] face + * \param [in,out] boundary + * \param [in] scheme + */ + void + element_get_boundary_face ([[maybe_unused]] const t8_element_t *elem, [[maybe_unused]] const int face, + [[maybe_unused]] t8_element_t *boundary, + [[maybe_unused]] const t8_scheme *scheme) const noexcept + { + SC_CHECK_ABORT (!element_is_subelement (elem), + "element_get_boundary_face is not implemented for subelements yet.\n"); + return derived ().underlying_scheme.element_get_boundary_face (element_to_standalone (elem), face, + element_to_standalone (boundary), scheme); + } + + // ################################################____LINEAR ID____################################################ + + /** Initialize the entries of an allocated element according to a given linear id in a uniform refinement. + * \note This is not implemented for subelements. + * \param [in,out] elem The element whose entries will be set. + * \param [in] level The level of the uniform refinement to consider. + * \param [in] id The linear id. id must fulfil 0 <= id < 'number of leaves in the uniform refinement' + */ + void + element_set_linear_id (t8_element_t *elem, const t8_element_level level, t8_linearidx_t id) const noexcept + { + SC_CHECK_ABORT (!element_is_subelement (elem), "element_set_linear_id is not implemented for subelements yet.\n"); + derived ().underlying_scheme.element_set_linear_id (element_to_standalone (elem), level, id); + } + + /** Compute the linear id of a given element in a hypothetical uniform refinement of a given level. + * \note that the id of a subelement equals the id of its parent. + * Therefore, the binary search (for example used in the leaf_face_neighbor function) will find a random subelement + * of the transition cell which might not be the desired neighbor of a given element. + * \param [in] elem The element whose id we compute. + * \param [in] level The level of the uniform refinement to consider. + * \return The linear id of the element. + */ + t8_linearidx_t + element_get_linear_id (const t8_element_t *elem, const t8_element_level level) const noexcept + { + T8_ASSERT (element_is_valid (elem)); + return derived ().underlying_scheme.element_get_linear_id (element_to_standalone (elem), level); + } + + /** Construct the successor in a uniform refinement of a given element. + * \note This is not implemented for subelements. + * \param [in] elem1 The element whose successor should be constructed. + * \param [in,out] elem2 The element whose entries will be set. + */ + void + element_construct_successor (const t8_element_t *elem1, t8_element_t *elem2) const noexcept + { + SC_CHECK_ABORT (!element_is_subelement (elem1), + "element_construct_successor is not implemented for subelements yet.\n"); + return derived ().underlying_scheme.element_construct_successor (element_to_standalone (elem1), + element_to_standalone (elem2)); + } + + /** Count how many leaf descendants of a given uniform level an element would produce. + * \note This is not implemented for subelements. + * \param [in] elem The element to be checked. + * \param [in] level A refinement level. + */ + t8_gloidx_t + element_count_leaves (const t8_element_t *elem, const t8_element_level level) const noexcept + { + SC_CHECK_ABORT (!element_is_subelement (elem), "element_count_leaves is not implemented for subelements yet.\n"); + return derived ().underlying_scheme.element_count_leaves (element_to_standalone (elem), level); + } + + /** Count how many leaf descendants of a given uniform level the root element will produce. + * \param [in] level A refinement level. + * \return The value of \ref t8_element_count_leaves if the input element + * is the root (level 0) element. + */ + t8_gloidx_t + count_leaves_from_root (const t8_element_level level) const noexcept + { + return derived ().underlying_scheme.count_leaves_from_root (level); + } + + /** Compare two elements. + * \note This is not implemented for subelements. + * \param [in] elem1 The first element. + * \param [in] elem2 The second element. + */ + int + element_compare (const t8_element_t *elem1, const t8_element_t *elem2) const noexcept + { + SC_CHECK_ABORT (!element_is_subelement (elem1) && !element_is_subelement (elem2), + "element_compare is not implemented for subelements yet.\n"); + return derived ().underlying_scheme.element_compare (element_to_standalone (elem1), element_to_standalone (elem2)); + } + + // ################################################____VISUALIZATION____############################################## + + /** Compute the coordinates of a given element vertex inside a reference tree + * that is embedded into [0,1]^d (d = dimension). + * \note This is not implemented for subelements. + * \param [in] elem The element to be considered. + * \param [in] vertex The id of the vertex whose coordinates shall be computed. + * \param [out] coords An array of at least as many doubles as the element's dimension + * whose entries will be filled with the coordinates of \a vertex. + */ + void + element_get_vertex_reference_coords (const t8_element_t *elem, const int vertex, double coords[]) const noexcept + { + SC_CHECK_ABORT (!element_is_subelement (elem), + "element_get_vertex_reference_coords is not implemented for subelements yet.\n"); + derived ().underlying_scheme.element_get_vertex_reference_coords (element_to_standalone (elem), vertex, coords); + } + + /** Convert a point in the reference space of an element to a point in the reference space of the tree. + * \param [in] elem The element. + * \param [in] ref_coords The coordinates \f$ [0,1]^\mathrm{dim} \f$ of the point in the reference space + * of the element. + * \param [in] num_coords The number of coordinates to evaluate. + * \param [out] out_coords The coordinates of the point in the reference space of the tree. + */ + void + element_get_reference_coords (const t8_element_t *elem, const double *ref_coords, const size_t num_coords, + double *out_coords) const noexcept + { + if (element_is_subelement (elem)) { + derived ().subelement_get_reference_coords (elem, ref_coords, num_coords, out_coords); + } + else { + derived ().underlying_scheme.element_get_reference_coords (element_to_standalone (elem), ref_coords, num_coords, + out_coords); + } + } + + // ################################################____MEMORY____################################################ + + /** Allocate memory for an array of elements of a given class and initialize them. + * \param [in] length The number of elements to be allocated. + * \param [in,out] elems On input an array of \b length many unallocated element pointers. + * On output all these pointers will point to an allocated and initialized element. + * \note Not every element that is created in t8code will be created by a call to this function. However, if an + * element is not created using \ref element_new, then it is guaranteed that \ref element_init is called on it. + * \note In debugging mode, an element that was created with \ref element_new must pass \ref element_is_valid. + * \note If an element was created by \ref element_new then \ref element_init may not be called for it. + * Thus, \ref element_new should initialize an element in the same way as a call to \ref element_init would. + * \see element_init + * \see element_is_valid + */ + void + element_new (const int length, t8_element_t **elems) const noexcept + { + T8_ASSERT (this->scheme_context != NULL); + T8_ASSERT (0 <= length); + T8_ASSERT (elems != NULL); + for (int i = 0; i < length; ++i) { + elems[i] = (t8_element_t *) sc_mempool_alloc ((sc_mempool_t *) this->scheme_context); + } + /* For other schemes, we only set sensible data in debug mode. For subelements, it is important that we always set + * the subelement id and type to zero for new elements. */ + for (int i = 0; i < length; i++) { + element_init (1, elems[i]); + } + } + + /** Initialize an array of allocated elements. + * \param [in] length The number of elements to be initialized. + * \param [in,out] elems On input an array of \b length many allocated elements. + * \note In debugging mode, an element that was passed to \ref element_init must pass \ref element_is_valid. + * \note If an element was created by \ref element_new then \ref element_init may not be called for it. + * Thus, \ref element_new should initialize an element in the same way as a call to \ref element_init would. + * \see element_new + * \see element_is_valid + */ + void + element_init ([[maybe_unused]] const int length, [[maybe_unused]] t8_element_t *elems) const noexcept + { + TSubelementType *subelement = (TSubelementType *) elems; + for (int ielem = 0; ielem < length; ielem++) { + reset_subelement_values (subelement + ielem); + derived ().underlying_scheme.element_init (1, subelement_to_standalone (subelement + ielem)); + T8_ASSERT (element_is_valid ((t8_element_t *) (subelement + ielem))); + } + } + + /** Deinitialize an array of allocated elements. + * \param [in] length The number of elements to be deinitialized. + * \param [in,out] elems On input an array of \a length many allocated and initialized elements, on output an array of + * \a length many allocated, but not initialized elements. + * \note Call this function if you called element_init on the element pointers. + * \see element_init + */ + static constexpr void + element_deinit ([[maybe_unused]] const int length, [[maybe_unused]] t8_element_t *elems) noexcept + { + } + + /** Deallocate an array of elements. + * \param [in] length The number of elements in the array. + * \param [in,out] elems On input an array of \b length many allocated element pointers. + * On output all these pointers will be freed. \b elems itself will not be freed by this function. + */ + void + element_destroy (const int length, t8_element_t **elems) const noexcept + { + T8_ASSERT (this->scheme_context != NULL); + T8_ASSERT (0 <= length); + T8_ASSERT (elems != NULL); + for (int i = 0; i < length; ++i) { + sc_mempool_free ((sc_mempool_t *) scheme_context, elems[i]); + } + } + + // ############################################____DEBUG____########################################################## + +#if T8_ENABLE_DEBUG + /** Query whether a given element can be considered as 'valid' and it is safe to perform any of the above + * algorithms on it. For example this could mean that all coordinates are in valid ranges and other member variables + * do have meaningful values. + * \param [in] elem The element to be checked. + * \return True if \a elem is safe to use. False otherwise. + * \note This function is used for debugging to catch certain errors. These can for example occur when + * an element points to a region of memory which should not be interpreted as an element. + * \note We recommend to use the assertion T8_ASSERT (element_is_valid (elem)) + * in the implementation of each of the functions in this file. + */ + int + element_is_valid (const t8_element_t *elem) const noexcept + { + T8_ASSERT (elem != NULL); + const auto *subelement = as_subelement (elem); + int element_valid = derived ().underlying_scheme.element_is_valid (subelement_to_standalone (subelement)); + if (!element_is_subelement (elem)) { + return element_valid; + } + // Subelement type 0 always means no subelement. + bool subelement_valid + = (subelement->subelement_type >= 1 + && subelement->subelement_type <= TSubelementSchemeSpecialization::subelement_get_number_of_valid_types ()) + && (subelement->subelement_id >= 0 + && subelement->subelement_id + < TSubelementSchemeSpecialization::element_get_number_of_subelements (subelement->subelement_type)); + + return subelement_valid && element_valid; + } + + /** Print a given element. This prints the subelement information and the standalone debug print. + * This function is only available in the debugging configuration. + * \param [in] elem The element to print + */ + void + element_debug_print (const t8_element_t *elem) const noexcept + { + const auto *subelement = as_subelement (elem); + t8_debugf ("Subelement type: %i\n", subelement->subelement_type); + t8_debugf ("Subelement id: %i\n", subelement->subelement_id); + derived ().underlying_scheme.element_debug_print (subelement_to_standalone (subelement)); + } + +#endif + /** Fill a string with readable information about the element + * \param[in] elem The element to translate into human-readable information + * \param[in, out] debug_string The string to fill. + * \param[in] string_size Buffer size of c-string + */ + static void + element_to_string ([[maybe_unused]] const t8_element_t *elem, [[maybe_unused]] char *debug_string, + [[maybe_unused]] const int string_size) noexcept + { + SC_ABORT ("Not implemented."); + } + + // ################################################____MPI____######################################################## + /** Pack multiple elements into contiguous memory, so they can be sent via MPI. + * \param [in] elements Array of elements that are to be packed + * \param [in] count Number of elements to pack + * \param [in,out] send_buffer Buffer in which to pack the elements + * \param [in] buffer_size size of the buffer (in order to check that we don't access out of range) + * \param [in, out] position the position of the first byte that is not already packed + * \param [in] comm MPI Communicator + */ + void + element_MPI_Pack (t8_element_t **const elements, const unsigned int count, void *send_buffer, const int buffer_size, + int *position, sc_MPI_Comm comm) const noexcept + + { + TSubelementType **els = (TSubelementType **) elements; + for (unsigned int ielem = 0; ielem < count; ielem++) { + t8_element_t *element = subelement_to_standalone (els[ielem]); + derived ().underlying_scheme.element_MPI_Pack (&element, 1, send_buffer, buffer_size, position, comm); + int mpiret = sc_MPI_Pack (&els[ielem]->subelement_type, 1, sc_MPI_INT, send_buffer, buffer_size, position, comm); + SC_CHECK_MPI (mpiret); + mpiret = sc_MPI_Pack (&els[ielem]->subelement_id, 1, sc_MPI_INT, send_buffer, buffer_size, position, comm); + SC_CHECK_MPI (mpiret); + } + } + + /** Determine an upper bound for the size of the packed message of \a count elements. + * \param [in] count Number of elements to pack + * \param [in] comm MPI Communicator + * \param [out] pack_size upper bound on the message size + */ + void + element_MPI_Pack_size (const unsigned int count, sc_MPI_Comm comm, int *pack_size) const noexcept + { + // Get single size from standalone scheme. + derived ().underlying_scheme.element_MPI_Pack_size (1, comm, pack_size); + int singlesize = *pack_size; + + /* Type and id are both of type int. */ + int datasize = 0; + int mpiret = sc_MPI_Pack_size (1, sc_MPI_INT, comm, &datasize); + SC_CHECK_MPI (mpiret); + singlesize += 2 * datasize; + + *pack_size = count * singlesize; + } + + /** Unpack multiple elements from contiguous memory that was received via MPI. + * \param [in] recvbuf Buffer from which to unpack the elements + * \param [in] buffer_size size of the buffer (in order to check that we don't access out of range) + * \param [in, out] position the position of the first byte that is not already packed + * \param [in] elements Array of initialised elements that is to be filled from the message + * \param [in] count Number of elements to unpack + * \param [in] comm MPI Communicator + */ + void + element_MPI_Unpack (void *recvbuf, const int buffer_size, int *position, t8_element_t **elements, + const unsigned int count, sc_MPI_Comm comm) const noexcept + { + TSubelementType **els = (TSubelementType **) elements; + for (unsigned int ielem = 0; ielem < count; ielem++) { + t8_element_t *single = subelement_to_standalone (els[ielem]); + derived ().underlying_scheme.element_MPI_Unpack (recvbuf, buffer_size, position, &single, 1, comm); + int mpiret = sc_MPI_Unpack (recvbuf, buffer_size, position, &els[ielem]->subelement_type, 1, sc_MPI_INT, comm); + SC_CHECK_MPI (mpiret); + mpiret = sc_MPI_Unpack (recvbuf, buffer_size, position, &els[ielem]->subelement_id, 1, sc_MPI_INT, comm); + SC_CHECK_MPI (mpiret); + } + } + + // ########################################____SUBELEMENTS____######################################################## + /** Check if \a elem is a subelement. + * \param [in] elem The elem to be checked. + */ + static bool + element_is_subelement (const t8_element_t *elem) + { + const auto *subelement = as_subelement (elem); + return (subelement->subelement_type != 0); + } + + /** Get the number of subelements an element is refined into for a specific type. + * \param [in] subelement_type The subelement type used for refinement. + */ + static int + element_get_number_of_subelements (int subelement_type) + { + return TSubelementSchemeSpecialization::element_get_number_of_subelements (subelement_type); + } + + /** This defines how an element is refined in subelements using a specified subelement type. + * \param [in] elem The element to be refined. + * \param [in] type The subelement type to be used for refinement. + * \param [in, out] c An array of allocated elements that will be filled with the subelements of \a elem. + * The number of subelements is determined by \ref element_get_number_of_subelements. + */ + void + refine_element_in_subelements (const t8_element_t *elem, int type, t8_element_t *c[]) const noexcept + { + derived ().refine_element_in_subelements (elem, type, c); + } + + protected: + // PRIVATE HELPER + /** Return the standalone element stored inside a subelement. + * Const version. + * \param[in] subelement The subelement for which the standalone element should be extracted. + */ + static const t8_element_t * + subelement_to_standalone (const TSubelementType *subelement) noexcept + { + return (const t8_element_t *) &subelement->element; + } + + /** Return the standalone element stored inside a subelement. + * Non-const version. + * \param[in] subelement The subelement for which the standalone element should be extracted. + */ + static t8_element_t * + subelement_to_standalone (TSubelementType *subelement) noexcept + { + return (t8_element_t *) &subelement->element; + } + + /** Same as \ref subelement_to_standalone but takes a general element type. + * Const version. + * \param[in] element Element of class t8_element_t that can be interpreted as subelement. + */ + static const t8_element_t * + element_to_standalone (const t8_element_t *element) noexcept + { + return subelement_to_standalone (as_subelement (element)); + } + + /** Same as \ref subelement_to_standalone but takes a general element type. + * Non-const version. + * \param[in] element Element of class t8_element_t that can be interpreted as subelement. + */ + static t8_element_t * + element_to_standalone (t8_element_t *element) noexcept + { + return subelement_to_standalone (as_subelement (element)); + } + + /** Interpret an element as a subelement. + * Const version. + * \param[in] element Element of class t8_element_t that can be interpreted as subelement. + */ + static const TSubelementType * + as_subelement (const t8_element_t *element) noexcept + { + return reinterpret_cast (element); + } + /** Interpret an element as a subelement. + * Non-const version. + * \param[in] element Element of class t8_element_t that can be interpreted as subelement. + */ + static TSubelementType * + as_subelement (t8_element_t *element) noexcept + { + return reinterpret_cast (element); + } + + /** Reset the subelement-specific data. + * \param [in,out] subelement The element that is filled with the root of the subelement. + */ + static void + reset_subelement_values (TSubelementType *subelement) noexcept + { + subelement->subelement_type = 0; + subelement->subelement_id = 0; + } + + /** Return the edge length of the parent standalone element. + * \param [in] subelement A subelement of the parent element. + * \return The length of the parent element in integer coordinates. + */ + t8_element_coord + parent_element_get_len (const TSubelementType *subelement) const noexcept + { + return 1 << (T8_ELEMENT_MAXLEVEL[TEclass] + - (derived ().underlying_scheme.element_get_level (subelement_to_standalone (subelement)))); + } + + /** Return the derived subelement scheme. + * Const Version. + */ + const TSubelementSchemeSpecialization & + derived () const noexcept + { + return static_cast (*this); + } + + /** Return the derived subelement scheme. + * Non-const Version. + */ + TSubelementSchemeSpecialization & + derived () noexcept + { + return static_cast (*this); + } +}; diff --git a/src/t8_schemes/t8_subelement/t8_subelement_traits.hxx b/src/t8_schemes/t8_subelement/t8_subelement_traits.hxx new file mode 100644 index 0000000000..4c80ea2f22 --- /dev/null +++ b/src/t8_schemes/t8_subelement/t8_subelement_traits.hxx @@ -0,0 +1,68 @@ +/* + 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_subelement_traits.hxx + * Traits for subelement schemes. + * Each trait defines the underlying scheme and subelement type associated with the subelement scheme. + * This allows the common subelement implementation in \ref t8_subelement_scheme.hxx to remain generic + * while using the types required by each specialization. + */ + +#pragma once + +#include +#include +#include +#include +#include + +/** Forward declaration of the quadrilateral subelement scheme. */ +struct t8_subelementquad_scheme; + +/** Forward declaration of the triangular subelement scheme. */ +struct t8_subelementtri_scheme; + +/** Traits associating a subelement scheme with its underlying scheme and subelement type. + * \tparam TScheme The subelement scheme for which the traits are defined. + */ +template +struct t8_subelement_traits; + +/** Traits specialization for quadrilateral subelements. */ +template <> +struct t8_subelement_traits +{ + /** Subelement class used for the quad scheme. See also \ref t8_subelement_type.hxx. */ + using SubelementType = t8_subelement_element>; + /** The standalone scheme for recursive refinement that gets extended by the hanging node resolution for quads. */ + using UnderlyingScheme = t8_standalone_scheme; +}; + +/** Traits specialization for triangular subelements. */ +template <> +struct t8_subelement_traits +{ + /** Subelement class used for the triangle scheme. See also \ref t8_subelement_type.hxx. */ + using SubelementType = t8_subelement_element; + /** The scheme for recursive refinement that gets extended by the hanging node resolution for triangles. */ + using UnderlyingScheme = t8_default_scheme_tri; +}; diff --git a/src/t8_schemes/t8_subelement/t8_subelement_type.hxx b/src/t8_schemes/t8_subelement/t8_subelement_type.hxx new file mode 100644 index 0000000000..0b6b951f35 --- /dev/null +++ b/src/t8_schemes/t8_subelement/t8_subelement_type.hxx @@ -0,0 +1,45 @@ +/* + 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_subelement_type.hxx + * Definition of the element class of a subelement. A subelement always consists of an underlying element and + * subelement type and id defining how the underlying element is transitioned into a subelement. + */ + +#pragma once + +/** Definition of the subelement class. A subelement always has an underlying element. + * With the type, it is defined if the element is further defined into subelements (e.g. for hanging node resolution). + * Type 0 means no subelement and the subelement is just the underlying element. + * For hanging node resolution, the type encodes which faces are hanging and therefore the number of subelements in + * which the element is transitioned. + * Accordingly, the subelement id is between 0 and num_subelement - 1. + * \tparam TUnderlyingElement The type of the underlying element. For example a standalone element. + */ +template +struct t8_subelement_element +{ + TUnderlyingElement element; /**< Standalone element of the subelement. */ + int subelement_type + = 0; /**< Type of the transition cell a subelement is associated to (default is 0, meaning no subelement). */ + int subelement_id = 0; /**< Id of the children subelement the given element is (default is 0). */ +}; diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index f1e2f85152..0c41dd2bc1 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -160,6 +160,7 @@ add_t8_cpp_test( NAME t8_gtest_partition_data_parallel SOURCES t8_for add_t8_cpp_test( NAME t8_gtest_set_partition_offset_parallel SOURCES t8_forest/t8_gtest_set_partition_offset.cxx ) add_t8_cpp_test( NAME t8_gtest_partition_for_coarsening_parallel SOURCES t8_forest/t8_gtest_partition_for_coarsening.cxx ) add_t8_cpp_test( NAME t8_gtest_weighted_partitioning_parallel SOURCES t8_forest/t8_gtest_weighted_partitioning.cxx ) +add_t8_cpp_test( NAME t8_gtest_subelement_parallel SOURCES t8_forest/t8_gtest_subelement.cxx t8_gtest_adapt_callbacks.cxx ) add_t8_cpp_test( NAME t8_gtest_element_is_leaf_parallel SOURCES t8_forest/t8_gtest_element_is_leaf.cxx t8_gtest_adapt_callbacks.cxx ) add_t8_cpp_test( NAME t8_gtest_permute_hole_serial SOURCES t8_forest_incomplete/t8_gtest_permute_hole.cxx ) diff --git a/test/t8_forest/t8_gtest_subelement.cxx b/test/t8_forest/t8_gtest_subelement.cxx new file mode 100644 index 0000000000..275b563b12 --- /dev/null +++ b/test/t8_forest/t8_gtest_subelement.cxx @@ -0,0 +1,91 @@ +/* + 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_subelement.cxx + * Unit test for the subelement scheme. This test checks the complete pipeline of hanging node resolution for + *´hybrid 2D meshes. Currently, we only test that the functions required for visualization work correctly. + */ + +#include +#include + +#include +#include +#include +#include +#include +#include + +/** Check that the hanging node resolution for 2D hybrid meshes works. At the moment we only check the functionality +* needed for visualization (so e.g. no connectivity). +*/ +TEST (t8_gtest_subelement, hybrid_hanging_nodes_visualization) +{ + /* Setup: Build hypercube cmesh and uniform forest with the subelement scheme. */ + const int level = 3; + t8_cmesh_t cmesh; + t8_cmesh_init (&cmesh); + t8_cmesh_new_2D_hypercube_hybrid (cmesh, sc_MPI_COMM_WORLD); + + t8_forest_t forest = t8_forest_new_uniform (cmesh, t8_scheme_new_subelement (), level, 0, sc_MPI_COMM_WORLD); + + /* Initial uniform forest should not have any subelements. */ + EXPECT_FALSE (t8_forest_has_global_subelements (forest)); + + /* Adapt the forest (refining every second element). */ + forest = t8_forest_new_adapt (forest, t8_test_adapt_even_global_id, 0, 0, NULL); + + /* Before resolving hanging nodes, subelements should not yet be introduced. */ + EXPECT_FALSE (t8_forest_has_global_subelements (forest)); + const t8_gloidx_t num_leaves_adapted = t8_forest_get_global_num_leaf_elements (forest); + + /* Remove hanging nodes by inserting subelements. The forest is already balanced as we only adapted once. */ + forest = t8_forest_remove_hanging_nodes (forest); + EXPECT_TRUE (t8_forest_is_committed (forest)); + + /* Hanging node resolution must introduce subelements into the forest. */ + EXPECT_TRUE (t8_forest_has_global_subelements (forest)); + + /* Adding transition subelements must increase (or equal) the total leaf count. */ + const t8_gloidx_t num_leaves_sub = t8_forest_get_global_num_leaf_elements (forest); + EXPECT_GT (num_leaves_sub, num_leaves_adapted); + + /* Repartition the forest containing subelements (exercises MPI_Pack / MPI_Unpack). */ + t8_forest_t forest_partitioned; + t8_forest_init (&forest_partitioned); + t8_forest_set_partition (forest_partitioned, forest, 0); + t8_forest_commit (forest_partitioned); + + /* Subelements and leaf count must remain consistent after repartitioning. */ + EXPECT_TRUE (t8_forest_has_global_subelements (forest_partitioned)); + EXPECT_EQ (t8_forest_get_global_num_leaf_elements (forest_partitioned), num_leaves_sub); + + /* Discard subelements from the partitioned forest. */ + forest = t8_forest_discard_subelements (forest_partitioned); + /* Subelements should now be completely removed. */ + EXPECT_FALSE (t8_forest_has_global_subelements (forest)); + /* Discarding subelements should restore the pre-resolution leaf count. */ + EXPECT_EQ (t8_forest_get_global_num_leaf_elements (forest), num_leaves_adapted); + + /* Clean up forest */ + t8_forest_unref (&forest); +} diff --git a/test/t8_gtest_adapt_callbacks.cxx b/test/t8_gtest_adapt_callbacks.cxx index bab9f26a08..b574f52418 100644 --- a/test/t8_gtest_adapt_callbacks.cxx +++ b/test/t8_gtest_adapt_callbacks.cxx @@ -56,3 +56,15 @@ t8_test_adapt_first_child (t8_forest_t forest, [[maybe_unused]] t8_forest_t fore } return 0; } + +int +t8_test_adapt_even_global_id ([[maybe_unused]] t8_forest_t forest, t8_forest_t forest_from, t8_locidx_t which_tree, + [[maybe_unused]] t8_eclass_t eclass, t8_locidx_t lelement_id, + [[maybe_unused]] const t8_scheme *scheme, [[maybe_unused]] const int is_family, + [[maybe_unused]] const int num_elements, [[maybe_unused]] t8_element_t *elements[]) +{ + if ((t8_forest_get_tree_element_offset (forest_from, which_tree) + lelement_id) % 2 == 0) { + return 1; + } + return 0; +} diff --git a/test/t8_gtest_adapt_callbacks.hxx b/test/t8_gtest_adapt_callbacks.hxx index 6198d0ab5e..9409c0bc34 100644 --- a/test/t8_gtest_adapt_callbacks.hxx +++ b/test/t8_gtest_adapt_callbacks.hxx @@ -40,7 +40,7 @@ * \param [in] forest The forest to which the new elements belong. * \param [in] forest_from The forest that is adapted. * \param [in] which_tree The local tree containing \a elements. - * \param [in] eclass The eclass of \a which_tree. + * \param [in] eclass The eclass of \a which_tree. * \param [in] lelement_id The local element id in \a forest_from in the tree of the current element. * \param [in] scheme The scheme of the forest. * \param [in] is_family If 1, the first \a num_elements entries in \a elements form a family. If 0, they do not. @@ -53,4 +53,22 @@ t8_test_adapt_first_child (t8_forest_t forest, t8_forest_t forest_from, t8_locid const t8_eclass_t eclass, t8_locidx_t lelement_id, const t8_scheme *scheme, const int is_family, const int num_elements, t8_element_t *elements[]); +/** Adapt callback for a forest to refine every second element, so every element with an even global id. + * It is not intended to be used as a recursive adaption callback and does not check the level of an element. + * + * \param [in] forest The forest to which the new elements belong. + * \param [in] forest_from The forest that is adapted. + * \param [in] which_tree The local tree containing \a elements. + * \param [in] eclass The eclass of \a which_tree. + * \param [in] lelement_id The local element id in \a forest_from in the tree of the current element. + * \param [in] scheme The scheme of the forest. + * \param [in] is_family If 1, the first \a num_elements entries in \a elements form a family. If 0, they do not. + * \param [in] num_elements The number of entries in \a elements that are defined + * \param [in] elements Pointers to a family or, if \a is_family is zero, pointer to one element. + */ +int +t8_test_adapt_even_global_id (t8_forest_t forest, t8_forest_t forest_from, t8_locidx_t which_tree, t8_eclass_t eclass, + t8_locidx_t lelement_id, const t8_scheme *scheme, const int is_family, + const int num_elements, t8_element_t *elements[]); + #endif /* T8_GTEST_ADAPT_CALLBACKS */