diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index f51e1f050..f3c581ce6 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -194,6 +194,12 @@ jobs: ${ATLAS_TOOLS}/install-fftw.sh --version 3.3.10 --prefix ${DEPS_DIR}/fftw echo "FFTW_ROOT=${DEPS_DIR}/fftw" >> $GITHUB_ENV + - name: Install Qhull + shell: bash -eux {0} + run: | + ${ATLAS_TOOLS}/install-qhull.sh --version 8.1-alpha3 --prefix ${DEPS_DIR}/qhull + echo "Qhull_ROOT=${DEPS_DIR}/qhull" >> $GITHUB_ENV + - name: Install LZ4 if: "!contains( matrix.compiler, 'nvhpc' )" run: | diff --git a/cmake/features/TESSELATION.cmake b/cmake/features/TESSELATION.cmake index face9c0e4..3768e4089 100644 --- a/cmake/features/TESSELATION.cmake +++ b/cmake/features/TESSELATION.cmake @@ -1,16 +1,32 @@ if( atlas_HAVE_ATLAS_FUNCTIONSPACE ) ### tesselation ... -set(Boost_USE_MULTITHREADED ON ) - ecbuild_add_option( FEATURE TESSELATION + DESCRIPTION "Support for unstructured mesh generation" + CONDITION atlas_HAVE_ATLAS_FUNCTIONSPACE + REQUIRED_PACKAGES "Qhull" ) +if(HAVE_TESSELATION) + set(QHULL_LIBRARIES Qhull::qhullcpp Qhull::qhull_r) + set(atlas_HAVE_QHULL 1) +else() + set(atlas_HAVE_QHULL 0) +endif() + +### NOTE +# +# CGAL is deprecated as TESSELATION backend. Qhull is to be used instead. +# To use CGAL regardless, turn ON CGAL feature (-DENABLE_CGAL=ON) + +set(Boost_USE_MULTITHREADED ON ) +ecbuild_add_option( FEATURE CGAL + DEFAULT OFF DESCRIPTION "Support for unstructured mesh generation" CONDITION atlas_HAVE_ATLAS_FUNCTIONSPACE REQUIRED_PACKAGES - "CGAL QUIET" + "CGAL" "Boost VERSION 1.45.0 QUIET" ) -if( HAVE_TESSELATION ) +if( HAVE_CGAL ) list( APPEND CGAL_INCLUDE_DIRS ${Boost_INCLUDE_DIRS} ) if ( TARGET CGAL::CGAL ) list( APPEND CGAL_LIBRARIES CGAL::CGAL ${CGAL_3RD_PARTY_LIBRARIES} ${GMP_LIBRARIES} ${MPFR_LIBRARIES} ${Boost_THREAD_LIBRARY} ${Boost_SYSTEM_LIBRARY} ) @@ -27,8 +43,8 @@ if( HAVE_TESSELATION ) endif() endif() -if( NOT HAVE_TESSELATION ) +if( NOT HAVE_CGAL ) unset( CGAL_LIBRARIES ) unset( CGAL_INCLUDE_DIRS ) endif() -endif() \ No newline at end of file +endif() diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 14a7087d4..66cfb26df 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -24,6 +24,18 @@ else() set( atlas_HAVE_TESSELATION 0 ) endif() +if( atlas_HAVE_QHULL ) + set( atlas_HAVE_QHULL 1 ) +else() + set( atlas_HAVE_QHULL 0 ) +endif() + +if( atlas_HAVE_CGAL ) + set( atlas_HAVE_CGAL 1 ) +else() + set( atlas_HAVE_CGAL 0 ) +endif() + if( atlas_HAVE_PROJ ) set( atlas_HAVE_PROJ 1 ) else() diff --git a/src/atlas/CMakeLists.txt b/src/atlas/CMakeLists.txt index 11fb5f756..dad416859 100644 --- a/src/atlas/CMakeLists.txt +++ b/src/atlas/CMakeLists.txt @@ -819,6 +819,10 @@ util/vector.h util/mdspan.h util/detail/mdspan/mdspan.hpp util/VectorOfAbstract.h +util/QhullSphericalTriangulation.h +util/QhullSphericalTriangulation.cc +util/CGALSphericalTriangulation.h +util/CGALSphericalTriangulation.cc util/detail/Cache.h util/detail/KDTree.h util/function/MDPI_functions.h @@ -932,6 +936,7 @@ ecbuild_add_library( TARGET atlas ${CGAL_LIBRARIES} ${FFTW_LIBRARIES} ${PROJ_LIBRARIES} + ${QHULL_LIBRARIES} PUBLIC_LIBS eckit diff --git a/src/atlas/library/defines.h.in b/src/atlas/library/defines.h.in index 51d55a733..ec71db406 100644 --- a/src/atlas/library/defines.h.in +++ b/src/atlas/library/defines.h.in @@ -18,6 +18,8 @@ #define ATLAS_OMP_TASK_SUPPORTED @ATLAS_OMP_TASK_SUPPORTED@ #define ATLAS_OMP_TASK_UNTIED_SUPPORTED @ATLAS_OMP_TASK_UNTIED_SUPPORTED@ #define ATLAS_HAVE_ACC @atlas_HAVE_ACC@ +#define ATLAS_HAVE_QHULL @atlas_HAVE_QHULL@ +#define ATLAS_HAVE_CGAL @atlas_HAVE_CGAL@ #define ATLAS_HAVE_TESSELATION @atlas_HAVE_TESSELATION@ #define ATLAS_HAVE_FORTRAN @atlas_HAVE_FORTRAN@ #define ATLAS_HAVE_EIGEN @atlas_HAVE_EIGEN@ diff --git a/src/atlas/mesh/actions/BuildConvexHull3D.cc b/src/atlas/mesh/actions/BuildConvexHull3D.cc index 1d24058cf..756f20bbc 100644 --- a/src/atlas/mesh/actions/BuildConvexHull3D.cc +++ b/src/atlas/mesh/actions/BuildConvexHull3D.cc @@ -13,32 +13,10 @@ #include #include #include "eckit/log/BigNum.h" - +#include "atlas/util/QhullSphericalTriangulation.h" +#include "atlas/util/CGALSphericalTriangulation.h" #include "atlas/library/config.h" -#if ATLAS_HAVE_TESSELATION -// CGAL needs -DCGAL_NDEBUG to reach peak performance ... -#define CGAL_NDEBUG - -#include -#include -#include -#include -#include -#include - -using K = CGAL::Exact_predicates_inexact_constructions_kernel; -using Polyhedron_3 = CGAL::Polyhedron_3; - -using Vector_3 = K::Vector_3; -using FT = K::FT; -using Segment_3 = K::Segment_3; -using Point_3 = K::Point_3; - -const Point_3 origin = Point_3(CGAL::ORIGIN); - -#endif - #include "atlas/array/ArrayView.h" #include "atlas/array/MakeView.h" #include "atlas/field/Field.h" @@ -53,9 +31,6 @@ const Point_3 origin = Point_3(CGAL::ORIGIN); #include "atlas/runtime/Trace.h" #include "atlas/util/CoordinateEnums.h" -using namespace eckit::geometry; - -using atlas::interpolation::method::PointIndex3; using atlas::interpolation::method::PointSet; namespace atlas { @@ -66,156 +41,92 @@ namespace actions { BuildConvexHull3D::BuildConvexHull3D(const eckit::Parametrisation& config) { config.get("remove_duplicate_points", remove_duplicate_points_ = true); - config.get("reshuffle", reshuffle_ = true); } //---------------------------------------------------------------------------------------------------------------------- -#if ATLAS_HAVE_TESSELATION - -static Polyhedron_3* create_convex_hull_from_points(const std::vector& pts) { - ATLAS_TRACE(); - - Polyhedron_3* poly = new Polyhedron_3(); - - // insertion from a vector : - - std::vector vertices(pts.size()); - for (idx_t i = 0, size = vertices.size(); i < size; ++i) { - vertices[i] = Point_3(pts[i](XX), pts[i](YY), pts[i](ZZ)); - } - - // compute convex hull of non-collinear points - - CGAL::convex_hull_3(vertices.begin(), vertices.end(), *poly); - - return poly; -} - -static void cgal_polyhedron_to_atlas_mesh(Mesh& mesh, Polyhedron_3& poly, PointSet& points) { - ATLAS_TRACE(); - - bool ensure_outward_normals = true; - - mesh::Nodes& nodes = mesh.nodes(); - - ATLAS_ASSERT(points.size() == size_t(nodes.size())); - - const idx_t nb_nodes = idx_t(points.size()); - - ATLAS_ASSERT(mesh.cells().size() == 0); - - /* triangles */ - - const idx_t nb_triags = poly.size_of_facets(); - mesh.cells().add(mesh::ElementType::create("Triangle"), nb_triags); - mesh::HybridElements::Connectivity& triag_nodes = mesh.cells().node_connectivity(); - array::ArrayView triag_gidx = array::make_view(mesh.cells().global_index()); - array::ArrayView triag_part = array::make_view(mesh.cells().partition()); - - Point3 pt; - idx_t idx[3]; - Polyhedron_3::Vertex_const_handle vts[3]; - - Log::debug() << "Inserting triags (" << eckit::BigNum(nb_triags) << ")" << std::endl; - - idx_t tidx = 0; - for (Polyhedron_3::Facet_const_iterator f = poly.facets_begin(); f != poly.facets_end(); ++f) { - // loop over half-edges and take each vertex() - - idx_t iedge = 0; - Polyhedron_3::Halfedge_around_facet_const_circulator edge = f->facet_begin(); - do { - Polyhedron_3::Vertex_const_handle vh = edge->vertex(); - const Polyhedron_3::Point_3& p = vh->point(); - - pt.assign(p); - - idx[iedge] = points.unique(pt); - - ATLAS_ASSERT(idx[iedge] < nb_nodes); - - vts[iedge] = vh; - - ++iedge; - ++edge; - } while (edge != f->facet_begin() && iedge < 3); - - ATLAS_ASSERT(iedge == 3); - - if (ensure_outward_normals) /* ensure outward pointing normal */ - { - Vector_3 p0(origin, vts[0]->point()); - Vector_3 n = CGAL::normal(vts[0]->point(), vts[1]->point(), vts[2]->point()); - - FT innerp = n * p0; - - if (innerp < 0) { // need to swap an edge of the triag - std::swap(vts[1], vts[2]); - } - } - - /* define the triag */ - - triag_nodes.set(tidx, idx); - - triag_gidx(tidx) = tidx + 1; - - triag_part(tidx) = 0; - - ++tidx; - } - - ATLAS_ASSERT(tidx == nb_triags); -} - -#else - -struct Polyhedron_3 { - idx_t size_of_vertices() const { return 0; } -}; - -static Polyhedron_3* create_convex_hull_from_points(const std::vector& pts) { - throw_NotImplemented("CGAL package not found -- Delaunay triangulation is disabled", Here()); -} - -static void cgal_polyhedron_to_atlas_mesh(Mesh& mesh, Polyhedron_3& poly, PointSet& points) { - throw_NotImplemented("CGAL package not found -- Delaunay triangulation is disabled", Here()); -} - -#endif - -//---------------------------------------------------------------------------------------------------------------------- - void BuildConvexHull3D::operator()(Mesh& mesh) const { // don't tesselate meshes already with triags or quads if (mesh.cells().size()) { return; } - ATLAS_TRACE(); - - // remove duplicate points - - PointSet points(mesh); - - std::vector ipts; - - points.list_unique_points(ipts); + std::string default_backend = (ATLAS_HAVE_CGAL ? "cgal" : "qhull"); + std::string backend = eckit::Resource("$ATLAS_DELAUNAY_BACKEND",default_backend); - // std::cout << "unique pts " << ipts.size() << std::endl; - // std::cout << "duplicates " << points.duplicates().size() << std::endl; + ATLAS_TRACE("BuildConvexHull3D [" + backend + "]"); - // define polyhedron to hold convex hull - - std::unique_ptr poly(create_convex_hull_from_points(ipts)); + std::vector local_index; + if (remove_duplicate_points_) { + PointSet points(mesh); + points.list_unique_points(local_index); + } + + auto add_triangles = [&]( const auto& triangles ) { + const idx_t nb_triags = triangles.size(); + mesh.cells().add(mesh::ElementType::create("Triangle"), nb_triags); + mesh::HybridElements::Connectivity& triag_nodes = mesh.cells().node_connectivity(); + auto triag_gidx = array::make_view(mesh.cells().global_index()); + auto triag_part = array::make_view(mesh.cells().partition()); + + Log::debug() << "Inserting triags (" << eckit::BigNum(nb_triags) << ")" << std::endl; + + for (idx_t tidx = 0; tidx idx{t[0],t[1],t[2]}; + if( local_index.size() ) { + idx[0] = local_index[idx[0]]; + idx[1] = local_index[idx[1]]; + idx[2] = local_index[idx[2]]; + } + triag_nodes.set(tidx, idx.data()); + triag_gidx(tidx) = tidx + 1; + triag_part(tidx) = 0; + } + }; - // std::cout << "convex hull " << poly->size_of_vertices() << " vertices" - // << std::endl; - ATLAS_ASSERT(poly->size_of_vertices() == static_cast(ipts.size())); + if( local_index.size() == mesh.nodes().size() or local_index.empty() ) { + local_index.clear(); + auto lonlat = array::make_view(mesh.nodes().lonlat()); + if (backend == "qhull") { + ATLAS_TRACE("qhull"); + auto triangles = util::QhullSphericalTriangulation{static_cast(lonlat.shape(0)),lonlat.data()}.triangles(); + add_triangles(triangles); + } + else if (backend == "cgal") { + ATLAS_TRACE("cgal"); + auto triangles = util::CGALSphericalTriangulation{static_cast(lonlat.shape(0)),lonlat.data()}.triangles(); + add_triangles(triangles); + } + else { + ATLAS_THROW_EXCEPTION("backend " << backend << " not supported"); + } + } + else { + auto lonlat_view = array::make_view(mesh.nodes().lonlat()); + + std::vector lonlat(local_index.size()); + size_t jnode = 0; + for( auto& ip: local_index ) { + lonlat[jnode] = {lonlat_view(ip,0),lonlat_view(ip,1)}; + ++jnode; + } + if (backend == "qhull") { + ATLAS_TRACE("qhull"); + auto triangles = util::QhullSphericalTriangulation{lonlat}.triangles(); + add_triangles(triangles); + } + else if (backend == "cgal") { + ATLAS_TRACE("cgal"); + auto triangles = util::CGALSphericalTriangulation{lonlat}.triangles(); + add_triangles(triangles); + } + else { + ATLAS_THROW_EXCEPTION("backend " << backend << " not supported"); + } + } - cgal_polyhedron_to_atlas_mesh(mesh, *poly, points); } //---------------------------------------------------------------------------------------------------------------------- diff --git a/src/atlas/util/CGALSphericalTriangulation.cc b/src/atlas/util/CGALSphericalTriangulation.cc new file mode 100644 index 000000000..f3d563cfa --- /dev/null +++ b/src/atlas/util/CGALSphericalTriangulation.cc @@ -0,0 +1,193 @@ +#include "CGALSphericalTriangulation.h" + +#include "atlas/library/defines.h" + +#include +#include +#include +#include + +#include "atlas/runtime/Exception.h" + +#if ATLAS_HAVE_CGAL + +// CGAL needs -DCGAL_NDEBUG to reach peak performance ... +#define CGAL_NDEBUG + +#include +#include +#include +#include +#include +#include + + +using K = ::CGAL::Exact_predicates_inexact_constructions_kernel; +using Polyhedron_3 = ::CGAL::Polyhedron_3; +using Point_3 = Polyhedron_3::Point_3; + +#endif + + +namespace atlas{ +namespace util{ + +#if ATLAS_HAVE_CGAL +struct CGALSphericalTriangulation::CGAL { + CGAL(const std::vector>& pts) { + std::vector vertices(pts.size()); + for (size_t i = 0, size = vertices.size(); i < size; ++i) { + vertices[i] = Point_3(pts[i][0], pts[i][1], pts[i][2]); + point_map_[vertices[i]] = i; + } + + // compute convex hull of non-collinear points + ::CGAL::convex_hull_3(vertices.begin(), vertices.end(), poly_); + } + std::unordered_map point_map_; + Polyhedron_3 poly_; +}; +#else +struct CGALSphericalTriangulation::CGAL { + template + CGAL(Args...) { + throw_Exception("Atlas has not been compiled with CGAL",Here()); + } +}; +#endif + +CGALSphericalTriangulation::~CGALSphericalTriangulation() = default; + +CGALSphericalTriangulation::CGALSphericalTriangulation(size_t N, const double lonlat[]) : + CGALSphericalTriangulation(N, lonlat, lonlat+1, 2, 2) { +} + +CGALSphericalTriangulation::CGALSphericalTriangulation(size_t N, const double lon[], const double lat[]) : + CGALSphericalTriangulation(N, lon, lat, 1, 1) { +} + +CGALSphericalTriangulation::CGALSphericalTriangulation(size_t N, const double lon[], const double lat[], int lon_stride, int lat_stride) { + auto lonlat2xyz = [](double lon, double lat, auto& xyz) { + constexpr double deg2rad = M_PI / 180.; + const double lambda = deg2rad * lon; + const double phi = deg2rad * lat; + + const double sin_phi = std::sin(phi); + const double cos_phi = std::cos(phi); + const double sin_lambda = std::sin(lambda); + const double cos_lambda = std::cos(lambda); + + xyz[0] = cos_phi * cos_lambda; + xyz[1] = cos_phi * sin_lambda; + xyz[2] = sin_phi; + }; + + points_xyz_.resize(N); + + for (size_t i = 0; i < N; ++i) { + lonlat2xyz(lon[i*lon_stride], lat[i*lat_stride], points_xyz_[i]); + } + + cgal_ = std::make_unique(points_xyz_); +} + +size_t CGALSphericalTriangulation::size() const { +#if ATLAS_HAVE_CGAL + return cgal_->poly_.size_of_facets(); +#else + return 0; +#endif +} + + +template +static inline void CGAL_get_triangles(const CGAL& cgal, const Points& points_xyz_, std::array triangles[]) { +#if ATLAS_HAVE_CGAL + auto ensure_outward_normal = [&](auto& tri) { + + auto dot = [](const auto& p1, const auto& p2) { + return p1[0]*p2[0] + p1[1]*p2[1] + p1[2]*p2[2]; + }; + auto cross = [](const auto& p1, const auto& p2) { + return std::array{ + p1[1] * p2[2] - p1[2] * p2[1], p1[2] * p2[0] - p1[0] * p2[2], + p1[0] * p2[1] - p1[1] * p2[0] + }; + }; + + const std::array& a = points_xyz_[tri[0]]; + const std::array& b = points_xyz_[tri[1]]; + const std::array& c = points_xyz_[tri[2]]; + std::array ba {a[0]-b[0], a[1]-b[1], a[2]-b[2]}; + std::array bc {c[0]-b[0], c[1]-b[1], c[2]-b[2]}; + + bool outward = dot(b, cross(bc,ba)) > 0; + + if (not outward) { + std::swap(tri[1], tri[2]); + } + }; + + size_t jtri{0}; + + const auto& poly = cgal.poly_; + const auto& point_map = cgal.point_map_; + const idx_t nb_points = idx_t(points_xyz_.size()); + + /* triangles */ + + const idx_t nb_triags = poly.size_of_facets(); + + idx_t idx[3]; + Polyhedron_3::Vertex_const_handle vts[3]; + + idx_t tidx = 0; + for (Polyhedron_3::Facet_const_iterator f = poly.facets_begin(); f != poly.facets_end(); ++f) { + // loop over half-edges and take each vertex() + + auto& tri = triangles[jtri++]; + + size_t jvrt{0}; + + idx_t iedge = 0; + Polyhedron_3::Halfedge_around_facet_const_circulator edge = f->facet_begin(); + do { + Polyhedron_3::Vertex_const_handle vh = edge->vertex(); + const Polyhedron_3::Point_3& p = vh->point(); + tri[jvrt++] = point_map.at(vh->point()); + + ++iedge; + ++edge; + } while (edge != f->facet_begin() && iedge < 3); + + ATLAS_ASSERT(iedge == 3); + + ensure_outward_normal(tri); + } + + ATLAS_ASSERT(jtri == nb_triags); +#endif +} + +void CGALSphericalTriangulation::triangles(std::array triangles[]) const { + CGAL_get_triangles(*cgal_,points_xyz_,triangles); +} + +void CGALSphericalTriangulation::triangles(std::array triangles[]) const { + CGAL_get_triangles(*cgal_,points_xyz_,triangles); +} + +void CGALSphericalTriangulation::triangles(std::array triangles[]) const { + CGAL_get_triangles(*cgal_,points_xyz_,triangles); +} + +void CGALSphericalTriangulation::triangles(std::array triangles[]) const { + CGAL_get_triangles(*cgal_,points_xyz_,triangles); +} + +std::vector> CGALSphericalTriangulation::triangles() const { + return triangles(); +} + +} // namespace util +} // namespace atlas diff --git a/src/atlas/util/CGALSphericalTriangulation.h b/src/atlas/util/CGALSphericalTriangulation.h new file mode 100644 index 000000000..c793689c4 --- /dev/null +++ b/src/atlas/util/CGALSphericalTriangulation.h @@ -0,0 +1,64 @@ +/* + * (C) Copyright 2023 ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation + * nor does it submit to any jurisdiction. + */ + +#pragma once + +#include +#include +#include + +#include "atlas/library/config.h" +namespace atlas { +namespace util { + + +class CGALSphericalTriangulation { +public: + template + CGALSphericalTriangulation(const LonLatVector& lonlat); + CGALSphericalTriangulation(size_t N, const double lonlat[]); + CGALSphericalTriangulation(size_t N, const double lon[], const double lat[]); + CGALSphericalTriangulation(size_t N, const double lon[], const double lat[], int lon_stride, int lat_stride); + + ~CGALSphericalTriangulation(); + + /// @return number of triangles + size_t size() const; + + void triangles(std::array[]) const; + void triangles(std::array[]) const; + void triangles(std::array[]) const; + void triangles(std::array[]) const; + + template + std::vector> triangles() const; + std::vector> triangles() const; + +private: + struct CGAL; + std::unique_ptr cgal_; + std::vector> points_xyz_; +}; + +template +inline CGALSphericalTriangulation::CGALSphericalTriangulation(const LonLatVector& lonlat) : + CGALSphericalTriangulation(lonlat.size(), reinterpret_cast(lonlat.data())) {} + +template +inline std::vector> CGALSphericalTriangulation::triangles() const { + std::vector> vector(size()); + triangles(vector.data()); + return vector; +} + + +} +} + diff --git a/src/atlas/util/QhullSphericalTriangulation.cc b/src/atlas/util/QhullSphericalTriangulation.cc new file mode 100644 index 000000000..d0ec609de --- /dev/null +++ b/src/atlas/util/QhullSphericalTriangulation.cc @@ -0,0 +1,146 @@ +#include "QhullSphericalTriangulation.h" + +#include "atlas/library/defines.h" + +#include +#include +#include +#include + +#if ATLAS_HAVE_QHULL +#include +#include +#include +#include +#include +#endif + +#include "atlas/runtime/Exception.h" + +namespace atlas{ +namespace util{ + +#if ATLAS_HAVE_QHULL +struct QhullSphericalTriangulation::Qhull : public orgQhull::Qhull { + using orgQhull::Qhull::Qhull; +}; +#else +struct QhullSphericalTriangulation::Qhull { + template + Qhull(Args...) { + throw_Exception("Atlas has not been compiled with Qhull",Here()); + } +}; +#endif + +QhullSphericalTriangulation::~QhullSphericalTriangulation() = default; + +QhullSphericalTriangulation::QhullSphericalTriangulation(size_t N, const double lonlat[]) : + QhullSphericalTriangulation(N, lonlat, lonlat+1, 2, 2) { +} + +QhullSphericalTriangulation::QhullSphericalTriangulation(size_t N, const double lon[], const double lat[]) : + QhullSphericalTriangulation(N, lon, lat, 1, 1) { +} + +QhullSphericalTriangulation::QhullSphericalTriangulation(size_t N, const double lon[], const double lat[], int lon_stride, int lat_stride) { + auto lonlat2xyz = [](double lon, double lat, auto& xyz) { + constexpr double deg2rad = M_PI / 180.; + const double lambda = deg2rad * lon; + const double phi = deg2rad * lat; + + const double sin_phi = std::sin(phi); + const double cos_phi = std::cos(phi); + const double sin_lambda = std::sin(lambda); + const double cos_lambda = std::cos(lambda); + + xyz[0] = cos_phi * cos_lambda; + xyz[1] = cos_phi * sin_lambda; + xyz[2] = sin_phi; + }; + + points_xyz_.resize(N); + + for (size_t i = 0; i < N; ++i) { + lonlat2xyz(lon[i*lon_stride], lat[i*lat_stride], points_xyz_[i]); + } + + constexpr int dim = 3; + constexpr const char* command = "Qt"; + constexpr const char* comment = ""; + const double* coordinates = reinterpret_cast(points_xyz_.data()); + qhull_ = std::make_unique(comment, dim, N, coordinates, command); +} + +size_t QhullSphericalTriangulation::size() const { +#if ATLAS_HAVE_QHULL + return qhull_->facetList().size(); +#else + return 0; +#endif +} + + +template +static inline void qhull_get_triangles(const Qhull& qhull, const Points& points_xyz_, std::array triangles[]) { +#if ATLAS_HAVE_QHULL + auto ensure_outward_normal = [&](auto& tri) { + + auto dot = [](const auto& p1, const auto& p2) { + return p1[0]*p2[0] + p1[1]*p2[1] + p1[2]*p2[2]; + }; + auto cross = [](const auto& p1, const auto& p2) { + return std::array{ + p1[1] * p2[2] - p1[2] * p2[1], p1[2] * p2[0] - p1[0] * p2[2], + p1[0] * p2[1] - p1[1] * p2[0] + }; + }; + + const std::array& a = points_xyz_[tri[0]]; + const std::array& b = points_xyz_[tri[1]]; + const std::array& c = points_xyz_[tri[2]]; + std::array ba {a[0]-b[0], a[1]-b[1], a[2]-b[2]}; + std::array bc {c[0]-b[0], c[1]-b[1], c[2]-b[2]}; + + bool outward = dot(b, cross(bc,ba)) > 0; + + if (not outward) { + std::swap(tri[1], tri[2]); + } + }; + + size_t jtri{0}; + for (const auto& facet : qhull.facetList()){ + auto& tri = triangles[jtri++]; + size_t jvrt{0}; + for( const auto& vertex: facet.vertices()){ + tri[jvrt++] = vertex.point().id(); + } + + ensure_outward_normal(tri); + } +#endif +} + +void QhullSphericalTriangulation::triangles(std::array triangles[]) const { + qhull_get_triangles(*qhull_,points_xyz_,triangles); +} + +void QhullSphericalTriangulation::triangles(std::array triangles[]) const { + qhull_get_triangles(*qhull_,points_xyz_,triangles); +} + +void QhullSphericalTriangulation::triangles(std::array triangles[]) const { + qhull_get_triangles(*qhull_,points_xyz_,triangles); +} + +void QhullSphericalTriangulation::triangles(std::array triangles[]) const { + qhull_get_triangles(*qhull_,points_xyz_,triangles); +} + +std::vector> QhullSphericalTriangulation::triangles() const { + return triangles(); +} + +} // namespace util +} // namespace atlas diff --git a/src/atlas/util/QhullSphericalTriangulation.h b/src/atlas/util/QhullSphericalTriangulation.h new file mode 100644 index 000000000..f5beff754 --- /dev/null +++ b/src/atlas/util/QhullSphericalTriangulation.h @@ -0,0 +1,63 @@ +/* + * (C) Copyright 2023 ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation + * nor does it submit to any jurisdiction. + */ + +#pragma once + +#include +#include +#include + +#include "atlas/library/config.h" +namespace atlas { +namespace util { + +class QhullSphericalTriangulation { +public: + template + QhullSphericalTriangulation(const LonLatVector& lonlat); + QhullSphericalTriangulation(size_t N, const double lonlat[]); + QhullSphericalTriangulation(size_t N, const double lon[], const double lat[]); + QhullSphericalTriangulation(size_t N, const double lon[], const double lat[], int lon_stride, int lat_stride); + + ~QhullSphericalTriangulation(); + + /// @return number of triangles + size_t size() const; + + void triangles(std::array[]) const; + void triangles(std::array[]) const; + void triangles(std::array[]) const; + void triangles(std::array[]) const; + + template + std::vector> triangles() const; + std::vector> triangles() const; + +private: + struct Qhull; + std::unique_ptr qhull_; + std::vector> points_xyz_; +}; + +template +inline QhullSphericalTriangulation::QhullSphericalTriangulation(const LonLatVector& lonlat) : + QhullSphericalTriangulation(lonlat.size(), reinterpret_cast(lonlat.data())) {} + +template +inline std::vector> QhullSphericalTriangulation::triangles() const { + std::vector> vector(size()); + triangles(vector.data()); + return vector; +} + + +} +} + diff --git a/tools/install-qhull.sh b/tools/install-qhull.sh new file mode 100755 index 000000000..eed700402 --- /dev/null +++ b/tools/install-qhull.sh @@ -0,0 +1,87 @@ +#! /usr/bin/env bash + +# (C) Copyright 2013 ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation nor +# does it submit to any jurisdiction. + +set +x +set -e -o pipefail + +SCRIPTDIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd )" +export PATH=$SCRIPTDIR:$PATH + +# Some defaults for the arguments +PREFIX=$(pwd)/install +qhull_version=8.1-alpha3 + +while [ $# != 0 ]; do + case "$1" in + "--prefix") + PREFIX="$2"; shift + ;; + "--version") + qhull_version="$2"; shift + ;; + *) + echo "Unrecognized argument '$1'" + exit 1 + ;; + esac + shift +done + +echo "Installing Qhull version ${qhull_version}" + +qhull_installed=${PREFIX}/qhull-${qhull_version}-installed +if [[ -f "${qhull_installed}" ]]; then + echo "Qhull ${qhull_version} is already installed at ${PREFIX}" + exit +fi + +os=$(uname) +case "$os" in + Darwin) + brew install qhull + exit + ;; + *) + ;; +esac + + +if [ -z "${TMPDIR+x}" ]; then + TMPDIR=${HOME}/tmp +fi +mkdir -p ${TMPDIR}/downloads + +qhull_tarball_url=https://github.com/qhull/qhull/archive/refs/tags/v${qhull_version}.tar.gz +qhull_tarball=$TMPDIR/downloads/qhull-v${qhull_version}.tar.gz +qhull_dir=$TMPDIR/downloads/qhull-${qhull_version} + +echo "+ curl -L ${qhull_tarball_url} > ${qhull_tarball}" +curl -L ${qhull_tarball_url} > ${qhull_tarball} +echo "+ tar xzf ${qhull_tarball} -C ${TMPDIR}/downloads" +tar xzf ${qhull_tarball} -C ${TMPDIR}/downloads +echo "+ cd ${qhull_dir}" +cd ${qhull_dir} +#echo "+ ./configure --prefix=${PREFIX} ${fftw_configure}" +#./configure --prefix=${PREFIX} ${fftw_configure} +echo "+ mkdir build-release" +mkdir build-release +echo "+ cd build-release" +cd build-release +cmake .. -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=${PREFIX} +echo "+ make -j8" +make -j8 +echo "+ make install" +make install + +echo "+ rm -rf \${qhull_tarball} \${qhull_dir}" +rm -rf ${qhull_tarball} ${qhull_dir} + +echo "+ touch ${qhull_installed}" +touch ${qhull_installed}