Skip to content

Commit

Permalink
Implement Delaunay triangulation using Qhull instead of CGAL
Browse files Browse the repository at this point in the history
  • Loading branch information
@wdeconinck
wdeconinck committed Sep 12, 2023
1 parent 2b51852 commit a922f4b
Show file tree
Hide file tree
Showing 11 changed files with 673 additions and 168 deletions.
6 changes: 6 additions & 0 deletions .github/workflows/build.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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: |
Expand Down
28 changes: 22 additions & 6 deletions cmake/features/TESSELATION.cmake
View file
Original file line number Diff line number Diff line change
@@ -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} )
Expand All @@ -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()
endif()
12 changes: 12 additions & 0 deletions src/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand Down
5 changes: 5 additions & 0 deletions src/atlas/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -932,6 +936,7 @@ ecbuild_add_library( TARGET atlas
${CGAL_LIBRARIES}
${FFTW_LIBRARIES}
${PROJ_LIBRARIES}
${QHULL_LIBRARIES}

PUBLIC_LIBS
eckit
Expand Down
2 changes: 2 additions & 0 deletions src/atlas/library/defines.h.in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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@
Expand Down
235 changes: 73 additions & 162 deletions src/atlas/mesh/actions/BuildConvexHull3D.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,32 +13,10 @@
#include <memory>
#include <vector>
#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 <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/Real_timer.h>
#include <CGAL/algorithm.h>
#include <CGAL/convex_hull_3.h>
#include <CGAL/point_generators_3.h>

using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using Polyhedron_3 = CGAL::Polyhedron_3<K>;

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"
Expand All @@ -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 {
Expand All @@ -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<Point3>& pts) {
ATLAS_TRACE();

Polyhedron_3* poly = new Polyhedron_3();

// insertion from a vector :

std::vector<Point_3> 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<gidx_t, 1> triag_gidx = array::make_view<gidx_t, 1>(mesh.cells().global_index());
array::ArrayView<int, 1> triag_part = array::make_view<int, 1>(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<Point3>& 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<Point3> ipts;

points.list_unique_points(ipts);
std::string default_backend = (ATLAS_HAVE_CGAL ? "cgal" : "qhull");
std::string backend = eckit::Resource<std::string>("$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<Polyhedron_3> poly(create_convex_hull_from_points(ipts));
std::vector<size_t> 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<gidx_t, 1>(mesh.cells().global_index());
auto triag_part = array::make_view<int, 1>(mesh.cells().partition());

Log::debug() << "Inserting triags (" << eckit::BigNum(nb_triags) << ")" << std::endl;

for (idx_t tidx = 0; tidx<nb_triags; ++tidx) {
auto& t = triangles[tidx];
std::array<idx_t,3> 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<idx_t>(ipts.size()));
if( local_index.size() == mesh.nodes().size() or local_index.empty() ) {
local_index.clear();
auto lonlat = array::make_view<double,2>(mesh.nodes().lonlat());
if (backend == "qhull") {
ATLAS_TRACE("qhull");
auto triangles = util::QhullSphericalTriangulation{static_cast<size_t>(lonlat.shape(0)),lonlat.data()}.triangles();
add_triangles(triangles);
}
else if (backend == "cgal") {
ATLAS_TRACE("cgal");
auto triangles = util::CGALSphericalTriangulation{static_cast<size_t>(lonlat.shape(0)),lonlat.data()}.triangles();
add_triangles(triangles);
}
else {
ATLAS_THROW_EXCEPTION("backend " << backend << " not supported");
}
}
else {
auto lonlat_view = array::make_view<double,2>(mesh.nodes().lonlat());

std::vector<PointLonLat> 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);
}

//----------------------------------------------------------------------------------------------------------------------
Expand Down
Loading

0 comments on commit a922f4b

Please sign in to comment.