feat: implements healpix geometry

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

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 <cmath>
+#include <vector>
+#include <t8_cmesh/t8_cmesh.hxx>
+#include <t8_cmesh/t8_cmesh_helpers.h>
+#include <t8_geometry/t8_geometry_implementations/t8_geometry_linear.hxx>
+#include <t8_cmesh/t8_cmesh_healpix/t8_geometry_healpix.hxx>
+
+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<double> all_verts;
+  all_verts.reserve (ntrees * nverts * 3);
+
+  /* Register geometry and retain the pointer */
+  t8_geometry_c *geom = t8_cmesh_register_geometry<t8_geometry_healpix> (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<double> 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;
+}

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 <t8_cmesh/t8_cmesh.h>
+/** 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);

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 <t8.h>
+#include <t8_geometry/t8_geometry_with_vertices.hxx>
+#include <t8_cmesh/t8_cmesh_healpix/t8_geometry_healpix.hxx>
+#include <numbers>
+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<double> (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);
+  }
+}

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 <t8.h>
+#include <t8_geometry/t8_geometry_with_vertices.hxx>
+/**
+* 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);

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)

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 <gtest/gtest.h>
+#include <t8_cmesh/t8_cmesh.h>
+#include <t8_types/t8_vec.hxx>
+#include <t8_cmesh/t8_cmesh_healpix/t8_geometry_healpix.hxx>
+#include <t8_cmesh/t8_cmesh_healpix/t8_cmesh_healpix.h>
+#include <test/t8_gtest_macros.hxx>
+#include <random>
+/* 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<double> dist (0.0, 1.0);
+
+  // A set of 2D reference coordinates within the element
+  std::vector<double> 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<double> 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<const double> 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);
+}