feat: implements healpix geometry

src/CMakeLists.txt

@@ -150,6 +150,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/t8_cmesh_save.cxx
     t8_cmesh/t8_cmesh_io/t8_cmesh_triangle.cxx

src/t8_cmesh/t8_cmesh_healpix/t8_cmesh_healpix.cxx

@@ -0,0 +1,163 @@
+/*
+  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>
+
+std::array<double, 3>
+getCoordsUpperRing (int side)
+{
+  if (side < 0 || side > 3) {
+    return {};
+  }
+
+  const double phi = M_PI / 2;
+  std::array<double, 3> coords;
+
+  coords[0] = cos (side * phi);
+  coords[1] = sin (side * phi);
+  coords[2] = 2.0 / 3;
+
+  return coords;
+}
+
+std::array<double, 3>
+getCoordsBottomRing (int side)
+{
+  if (side < 0 || side > 3) {
+    return {};
+  }
+  std::array<double, 3> coords;
+  const double phi = M_PI / 2;
+  coords[0] = cos (side * phi);
+  coords[1] = sin (side * phi);
+  coords[2] = -2.0 / 3;
+  return coords;
+}
+
+std::array<double, 3>
+getCoordsEquator (int side)
+{
+  if (side < 0 || side > 3) {
+    return {};
+  }
+  std::array<double, 3> coords;
+  const double phi = M_PI / 2;
+  const double shift = M_PI / 4;
+
+  coords[0] = cos (side * phi - shift);
+  coords[1] = sin (side * phi - shift);
+  coords[2] = 0;
+  return coords;
+}
+
+t8_cmesh_t
+t8_cmesh_new_healpix (sc_MPI_Comm comm)
+{
+  /* Initialization of the mesh */
+  t8_cmesh_t cmesh;
+  t8_cmesh_init (&cmesh);
+
+  // t8_cmesh_register_geometry<t8_geometry_linear> (cmesh); /* Use spherical geometry. */
+
+  const int ntrees = 12;
+  const int nverts = 4; /* Number of vertices per cmesh element. */
+  std::vector<double> verts;
+  t8_eclass_t all_eclasses[ntrees];
+  std::vector<double> all_verts;
+  t8_cmesh_register_geometry<t8_geometry_healpix> (cmesh);
+  /* Defitition of the tree class. */
+  for (int itree = 0; itree < ntrees; itree++) {
+    t8_cmesh_set_tree_class (cmesh, itree, T8_ECLASS_QUAD);
+    all_eclasses[itree] = T8_ECLASS_QUAD;
+  }
+  std::array<double, 3> northPole = { 0, 0, 1 };
+  std::array<double, 3> southPole = { 0, 0, -1 };
+  int itree = 0;
+  // every side section has 4 quads, this builds the upper layer
+  // build upper section
+  for (int side = 0; side < 4; side++) {
+    verts.clear ();
+
+    auto equator = getCoordsEquator ((side + 1) % 4);
+    auto upper0 = getCoordsUpperRing (side % 4);
+    auto upper1 = getCoordsUpperRing ((side + 1) % 4);
+
+    verts.insert (verts.end (), northPole.begin (), northPole.end ());
+    verts.insert (verts.end (), upper1.begin (), upper1.end ());
+    verts.insert (verts.end (), upper0.begin (), upper0.end ());
+    verts.insert (verts.end (), equator.begin (), equator.end ());
+
+    all_verts.insert (std::end (all_verts), std::begin (verts), std::end (verts));
+    t8_cmesh_set_tree_vertices (cmesh, itree, verts.data (), nverts);
+
+    itree++;
+  }
+
+  // build middle section
+  for (int side = 0; side < 4; side++) {
+    verts.clear ();
+
+    auto bottom0 = getCoordsBottomRing (side % 4);
+    auto equator1 = getCoordsEquator ((side + 1) % 4);
+    auto equator0 = getCoordsEquator (side % 4);
+    auto upper0 = getCoordsUpperRing (side % 4);
+
+    verts.insert (verts.end (), bottom0.begin (), bottom0.end ());
+    verts.insert (verts.end (), equator0.begin (), equator0.end ());
+    verts.insert (verts.end (), equator1.begin (), equator1.end ());
+    verts.insert (verts.end (), upper0.begin (), upper0.end ());
+
+    all_verts.insert (std::end (all_verts), std::begin (verts), std::end (verts));
+    t8_cmesh_set_tree_vertices (cmesh, itree, verts.data (), nverts);
+
+    itree++;
+  }
+
+  // build lower section
+  for (int side = 0; side < 4; side++) {
+    verts.clear ();
+
+    auto bottom1 = getCoordsBottomRing ((side + 1) % 4);
+    auto bottom0 = getCoordsBottomRing (side % 4);
+    auto equator = getCoordsEquator ((side + 1) % 4);
+
+    verts.insert (verts.end (), southPole.begin (), southPole.end ());
+    verts.insert (verts.end (), bottom0.begin (), bottom0.end ());
+    verts.insert (verts.end (), bottom1.begin (), bottom1.end ());
+    verts.insert (verts.end (), equator.begin (), equator.end ());
+
+    all_verts.insert (std::end (all_verts), std::begin (verts), std::end (verts));
+    t8_cmesh_set_tree_vertices (cmesh, itree, verts.data (), nverts);
+
+    itree++;
+  }
+
+  t8_cmesh_set_join_by_vertices (cmesh, 12, all_eclasses, all_verts.data (), nullptr, 0);
+  T8_FREE (all_eclasses);
+  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,76 @@
+/*
+  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>
+
+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;
+
+  // Calculate the offsets based on the specific layer and face of this tree.
+  const double offset_y = static_cast<double> (layer) - 1.0;
+  const double offset_x = 2.0 * face + ((layer + 1) % 2);
+
+  for (size_t i = 0; i < num_coords; ++i) {
+    size_t base_idx_2d = i * 2;
+    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);
+
+    // Transform reference coordinates into the HEALPix intermediate project space (dX, dY)
+    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_coords[base_idx_3d + 0] = sin_theta * std::cos (phi);
+    out_coords[base_idx_3d + 1] = sin_theta * std::sin (phi);
+    out_coords[base_idx_3d + 2] = z;
+  }
+}

src/t8_cmesh/t8_cmesh_healpix/t8_geometry_healpix.hxx

@@ -0,0 +1,93 @@
+/*
+  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. */
+};