diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 6d5197d07e..a130355ef0 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -154,6 +154,8 @@ target_sources( T8 PRIVATE t8_cmesh/t8_cmesh_vertex_connectivity/t8_cmesh_vertex_conn_tree_to_vertex.cxx t8_cmesh/t8_cmesh_vertex_connectivity/t8_cmesh_vertex_conn_vertex_to_tree.cxx t8_cmesh/t8_cmesh_vertex_connectivity/t8_cmesh_vertex_connectivity.cxx + t8_cmesh/t8_cmesh_healpix/t8_cmesh_healpix.cxx + t8_cmesh/t8_cmesh_healpix/t8_geometry_healpix.cxx t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx t8_cmesh/t8_cmesh_io/deprecated/deprecated_t8_cmesh_save.cxx t8_cmesh/t8_cmesh_io/deprecated/deprecated_t8_cmesh_triangle.cxx diff --git a/src/t8_cmesh/t8_cmesh_healpix/t8_cmesh_healpix.cxx b/src/t8_cmesh/t8_cmesh_healpix/t8_cmesh_healpix.cxx new file mode 100644 index 0000000000..d9d4e621c5 --- /dev/null +++ b/src/t8_cmesh/t8_cmesh_healpix/t8_cmesh_healpix.cxx @@ -0,0 +1,83 @@ +/* + 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. +*/ + +#include +#include +#include +#include +#include +#include + +t8_cmesh_t +t8_cmesh_new_healpix (sc_MPI_Comm comm) +{ + /* Initialization of the mesh */ + t8_cmesh_t cmesh; + t8_cmesh_init (&cmesh); + + const int ntrees = 12; + const int nverts = 4; /* Number of vertices per cmesh element. */ + t8_eclass_t all_eclasses[ntrees]; + std::vector all_verts; + all_verts.reserve (ntrees * nverts * 3); + + /* Register geometry and retain the pointer */ + t8_geometry_c *geom = t8_cmesh_register_geometry (cmesh); + + /* Reference coordinates for the 4 corners of a quad element */ + double ref_corners[4][2] = { { 0.0, 0.0 }, { 1.0, 0.0 }, { 1.0, 1.0 }, { 0.0, 1.0 } }; + + /* Build trees for all 3 layers (upper, middle, lower) */ + for (int itree = 0; itree < ntrees; itree++) { + const t8_gloidx_t layer = itree / 4; + const t8_gloidx_t face = itree % 4; + + t8_cmesh_set_tree_class (cmesh, itree, T8_ECLASS_QUAD); + all_eclasses[itree] = T8_ECLASS_QUAD; + + /* Associate the custom HEALPix geometry with this tree */ + t8_cmesh_set_tree_geometry (cmesh, itree, geom); + + std::vector verts; + verts.reserve (nverts * 3); + + for (int i = 0; i < nverts; i++) { + double coord[3]; + const double xi = std::clamp (ref_corners[i][0], 1e-10, 1.0 - 1e-10); + const double eta = std::clamp (ref_corners[i][1], 1e-10, 1.0 - 1e-10); + + t8_eval_geom_point (layer, face, xi, eta, coord); + verts.push_back (coord[0]); + verts.push_back (coord[1]); + verts.push_back (coord[2]); + } + + all_verts.insert (all_verts.end (), verts.begin (), verts.end ()); + t8_cmesh_set_tree_vertices (cmesh, itree, verts.data (), nverts); + } + + /* Compute face connectivity using topological vertices */ + t8_cmesh_set_join_by_vertices (cmesh, 12, all_eclasses, all_verts.data (), nullptr, 0); + + t8_cmesh_commit (cmesh, comm); + return cmesh; +} diff --git a/src/t8_cmesh/t8_cmesh_healpix/t8_cmesh_healpix.h b/src/t8_cmesh/t8_cmesh_healpix/t8_cmesh_healpix.h new file mode 100644 index 0000000000..55c35cd9b9 --- /dev/null +++ b/src/t8_cmesh/t8_cmesh_healpix/t8_cmesh_healpix.h @@ -0,0 +1,30 @@ +/* + 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. +*/ + +#include +/** Creates a t8code cmesh consisting of 12 quadrilateral trees arranged + * according to the HEALPix base decomposition + * \param [in] comm MPI Communicator to use. + * \return A coarse mesh representing the HEALPix base grid. + */ +t8_cmesh_t +t8_cmesh_new_healpix (sc_MPI_Comm comm); diff --git a/src/t8_cmesh/t8_cmesh_healpix/t8_geometry_healpix.cxx b/src/t8_cmesh/t8_cmesh_healpix/t8_geometry_healpix.cxx new file mode 100644 index 0000000000..2aba93ed08 --- /dev/null +++ b/src/t8_cmesh/t8_cmesh_healpix/t8_geometry_healpix.cxx @@ -0,0 +1,81 @@ +/* + 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. +*/ + +#include +#include +#include +#include +void +t8_eval_geom_point (t8_gloidx_t layer, t8_gloidx_t face, double xi, double eta, double *out_coord) +{ + const double offset_y = static_cast (layer) - 1.0; + const double offset_x = 2.0 * face + ((layer + 1) % 2); + + const double dX = offset_x + (xi - eta); + const double dY = offset_y + (xi + eta - 1.0); + + double z = 0.0; + double phi = 0.0; + + if (std::abs (dY) <= 1.0) { + z = (2.0 / 3.0) * dY; + phi = (std::numbers::pi / 4.0) * dX; + } + else if (dY > 1.0) { + z = 1.0 - (1.0 / 3.0) * std::pow (2.0 - dY, 2.0); + phi = (std::numbers::pi / 4.0) * (offset_x + (dX - offset_x) / (2.0 - dY)); + } + else { + z = -1.0 + (1.0 / 3.0) * std::pow (2.0 + dY, 2.0); + phi = (std::numbers::pi / 4.0) * (offset_x + (dX - offset_x) / (2.0 + dY)); + } + + z = std::clamp (z, -1.0, 1.0); + + // Convert Spherical coordinates (z, phi) to 3D Cartesian coordinates (x, y, z) + const double theta = std::acos (z); + const double sin_theta = std::sin (theta); + + out_coord[0] = sin_theta * std::cos (phi); + out_coord[1] = sin_theta * std::sin (phi); + out_coord[2] = z; +} + +void +t8_geometry_healpix::t8_geom_evaluate ([[maybe_unused]] t8_cmesh_t cmesh, t8_gloidx_t gtreeid, const double *ref_coords, + const size_t num_coords, double *out_coords) const +{ + // This geometry is inspired by HEALPix + const t8_gloidx_t layer = gtreeid / 4; + const t8_gloidx_t face = gtreeid % 4; + + for (size_t i = 0; i < num_coords; ++i) { + const size_t base_idx_2d = i * 2; + const size_t base_idx_3d = i * 3; // t8code outputs 3D Cartesian coordinates (x, y, z) + + const double xi = std::clamp (ref_coords[base_idx_2d + 0], 1e-10, 1.0 - 1e-10); + const double eta = std::clamp (ref_coords[base_idx_2d + 1], 1e-10, 1.0 - 1e-10); + + // Pass the pointer offset for 3D coordinates: out_coords + base_idx_3d + t8_eval_geom_point (layer, face, xi, eta, out_coords + base_idx_3d); + } +} diff --git a/src/t8_cmesh/t8_cmesh_healpix/t8_geometry_healpix.hxx b/src/t8_cmesh/t8_cmesh_healpix/t8_geometry_healpix.hxx new file mode 100644 index 0000000000..deb04077eb --- /dev/null +++ b/src/t8_cmesh/t8_cmesh_healpix/t8_geometry_healpix.hxx @@ -0,0 +1,96 @@ +/* + 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. +*/ + +#include +#include +/** +* The t8_geometry_healpix struct implements a spherical geometry mapping +* based on the HEALPix (Hierarchical Equal Area isoLatitude Pixelization) +* base-mesh configuration. It inherits from t8_geometry_with_vertices and +* is designed to map planar quadrilateral faces onto a spherical surface. +*/ +struct t8_geometry_healpix: public t8_geometry_with_vertices +{ + public: + /* Basic constructor that sets the dimension and the name. */ + t8_geometry_healpix (): t8_geometry_with_vertices ("t8_geometry_healpix") + { + } + + /* The destructor. */ + virtual ~t8_geometry_healpix () + { + } + + /** + * Map the faces of a cube to a spherical surface. + * \param [in] cmesh The cmesh in which the point lies. + * \param [in] gtreeid The global tree (of the cmesh) in which the reference point is. + * \param [in] ref_coords Array of tree dimension x \a num_coords many entries, specifying a point in \f$ [0,1]^\mathrm{dim} \f$. + * \param [in] num_coords The number of points to map. + * \param [out] out_coords The mapped coordinates in physical space of \a ref_coords. The length is \a num_coords * 3. + * + * This routine expects an input mesh of twelve quadrilaterals in a configuration of the healpix base-mesh + * + */ + void + t8_geom_evaluate (t8_cmesh_t cmesh, t8_gloidx_t gtreeid, const double *ref_coords, const size_t num_coords, + double *out_coords) const; + + /** + * Jacobian, not implemented. + * \param[in] cmesh The cmesh in which the point lies. + * \param[in] gtreeid The global tree (of the cmesh) in which the reference point is. + * \param[in] ref_coords Array of tree dimension x \a num_coords many entries, specifying a point in \f$ [0,1]^\mathrm{dim} \f$. + * \param[in] num_coords The number of points to map. + * \param[in] jacobian The Jacobian matrix to be filled. + */ + void + t8_geom_evaluate_jacobian ([[maybe_unused]] t8_cmesh_t cmesh, [[maybe_unused]] t8_gloidx_t gtreeid, + [[maybe_unused]] const double *ref_coords, [[maybe_unused]] const size_t num_coords, + [[maybe_unused]] double *jacobian) const + { + SC_ABORT_NOT_REACHED (); + } + + /** + * Check for compatibility of the currently loaded tree with the geometry. + * Only hex elements are supported by this geometry. + * \return True if the geometry is compatible with the tree. + */ + bool + t8_geom_check_tree_compatibility () const + { + if (active_tree_class != T8_ECLASS_QUAD) { + t8_productionf ("t8_geometry_healpix is not compatible with tree type %s\n" + "It is only compatible with hex elements.\n", + t8_eclass_to_string[active_tree_class]); + return false; + } + return true; + } + + /* Load tree data is inherited from t8_geometry_with_vertices. */ +}; + +void +t8_eval_geom_point (t8_gloidx_t layer, t8_gloidx_t face, double xi, double eta, double *out_coord); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index f1e2f85152..c24438bc21 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -138,6 +138,7 @@ add_t8_cpp_test( NAME t8_gtest_cmesh_vertex_conn_serial SOUR add_t8_cpp_test( NAME t8_gtest_multiple_attributes_parallel SOURCES t8_cmesh/t8_gtest_multiple_attributes.cxx ) add_t8_cpp_test( NAME t8_gtest_attribute_gloidx_array_serial SOURCES t8_cmesh/t8_gtest_attribute_gloidx_array.cxx ) add_t8_cpp_test( NAME t8_gtest_cmesh_bounding_box_serial SOURCES t8_cmesh/t8_gtest_cmesh_bounding_box.cxx ) +add_t8_cpp_test( NAME t8_gtest_cmesh_healpix SOURCES t8_cmesh/t8_gtest_cmesh_healpix.cxx ) add_t8_cpp_test( NAME t8_gtest_shmem_parallel SOURCES t8_data/t8_gtest_shmem.cxx ) add_t8_cpp_test( NAME t8_gtest_data_pack_parallel SOURCES t8_data/t8_gtest_data_handler.cxx t8_data/t8_data_handler_specs.cxx) diff --git a/test/t8_cmesh/t8_gtest_cmesh_healpix.cxx b/test/t8_cmesh/t8_gtest_cmesh_healpix.cxx new file mode 100644 index 0000000000..60bd8d369b --- /dev/null +++ b/test/t8_cmesh/t8_gtest_cmesh_healpix.cxx @@ -0,0 +1,79 @@ +/* + 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. +*/ + +/* This gtest checks whether a cmesh of a healpix is committed and whether the points have the same radius*/ + +#include +#include +#include +#include +#include +#include +#include +/* Test whether the healpix geometry has an equal radius all across the geometry. */ +TEST (T8GeometryHealpixTest, AllTreesAreOnUnitSphere) +{ + t8_geometry_healpix geom; + t8_cmesh_t cmesh; + t8_cmesh_init (&cmesh); + + std::random_device rd; + std::mt19937 gen (rd ()); + std::uniform_real_distribution dist (0.0, 1.0); + + // A set of 2D reference coordinates within the element + std::vector ref_coords = { + 0.0, 0.0, // Corners + 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.5, 0.5, // Center + }; + int num_gen_coords = 0; +// in addition to the basic test cases generate random coordinates to test out implementation +#if T8_TEST_LEVEL_INT >= 2 + num_gen_coords = 10; // Test level basic +#else + num_gen_coords = 20; +#endif + + for (int i = 0; i < num_gen_coords; i++) { + std::cout << i << std::endl; + double point_x = dist (gen); + double point_y = dist (gen); + ref_coords.push_back (point_x); + ref_coords.push_back (point_y); + } + + const size_t num_coords = ref_coords.size () / 2; + std::vector out_coords (num_coords * 3, 0.0); + // Loop through all 12 base trees of the HEALPix projection + for (t8_gloidx_t gtreeid = 0; gtreeid < 12; ++gtreeid) { + geom.t8_geom_evaluate (cmesh, gtreeid, ref_coords.data (), num_coords, out_coords.data ()); + + for (size_t i = 0; i < num_coords; ++i) { + std::span point (out_coords.data () + 3 * i, 3); + double distance = t8_norm (point); + // Expect distance to be 1.0 (Unit Sphere) + EXPECT_NEAR (distance, 1.0, 1e-7) << "Radius check failed at Tree ID: " << gtreeid << " for point index: " << i + << " (xi=" << ref_coords[i * 2] << ", eta=" << ref_coords[i * 2 + 1] << ")"; + } + } + t8_cmesh_destroy (&cmesh); +}