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,158 @@
+/*
+  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) 2015 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::vector<double> getCoordsUpperRing(int k) {
+    if (k < 0 || k > 3) {
+        return {};
+    }
+
+    const double phi = M_PI / 2;
+    std::vector<double> coords(3);
+
+    coords[0] = cos(k * phi);
+    coords[1] = sin(k * phi);
+    coords[2] = 2.0 / 3;
+
+    return coords;
+}
+
+std::vector<double> getCoordsBottomRing (int k){
+    if (k < 0 || k > 3) {
+        return {};
+    }
+    std::vector <double> coords(3);
+    const double phi = M_PI / 2 ;
+    coords[0]= cos(k*phi);
+    coords[1]= sin(k*phi);
+    coords[2]= -2.0/3;
+    return coords;
+}
+
+
+std::vector<double> getCoordsEquator (int k){
+    if(k < 0 || k > 3) {
+        return {};
+    }
+    std::vector <double> coords(3);
+    const double phi = M_PI / 2 ;
+    const double shift = M_PI / 4;
+
+    coords[0]= cos(k*phi - shift);
+    coords[1]= sin(k*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 = T8_ALLOC (t8_eclass_t, 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::vector <double> northPole = {0,0,1};
+    std::vector <double> southPole = {0,0,-1};
+    int itree = 0;
+    // every side section has 4 quads, this builds the upper layer
+    // build upper section
+for (int k = 0; k < 4; k++) {
+    verts.clear();
+
+    auto equator = getCoordsEquator((k+1) % 4);
+    auto upper0  = getCoordsUpperRing(k % 4);
+    auto upper1  = getCoordsUpperRing((k + 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 k = 0; k < 4; k++) {
+    verts.clear();
+
+    auto bottom0  = getCoordsBottomRing(k % 4);
+    auto equator1 = getCoordsEquator((k + 1) % 4);
+    auto equator0 = getCoordsEquator(k % 4);
+    auto upper0   = getCoordsUpperRing(k % 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 k = 0; k < 4; k++) {
+    verts.clear();
+
+    auto bottom1 = getCoordsBottomRing((k + 1) % 4);
+    auto bottom0 = getCoordsBottomRing(k % 4);
+    auto equator = getCoordsEquator((k+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,26 @@
+/*
+  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) 2015 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>
+
+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) 2015 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,
+                                       [[maybe_unused]] t8_gloidx_t gtreeid, const double *ref_coords,
+                                       const size_t num_coords, double *out_coords) const
+{
+    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,90 @@
+/*
+  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) 2023 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>
+
+
+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")
+  {
+    std::cout<<"done1"<<std::endl;
+  }
+
+  /* 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 six hexaeders arranged into a cube.
+   *
+   */
+  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. */
+};