diff --git a/.github/ci-config.yml b/.github/ci-config.yml index 004e8ab62..b894b8d4d 100644 --- a/.github/ci-config.yml +++ b/.github/ci-config.yml @@ -1,5 +1,6 @@ dependencies: | ecmwf/ecbuild ecmwf/eckit + ecmwf/fckit dependency_branch: develop parallelism_factor: 8 diff --git a/.github/ci-hpc-config.yml b/.github/ci-hpc-config.yml index eab4a6cd2..3689b72e3 100644 --- a/.github/ci-hpc-config.yml +++ b/.github/ci-hpc-config.yml @@ -12,6 +12,7 @@ mpi_on: dependencies: - ecmwf/ecbuild@develop - ecmwf/eckit@develop + - ecmwf/fckit@develop parallel: 64 ntasks: 16 @@ -24,4 +25,5 @@ mpi_off: dependencies: - ecmwf/ecbuild@develop - ecmwf/eckit@develop + - ecmwf/fckit@develop parallel: 64 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/ECTRANS.cmake b/cmake/features/ECTRANS.cmake index 0d3a5355c..3d8ab3f08 100644 --- a/cmake/features/ECTRANS.cmake +++ b/cmake/features/ECTRANS.cmake @@ -2,13 +2,16 @@ if( atlas_HAVE_ATLAS_TRANS ) ### trans ... +if( NOT DEFINED ATLAS_ENABLE_ECTRANS AND DEFINED ATLAS_ENABLE_TRANS ) + ecbuild_warn("Atlas option ATLAS_ENABLE_TRANS is deprecated in favour of ATLAS_ENABLE_ECTRANS") + set( ATLAS_ENABLE_ECTRANS ${ATLAS_ENABLE_TRANS} ) +endif() if( NOT DEFINED ENABLE_ECTRANS AND DEFINED ENABLE_TRANS ) ecbuild_warn("Atlas option ENABLE_TRANS is deprecated in favour of ENABLE_ECTRANS") set( ENABLE_ECTRANS ${ENABLE_TRANS} ) endif() -if( NOT DEFINED ATLAS_ENABLE_ECTRANS AND DEFINED ATLAS_ENABLE_TRANS ) - ecbuild_warn("Atlas option ATLAS_ENABLE_TRANS is deprecated in favour of ATLAS_ENABLE_ECTRANS") - set( ATLAS_ENABLE_TRANS ${ENABLE_TRANS} ) +if( DEFINED ATLAS_ENABLE_ECTRANS ) + set( ENABLE_ECTRANS ${ATLAS_ENABLE_ECTRANS} ) endif() set( atlas_HAVE_PACKAGE_ECTRANS 0 ) @@ -32,4 +35,4 @@ ecbuild_add_option( FEATURE ECTRANS DESCRIPTION "Support for IFS spectral transforms" CONDITION atlas_HAVE_ATLAS_FUNCTIONSPACE AND transi_FOUND ) -endif() \ No newline at end of file +endif() 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 2f2fc213c..dad416859 100644 --- a/src/atlas/CMakeLists.txt +++ b/src/atlas/CMakeLists.txt @@ -341,6 +341,10 @@ mesh/Connectivity.cc mesh/Connectivity.h mesh/ElementType.cc mesh/ElementType.h +mesh/elementtypes/Classification.h +mesh/elementtypes/Triangle.h +mesh/elementtypes/Quadrilateral.h +mesh/elementtypes/Pentagon.h mesh/Elements.cc mesh/Elements.h mesh/Halo.cc @@ -815,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 @@ -928,6 +936,7 @@ ecbuild_add_library( TARGET atlas ${CGAL_LIBRARIES} ${FFTW_LIBRARIES} ${PROJ_LIBRARIES} + ${QHULL_LIBRARIES} PUBLIC_LIBS eckit diff --git a/src/atlas/functionspace/PointCloud.cc b/src/atlas/functionspace/PointCloud.cc index 1e6da5a58..2ce1d64c7 100644 --- a/src/atlas/functionspace/PointCloud.cc +++ b/src/atlas/functionspace/PointCloud.cc @@ -70,9 +70,22 @@ PointCloud::PointCloud(const Field& lonlat, const Field& ghost): lonlat_(lonlat) setupHaloExchange(); } -PointCloud::PointCloud(const FieldSet & flds): lonlat_(flds["lonlat"]), - ghost_(flds["ghost"]), remote_index_(flds["remote_index"]), partition_(flds["partition"]) { - setupHaloExchange(); +PointCloud::PointCloud(const FieldSet & flds): lonlat_(flds["lonlat"]) { + if (flds.has("ghost")) { + ghost_ = flds["ghost"]; + } + if (flds.has("remote_index")) { + remote_index_ = flds["remote_index"]; + } + if (flds.has("partition")) { + partition_ = flds["partition"]; + } + if (flds.has("global_index")) { + global_index_ = flds["global_index"]; + } + if( ghost_ && remote_index_ && partition_ ) { + setupHaloExchange(); + } } PointCloud::PointCloud(const Grid& grid) { diff --git a/src/atlas/functionspace/PointCloud.h b/src/atlas/functionspace/PointCloud.h index 376388cac..7c0e14ccd 100644 --- a/src/atlas/functionspace/PointCloud.h +++ b/src/atlas/functionspace/PointCloud.h @@ -56,6 +56,7 @@ class PointCloud : public functionspace::FunctionSpaceImpl { const Field& vertical() const { return vertical_; } Field ghost() const override; Field remote_index() const override { return remote_index_; } + Field global_index() const override { return global_index_; } virtual idx_t size() const override { return lonlat_.shape(0); } using FunctionSpaceImpl::createField; @@ -154,6 +155,7 @@ class PointCloud : public functionspace::FunctionSpaceImpl { Field vertical_; mutable Field ghost_; Field remote_index_; + Field global_index_; Field partition_; std::unique_ptr halo_exchange_; idx_t levels_{0}; diff --git a/src/atlas/interpolation/Cache.h b/src/atlas/interpolation/Cache.h index 3fdaebdfd..412957eb3 100644 --- a/src/atlas/interpolation/Cache.h +++ b/src/atlas/interpolation/Cache.h @@ -108,6 +108,7 @@ class MatrixCache final : public Cache { public: MatrixCache() = default; MatrixCache(const Cache& c); + MatrixCache(const MatrixCache& c) : MatrixCache(Cache(c)) {} MatrixCache(Matrix&& m); MatrixCache(std::shared_ptr m, const std::string& uid = ""); MatrixCache(const Matrix* m); @@ -146,6 +147,7 @@ class IndexKDTreeCache final : public Cache { public: IndexKDTreeCache() = default; IndexKDTreeCache(const Cache& c); + IndexKDTreeCache(const IndexKDTreeCache& c) : IndexKDTreeCache(Cache(c)) {} IndexKDTreeCache(const IndexKDTree&); IndexKDTreeCache(const Interpolation&); operator bool() const; diff --git a/src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc b/src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc index ce0e33b2b..9d4c3c589 100644 --- a/src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc +++ b/src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc @@ -105,18 +105,6 @@ namespace { return out.str(); }; - inline auto compute_src_west_east(functionspace::StructuredColumns src_fs) { - const auto& src_polygon = src_fs.polygon(0); - double src_west = std::numeric_limits::max(); - double src_east = std::numeric_limits::min(); - for( auto& p : src_polygon.lonlat() ) { - src_west = std::min(src_west, p[0]); - src_east = std::max(src_east, p[0]); - } - return std::make_tuple(src_west,src_east); - }; - - inline void output_source_partition_polygon(const std::string& path, functionspace::StructuredColumns src_fs, idx_t halo) { auto& polygon = src_fs.polygon(halo); std::ofstream f(filename(path)); @@ -182,13 +170,11 @@ namespace { output_source_partition_polygon("atlas_source_partition_polygons_halo_0_p%p.json",src_fs,0); output_source_partition_polygon("atlas_source_partition_polygons_halo_" + halo_str + "_p%p.json",src_fs,src_fs.halo()); output_target_points("atlas_target_failed_points_p%p.json", failed_points, lonlat); - auto [src_west, src_east] = compute_src_west_east(src_fs); idx_t my_rank = mpi::rank(); for( idx_t p=0; p::setup( const FunctionSpace& source ) { auto triplets = kernel_->allocate_triplets( out_npts_ ); using WorkSpace = typename Kernel::WorkSpace; - auto [west, east] = compute_src_west_east(source); - auto normalise = util::NormaliseLongitude(west); auto interpolate_point = [&]( idx_t n, PointLonLat&& p, WorkSpace& workspace ) -> int { try { - p.lon() = normalise(p.lon()); kernel_->insert_triplets( n, p, triplets, workspace ); return 0; } @@ -518,13 +501,11 @@ void StructuredInterpolation2D::execute_impl( const Kernel& kernel, cons using WorkSpace = typename Kernel::WorkSpace; - auto [west, east] = compute_src_west_east(source()); - util::NormaliseLongitude normalise{west}; auto interpolate_point = [&]( idx_t n, PointLonLat&& p, WorkSpace& workspace ) -> int { try { - p.lon() = normalise(p.lon()); kernel.compute_stencil( p.lon(), p.lat(), workspace.stencil ); kernel.compute_weights( p.lon(), p.lat(), workspace.stencil, workspace.weights ); + kernel.make_valid_stencil( p.lon(), p.lat(), workspace.stencil ); for ( idx_t i = 0; i < N; ++i ) { kernel.interpolate( workspace.stencil, workspace.weights, src_view[i], tgt_view[i], n ); } diff --git a/src/atlas/interpolation/method/structured/kernels/CubicHorizontalKernel.h b/src/atlas/interpolation/method/structured/kernels/CubicHorizontalKernel.h index b61b0cc51..2d5fa016b 100644 --- a/src/atlas/interpolation/method/structured/kernels/CubicHorizontalKernel.h +++ b/src/atlas/interpolation/method/structured/kernels/CubicHorizontalKernel.h @@ -72,15 +72,20 @@ class CubicHorizontalKernel { }; template - void compute_stencil(double& x, const double y, stencil_t& stencil, bool retry = true) const { + void compute_stencil(const double x, const double y, stencil_t& stencil) const { compute_horizontal_stencil_(x, y, stencil); + } + + template + void make_valid_stencil(double& x, const double y, stencil_t& stencil, bool retry = true) const { for (idx_t j = 0; j < stencil_width(); ++j) { idx_t imin = stencil.i(0, j); idx_t imax = stencil.i(stencil_width()-1, j); if (imin < src_.i_begin_halo(stencil.j(j))) { if (retry ) { x += 360.; - compute_stencil(x, y, stencil, false); + compute_stencil(x, y, stencil); + return make_valid_stencil(x, y, stencil, false); } else { Log::error() << "Stencil out of bounds" << std::endl; @@ -90,7 +95,8 @@ class CubicHorizontalKernel { if (imax >= src_.i_end_halo(stencil.j(j))) { if (retry ) { x -= 360.; - compute_stencil(x, y, stencil, false); + compute_stencil(x, y, stencil); + return make_valid_stencil(x, y, stencil, false); } else { Log::error() << "Stencil out of bounds" << std::endl; @@ -102,7 +108,7 @@ class CubicHorizontalKernel { } template - void compute_weights(double x, const double y, weights_t& weights) const { + void compute_weights(const double x, const double y, weights_t& weights) const { Stencil stencil; compute_stencil(x, y, stencil); compute_weights(x, y, stencil, weights); @@ -234,6 +240,7 @@ class CubicHorizontalKernel { compute_stencil(x, y, stencil); Weights weights; compute_weights(x, y, stencil, weights); + make_valid_stencil(x, y, stencil); return interpolate(stencil, weights, input); } @@ -241,6 +248,7 @@ class CubicHorizontalKernel { typename array_t::value_type interpolate(PointXY p, const array_t& input, WorkSpace& ws) const { compute_stencil(p.x(), p.y(), ws.stencil); compute_weights(p.x(), p.y(), ws.stencil, ws.weights); + make_valid_stencil(p.x(), p.y(), ws.stencil); return interpolate(ws.stencil, ws.weights, input); } @@ -267,6 +275,8 @@ class CubicHorizontalKernel { void insert_triplets(const idx_t row, double x, double y, Triplets& triplets, WorkSpace& ws) const { compute_stencil(x, y, ws.stencil); compute_weights(x, y, ws.stencil, ws.weights); + + make_valid_stencil(x, y, ws.stencil); const auto& wj = ws.weights.weights_j; idx_t pos = row * stencil_size(); diff --git a/src/atlas/interpolation/method/structured/kernels/LinearHorizontalKernel.h b/src/atlas/interpolation/method/structured/kernels/LinearHorizontalKernel.h index cccab6d0c..d0815c7c2 100644 --- a/src/atlas/interpolation/method/structured/kernels/LinearHorizontalKernel.h +++ b/src/atlas/interpolation/method/structured/kernels/LinearHorizontalKernel.h @@ -64,15 +64,20 @@ class LinearHorizontalKernel { }; template - void compute_stencil(double& x, double y, stencil_t& stencil, bool retry = true) const { + void compute_stencil(const double x, const double y, stencil_t& stencil) const { compute_horizontal_stencil_(x, y, stencil); + } + + template + void make_valid_stencil(double& x, const double y, stencil_t& stencil, bool retry = true) const { for (idx_t j = 0; j < stencil_width(); ++j) { idx_t imin = stencil.i(0, j); idx_t imax = stencil.i(stencil_width()-1, j); if (imin < src_.i_begin_halo(stencil.j(j))) { if (retry ) { x += 360.; - compute_stencil(x, y, stencil, false); + compute_stencil(x, y, stencil); + return make_valid_stencil(x, y, stencil, false); } else { Log::error() << "Stencil out of bounds" << std::endl; @@ -82,7 +87,8 @@ class LinearHorizontalKernel { if (imax >= src_.i_end_halo(stencil.j(j))) { if (retry ) { x -= 360.; - compute_stencil(x, y, stencil, false); + compute_stencil(x, y, stencil); + return make_valid_stencil(x, y, stencil, false); } else { Log::error() << "Stencil out of bounds" << std::endl; @@ -94,7 +100,7 @@ class LinearHorizontalKernel { template - void compute_weights(double x, double y, weights_t& weights) const { + void compute_weights(const double x, const double y, weights_t& weights) const { Stencil stencil; compute_stencil(x, y, stencil); compute_weights(x, y, stencil, weights); @@ -186,6 +192,7 @@ class LinearHorizontalKernel { Weights weights; compute_stencil(x, y, stencil); compute_weights(x, y, stencil, weights); + make_valid_stencil(x, y, stencil); return interpolate(stencil, weights, input); } @@ -193,6 +200,7 @@ class LinearHorizontalKernel { typename array_t::value_type interpolate(PointXY p, const array_t& input, WorkSpace& ws) const { compute_stencil(p.x(), p.y(), ws.stencil); compute_weights(p.x(), p.y(), ws.stencil, ws.weights); + make_valid_stencil(p.x(), p.y(), ws.stencil); return interpolate(ws.stencil, ws.weights, input); } @@ -217,8 +225,10 @@ class LinearHorizontalKernel { } void insert_triplets(const idx_t row, double x, double y, Triplets& triplets, WorkSpace& ws) const { - compute_stencil(x, y, ws.stencil, true); + compute_stencil(x, y, ws.stencil); compute_weights(x, y, ws.stencil, ws.weights); + + make_valid_stencil(x, y, ws.stencil); const auto& wj = ws.weights.weights_j; idx_t pos = row * stencil_size(); diff --git a/src/atlas/interpolation/method/structured/kernels/QuasiCubicHorizontalKernel.h b/src/atlas/interpolation/method/structured/kernels/QuasiCubicHorizontalKernel.h index e9280bd76..463cb9765 100644 --- a/src/atlas/interpolation/method/structured/kernels/QuasiCubicHorizontalKernel.h +++ b/src/atlas/interpolation/method/structured/kernels/QuasiCubicHorizontalKernel.h @@ -74,15 +74,20 @@ class QuasiCubicHorizontalKernel { }; template - void compute_stencil(double& x, const double y, stencil_t& stencil, bool retry = true) const { + void compute_stencil(const double x, const double y, stencil_t& stencil) const { compute_horizontal_stencil_(x, y, stencil); + } + + template + void make_valid_stencil(double& x, double y, stencil_t& stencil, bool retry = true) const { for (idx_t j = 0; j < stencil_width(); ++j) { idx_t imin = stencil.i(0, j); idx_t imax = stencil.i(stencil_width()-1, j); if (imin < src_.i_begin_halo(stencil.j(j))) { if (retry ) { x += 360.; - compute_stencil(x, y, stencil, false); + compute_stencil(x, y, stencil); + return make_valid_stencil(x, y, stencil, false); } else { Log::error() << "Stencil out of bounds" << std::endl; @@ -92,7 +97,8 @@ class QuasiCubicHorizontalKernel { if (imax >= src_.i_end_halo(stencil.j(j))) { if (retry ) { x -= 360.; - compute_stencil(x, y, stencil, false); + compute_stencil(x, y, stencil); + return make_valid_stencil(x, y, stencil, false); } else { Log::error() << "Stencil out of bounds" << std::endl; @@ -285,6 +291,7 @@ class QuasiCubicHorizontalKernel { compute_stencil(x, y, stencil); Weights weights; compute_weights(x, y, stencil, weights); + make_valid_stencil(x, y, stencil); return interpolate(stencil, weights, input); } @@ -292,6 +299,7 @@ class QuasiCubicHorizontalKernel { typename array_t::value_type interpolate(PointXY p, const array_t& input, WorkSpace& ws) const { compute_stencil(p.x(), p.y(), ws.stencil); compute_weights(p.x(), p.y(), ws.stencil, ws.weights); + make_valid_stencil(p.x(), p.y(), ws.stencil); return interpolate(ws.stencil, ws.weights, input); } @@ -318,6 +326,8 @@ class QuasiCubicHorizontalKernel { void insert_triplets(const idx_t row, double x, double y, Triplets& triplets, WorkSpace& ws) const { compute_stencil(x, y, ws.stencil); compute_weights(x, y, ws.stencil, ws.weights); + + make_valid_stencil(x, y, ws.stencil); const auto& wj = ws.weights.weights_j; idx_t pos = row * stencil_size(); 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/library/detail/Debug.h b/src/atlas/library/detail/Debug.h index 1e2d4f9d7..288b14288 100644 --- a/src/atlas/library/detail/Debug.h +++ b/src/atlas/library/detail/Debug.h @@ -20,6 +20,12 @@ namespace atlas { namespace debug { +inline gidx_t global_index(int i = 0) { + static std::vector g = + eckit::Resource>("$ATLAS_DEBUG_GLOBAL_INDEX", std::vector{-1}); + return g[i]; +} + inline gidx_t node_global_index(int i = 0) { static std::vector g = eckit::Resource>("$ATLAS_DEBUG_NODE_GLOBAL_INDEX", std::vector{-1}); @@ -54,6 +60,16 @@ inline bool is_node_global_index(gidx_t x) { return false; } +inline bool is_global_index(gidx_t x) { + static std::vector v = + eckit::Resource>("$ATLAS_GLOBAL_INDEX", std::vector()); + for (gidx_t g : v) { + if (x == g) + return true; + } + return false; +} + inline bool is_edge_global_index(gidx_t x) { static std::vector v = eckit::Resource>("$ATLAS_DEBUG_EDGE_GLOBAL_INDEX", std::vector()); diff --git a/src/atlas/mesh/ElementType.cc b/src/atlas/mesh/ElementType.cc index f301f09d0..60c033d3d 100644 --- a/src/atlas/mesh/ElementType.cc +++ b/src/atlas/mesh/ElementType.cc @@ -9,6 +9,10 @@ */ #include "atlas/mesh/ElementType.h" +#include "elementtypes/Line.h" +#include "elementtypes/Triangle.h" +#include "elementtypes/Quadrilateral.h" +#include "elementtypes/Pentagon.h" #include "atlas/runtime/Exception.h" #include "atlas/util/CoordinateEnums.h" @@ -17,8 +21,22 @@ namespace mesh { //------------------------------------------------------------------------------ -ElementType* ElementType::create(const std::string&) { - ATLAS_NOTIMPLEMENTED; +ElementType* ElementType::create( const std::string& name ) +{ + // currently naive implementation. This has to be replaced + // with self-registration and factory mechanism + if( name == "Triangle" ) + return new elementtypes::Triangle(); + else if( name == "Quadrilateral") + return new elementtypes::Quadrilateral(); + else if( name == "Line" ) + return new elementtypes::Line(); + else if( name == "Pentagon") + return new elementtypes::Pentagon(); + else { + throw_NotImplemented("Element type ["+name+"] does not exist", Here()); + } + return 0; } ElementType::ElementType() = default; @@ -30,14 +48,10 @@ extern "C" { void atlas__mesh__ElementType__delete(ElementType* This) { delete This; } -ElementType* atlas__mesh__Triangle__create() { - return new temporary::Triangle(); -} -ElementType* atlas__mesh__Quadrilateral__create() { - return new temporary::Quadrilateral(); -} -ElementType* atlas__mesh__Line__create() { - return new temporary::Line(); + +ElementType* atlas__mesh__ElementType__create(const char* name) +{ + return ElementType::create( std::string(name) ); } const char* atlas__mesh__ElementType__name(const ElementType* This) { diff --git a/src/atlas/mesh/ElementType.h b/src/atlas/mesh/ElementType.h index c96464d3c..3d3d6a06c 100644 --- a/src/atlas/mesh/ElementType.h +++ b/src/atlas/mesh/ElementType.h @@ -33,212 +33,30 @@ class ElementType : public util::Object { //-- Constructors ElementType(); - ~ElementType() = 0; - //-- Accessors + ~ElementType() = 0; virtual const std::string& name() const = 0; - // virtual idx_t dimensionality() const = 0; - - // virtual idx_t nb_vertices() const = 0; - virtual idx_t nb_edges() const = 0; - // virtual idx_t nb_faces() const = 0; - virtual idx_t nb_nodes() const = 0; - - virtual bool parametric() const = 0; -}; + virtual size_t dimensionality() const = 0; -namespace temporary { + virtual size_t nb_vertices() const = 0; -class Volume : public ElementType { -public: - enum - { - DIMENSIONALITY = 3 - }; -}; + virtual size_t nb_edges() const = 0; -class Face : public ElementType { -public: - enum - { - DIMENSIONALITY = 2 - }; - enum - { - FACES = 1 - }; - virtual idx_t nb_faces() const { return FACES; } -}; + virtual size_t nb_faces() const = 0; -class Edge : public ElementType { -public: - enum - { - DIMENSIONALITY = 1 - }; - enum - { - FACES = 0 - }; - enum - { - EDGES = 1 - }; - virtual idx_t nb_faces() const { return FACES; } - virtual idx_t nb_edges() const { return EDGES; } -}; + virtual size_t nb_nodes() const = 0; -class Vertex : public ElementType { -public: - enum - { - DIMENSIONALITY = 0 - }; - enum - { - FACES = 0 - }; - enum - { - EDGES = 0 - }; - enum - { - VERTICES = 1 - }; - virtual idx_t nb_faces() const { return FACES; } - virtual idx_t nb_edges() const { return EDGES; } - virtual idx_t nb_vertices() const { return VERTICES; } -}; + virtual bool parametric() const = 0; -class Quadrilateral : public Face { -public: - enum - { - EDGES = 4 - }; - enum - { - VERTICES = 4 - }; - enum - { - FACETS = EDGES - }; - enum - { - RIDGES = VERTICES - }; - virtual ~Quadrilateral() {} - virtual bool parametric() const { return true; } - virtual idx_t nb_vertices() const { return VERTICES; } - virtual idx_t nb_edges() const { return EDGES; } - virtual idx_t nb_nodes() const { return VERTICES; } - virtual idx_t nb_facets() const { return FACETS; } - virtual idx_t nb_ridges() const { return RIDGES; } - virtual const std::string& name() const { - static std::string s("Quadrilateral"); - return s; - } -}; + virtual bool simplex() const = 0; -class Triangle : public Face { -public: - enum - { - EDGES = 3 - }; - enum - { - VERTICES = 3 - }; - enum - { - FACETS = EDGES - }; - enum - { - RIDGES = VERTICES - }; - virtual ~Triangle() {} - virtual bool parametric() const { return true; } - virtual idx_t nb_vertices() const { return VERTICES; } - virtual idx_t nb_edges() const { return EDGES; } - virtual idx_t nb_nodes() const { return VERTICES; } - virtual idx_t nb_facets() const { return FACETS; } - virtual idx_t nb_ridges() const { return RIDGES; } - virtual const std::string& name() const { - static std::string s("Triangle"); - return s; - } }; -class Line : public Edge { -public: - enum - { - VERTICES = 2 - }; - enum - { - FACETS = VERTICES - }; - enum - { - RIDGES = 0 - }; - virtual ~Line() {} - virtual bool parametric() const { return true; } - virtual idx_t nb_vertices() const { return VERTICES; } - virtual idx_t nb_edges() const { return EDGES; } - virtual idx_t nb_nodes() const { return VERTICES; } - virtual idx_t nb_facets() const { return FACETS; } - virtual const std::string& name() const { - static std::string s("Line"); - return s; - } -}; - -class Pentagon : public Face { -public: - enum - { - EDGES = 5 - }; - enum - { - VERTICES = 5 - }; - enum - { - FACETS = EDGES - }; - enum - { - RIDGES = VERTICES - }; - virtual ~Pentagon() {} - virtual bool parametric() const { return false; } - virtual idx_t nb_vertices() const { return VERTICES; } - virtual idx_t nb_edges() const { return EDGES; } - virtual idx_t nb_nodes() const { return VERTICES; } - virtual idx_t nb_facets() const { return FACETS; } - virtual idx_t nb_ridges() const { return RIDGES; } - virtual const std::string& name() const { - static std::string s("Pentagon"); - return s; - } -}; - -} // namespace temporary - -extern "C" { -ElementType* atlas__mesh__Triangle__create(); -ElementType* atlas__mesh__Quadrilateral__create(); -ElementType* atlas__mesh__Line__create(); - +extern "C" +{ +ElementType* atlas__mesh__ElementType__create(const char* name); void atlas__mesh__ElementType__delete(ElementType* This); idx_t atlas__mesh__ElementType__nb_nodes(const ElementType* This); idx_t atlas__mesh__ElementType__nb_edges(const ElementType* This); @@ -250,3 +68,10 @@ const char* atlas__mesh__ElementType__name(const ElementType* This); } // namespace mesh } // namespace atlas + +// Following includes are added for backwards compatibility with temporary elementtypes. +// They contain deprecation warnings. +// When deprecation is final, these lines can be removed. +#include "atlas/mesh/elementtypes/Quadrilateral.h" +#include "atlas/mesh/elementtypes/Triangle.h" +#include "atlas/mesh/elementtypes/Line.h" diff --git a/src/atlas/mesh/MeshBuilder.cc b/src/atlas/mesh/MeshBuilder.cc index 7d1cf2f15..ca35b0d4b 100644 --- a/src/atlas/mesh/MeshBuilder.cc +++ b/src/atlas/mesh/MeshBuilder.cc @@ -223,10 +223,10 @@ Mesh MeshBuilder::operator()(size_t nb_nodes, const double lons[], const double const bool add_quads = (sum_nb_quads > 0); if (add_tris) { - mesh.cells().add(new mesh::temporary::Triangle(), nb_tris); + mesh.cells().add(mesh::ElementType::create("Triangle"), nb_tris); } if (add_quads) { - mesh.cells().add(new mesh::temporary::Quadrilateral(), nb_quads); + mesh.cells().add(mesh::ElementType::create("Quadrilateral"), nb_quads); } atlas::mesh::HybridElements::Connectivity& node_connectivity = mesh.cells().node_connectivity(); diff --git a/src/atlas/mesh/actions/BuildConvexHull3D.cc b/src/atlas/mesh/actions/BuildConvexHull3D.cc index 1aeb0cf9e..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(new mesh::temporary::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/mesh/actions/BuildEdges.cc b/src/atlas/mesh/actions/BuildEdges.cc index ac773bb01..1c7ccda1e 100644 --- a/src/atlas/mesh/actions/BuildEdges.cc +++ b/src/atlas/mesh/actions/BuildEdges.cc @@ -415,7 +415,7 @@ void build_edges(Mesh& mesh, const eckit::Configuration& config) { } // Build edges - mesh.edges().add(new mesh::temporary::Line(), (edge_end - edge_start), + mesh.edges().add(mesh::ElementType::create("Line"), (edge_end - edge_start), edge_nodes_data.data() + edge_halo_offsets[halo] * 2); auto& edge_nodes = mesh.edges().node_connectivity(); const auto& cell_nodes = mesh.cells().node_connectivity(); @@ -474,7 +474,7 @@ void build_edges(Mesh& mesh, const eckit::Configuration& config) { edge_start = edge_end; edge_end += nb_pole_edges; - mesh.edges().add(new mesh::temporary::Line(), nb_pole_edges, pole_edge_nodes.data()); + mesh.edges().add(mesh::ElementType::create("Line"), nb_pole_edges, pole_edge_nodes.data()); auto edge_ridx = array::make_indexview(mesh.edges().remote_index()); auto edge_part = array::make_view(mesh.edges().partition()); diff --git a/src/atlas/mesh/elementtypes/Classification.h b/src/atlas/mesh/elementtypes/Classification.h new file mode 100644 index 000000000..48e951b4b --- /dev/null +++ b/src/atlas/mesh/elementtypes/Classification.h @@ -0,0 +1,74 @@ +/* + * (C) Copyright 1996-2016 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. + */ + +/// @author Willem Deconinck +/// @date October 2015 + +#pragma once + +#include "atlas/mesh/ElementType.h" + +namespace atlas { +namespace mesh { +namespace elementtypes { + +//------------------------------------------------------------------- + +class Volume : public ElementType +{ +public: + enum { DIMENSIONALITY=3 }; + virtual size_t dimensionality() const { return DIMENSIONALITY; } +}; + +//------------------------------------------------------------------- + +class Face : public ElementType +{ +public: + enum { DIMENSIONALITY=2 }; + enum { FACES=1 }; + virtual size_t dimensionality() const { return DIMENSIONALITY; } + virtual size_t nb_faces() const { return FACES; } +}; + +//------------------------------------------------------------------- + +class Edge: public ElementType +{ +public: + enum { DIMENSIONALITY=1 }; + enum { FACES=0 }; + enum { EDGES=1 }; + virtual size_t dimensionality() const { return DIMENSIONALITY; } + virtual size_t nb_faces() const { return FACES; } + virtual size_t nb_edges() const { return EDGES; } +}; + +//------------------------------------------------------------------- + +class Vertex: public ElementType +{ +public: + enum { DIMENSIONALITY=0 }; + enum { FACES=0 }; + enum { EDGES=0 }; + enum { VERTICES=1 }; + virtual size_t dimensionality() const { return DIMENSIONALITY; } + virtual size_t nb_faces() const { return FACES; } + virtual size_t nb_edges() const { return EDGES; } + virtual size_t nb_vertices() const { return VERTICES; } +}; + +//------------------------------------------------------------------- + +} // namespace elementtypes +} // namespace mesh +} // namespace atlas diff --git a/src/atlas/mesh/elementtypes/Line.h b/src/atlas/mesh/elementtypes/Line.h new file mode 100644 index 000000000..8752b33cd --- /dev/null +++ b/src/atlas/mesh/elementtypes/Line.h @@ -0,0 +1,52 @@ +/* + * (C) Copyright 1996-2016 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. + */ + +/// @author Willem Deconinck +/// @date October 2015 + +#pragma once + +#include "atlas/mesh/elementtypes/Classification.h" + +namespace atlas { +namespace mesh { +namespace elementtypes { + +//------------------------------------------------------------------------------------------------------ + +class Line : public Edge +{ +public: + enum { VERTICES = 2 }; + enum { FACETS = VERTICES }; + enum { RIDGES = 0 }; + virtual ~Line() {} + virtual bool parametric() const { return true; } + virtual bool simplex() const { return false; } + virtual size_t nb_vertices() const { return VERTICES; } + virtual size_t nb_edges() const { return EDGES; } + virtual size_t nb_nodes() const { return VERTICES; } + virtual size_t nb_facets() const { return FACETS; } + virtual const std::string& name() const { static std::string s("Line"); return s; } +}; + +//------------------------------------------------------------------------------------------------------ + +} // namespace elementtypes + +namespace temporary { + class Line : public elementtypes::Line { + public: + [[deprecated("Use 'atlas::mesh::ElementType::create(\"Line\")' instead")]] Line() {} + }; +} + +} // namespace mesh +} // namespace atlas diff --git a/src/atlas/mesh/elementtypes/Pentagon.h b/src/atlas/mesh/elementtypes/Pentagon.h new file mode 100644 index 000000000..41199accf --- /dev/null +++ b/src/atlas/mesh/elementtypes/Pentagon.h @@ -0,0 +1,61 @@ +/* + * (C) Copyright 1996-2016 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. + */ + +/// @author Willem Deconinck +/// @date October 2015 + +#pragma once + +#include "atlas/mesh/elementtypes/Classification.h" + +namespace atlas { +namespace mesh { +namespace elementtypes { + +//------------------------------------------------------------------------------------------------------ + +class Pentagon : public Face { +public: + enum + { + EDGES = 5 + }; + enum + { + VERTICES = 5 + }; + enum + { + FACETS = EDGES + }; + enum + { + RIDGES = VERTICES + }; + virtual ~Pentagon() {} + virtual bool parametric() const { return false; } + virtual bool simplex() const { return false; } + virtual size_t nb_vertices() const { return VERTICES; } + virtual size_t nb_edges() const { return EDGES; } + virtual size_t nb_nodes() const { return VERTICES; } + virtual size_t nb_facets() const { return FACETS; } + virtual size_t nb_ridges() const { return RIDGES; } + virtual const std::string& name() const { + static std::string s("Pentagon"); + return s; + } +}; + +//------------------------------------------------------------------------------------------------------ + +} // namespace elementtypes +} // namespace mesh +} // namespace atlas + diff --git a/src/atlas/mesh/elementtypes/Quadrilateral.h b/src/atlas/mesh/elementtypes/Quadrilateral.h new file mode 100644 index 000000000..42fe12eac --- /dev/null +++ b/src/atlas/mesh/elementtypes/Quadrilateral.h @@ -0,0 +1,55 @@ +/* + * (C) Copyright 1996-2016 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. + */ + +/// @author Willem Deconinck +/// @date October 2015 + +#pragma once + +#include "atlas/mesh/elementtypes/Classification.h" + +namespace atlas { +namespace mesh { +namespace elementtypes { + +//------------------------------------------------------------------------------------------------------ + +class Quadrilateral : public Face +{ +public: + enum { EDGES=4 }; + enum { VERTICES=4 }; + enum { FACETS=EDGES }; + enum { RIDGES=VERTICES }; + virtual ~Quadrilateral() {} + virtual bool parametric() const { return true; } + virtual bool simplex() const { return false; } + virtual size_t nb_vertices() const { return VERTICES; } + virtual size_t nb_edges() const { return EDGES; } + virtual size_t nb_nodes() const { return VERTICES; } + virtual size_t nb_facets() const { return FACETS; } + virtual size_t nb_ridges() const { return RIDGES; } + virtual const std::string& name() const { static std::string s("Quadrilateral"); return s; } +}; + +//------------------------------------------------------------------------------------------------------ + +} // namespace elementtypes + +namespace temporary { + class Quadrilateral : public elementtypes::Quadrilateral { + public: + [[deprecated("Use 'atlas::mesh::ElementType::create(\"Quadrilateral\")' instead")]] Quadrilateral() {} + }; +} + +} // namespace mesh +} // namespace atlas + diff --git a/src/atlas/mesh/elementtypes/Triangle.h b/src/atlas/mesh/elementtypes/Triangle.h new file mode 100644 index 000000000..04cea1c2d --- /dev/null +++ b/src/atlas/mesh/elementtypes/Triangle.h @@ -0,0 +1,54 @@ +/* + * (C) Copyright 1996-2016 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. + */ + +/// @author Willem Deconinck +/// @date October 2015 + +#pragma once + +#include "atlas/mesh/elementtypes/Classification.h" + +namespace atlas { +namespace mesh { +namespace elementtypes { + +//------------------------------------------------------------------------------------------------------ + +class Triangle : public Face +{ +public: + enum { EDGES=3 }; + enum { VERTICES=3 }; + enum { FACETS=EDGES }; + enum { RIDGES=VERTICES }; + virtual ~Triangle() {} + virtual bool parametric() const { return true; } + virtual bool simplex() const { return true; } + virtual size_t nb_vertices() const { return VERTICES; } + virtual size_t nb_edges() const { return EDGES; } + virtual size_t nb_nodes() const { return VERTICES; } + virtual size_t nb_facets() const { return FACETS; } + virtual size_t nb_ridges() const { return RIDGES; } + virtual const std::string& name() const { static std::string s("Triangle"); return s; } +}; + +//------------------------------------------------------------------------------------------------------ + +} // namespace elementtypes + +namespace temporary { + class Triangle : public elementtypes::Triangle { + public: + [[deprecated("Use 'atlas::mesh::ElementType::create(\"Triangle\")' instead")]] Triangle() {} + }; +} + +} // namespace mesh +} // namespace atlas diff --git a/src/atlas/meshgenerator/detail/CubedSphereDualMeshGenerator.cc b/src/atlas/meshgenerator/detail/CubedSphereDualMeshGenerator.cc index 235b94e4e..66beb592e 100644 --- a/src/atlas/meshgenerator/detail/CubedSphereDualMeshGenerator.cc +++ b/src/atlas/meshgenerator/detail/CubedSphereDualMeshGenerator.cc @@ -388,17 +388,17 @@ void CubedSphereDualMeshGenerator::generate_mesh(const CubedSphereGrid& csGrid, switch (typeCount.first) { // Add a block of lines. case ElemType::LINE: { - cells.add(new mesh::temporary::Line(), typeCount.second); + cells.add(mesh::ElementType::create("Line"), typeCount.second); break; } // Add a block of triangles. case ElemType::TRIANGLE: { - cells.add(new mesh::temporary::Triangle(), typeCount.second); + cells.add(mesh::ElementType::create("Triangle"), typeCount.second); break; } // Add a block of quadrilaterals. case ElemType::QUADRILATERAL: { - cells.add(new mesh::temporary::Quadrilateral(), typeCount.second); + cells.add(mesh::ElementType::create("Quadrilateral"), typeCount.second); break; } default: { diff --git a/src/atlas/meshgenerator/detail/CubedSphereMeshGenerator.cc b/src/atlas/meshgenerator/detail/CubedSphereMeshGenerator.cc index bae86c1bd..f719dc461 100644 --- a/src/atlas/meshgenerator/detail/CubedSphereMeshGenerator.cc +++ b/src/atlas/meshgenerator/detail/CubedSphereMeshGenerator.cc @@ -818,7 +818,7 @@ void CubedSphereMeshGenerator::generate_mesh(const CubedSphereGrid& csGrid, cons auto& cells = mesh.cells(); // Resize cells. - cells.add(new mesh::temporary::Quadrilateral(), static_cast(localCells.size())); + cells.add(mesh::ElementType::create("Quadrilateral"), static_cast(localCells.size())); // Add extra fields. tijField = cells.add(Field("tij", make_datatype(), make_shape(cells.size(), 3))); diff --git a/src/atlas/meshgenerator/detail/DelaunayMeshGenerator.cc b/src/atlas/meshgenerator/detail/DelaunayMeshGenerator.cc index 641c07694..f8a3b0a2e 100644 --- a/src/atlas/meshgenerator/detail/DelaunayMeshGenerator.cc +++ b/src/atlas/meshgenerator/detail/DelaunayMeshGenerator.cc @@ -336,7 +336,7 @@ void DelaunayMeshGenerator::generate(const Grid& grid, const grid::Distribution& } - mesh.cells().add(new mesh::temporary::Triangle(), owned_elements.size()); + mesh.cells().add(mesh::ElementType::create("Triangle"), owned_elements.size()); auto& node_connectivity = mesh.cells().node_connectivity(); auto cell_gidx = array::make_view(mesh.cells().global_index()); auto cell_part = array::make_view(mesh.cells().partition()); diff --git a/src/atlas/meshgenerator/detail/HealpixMeshGenerator.cc b/src/atlas/meshgenerator/detail/HealpixMeshGenerator.cc index 189787e96..d77af3652 100644 --- a/src/atlas/meshgenerator/detail/HealpixMeshGenerator.cc +++ b/src/atlas/meshgenerator/detail/HealpixMeshGenerator.cc @@ -525,8 +525,8 @@ void HealpixMeshGenerator::generate_mesh(const StructuredGrid& grid, const grid: // array::ArrayView ghost ( nodes.ghost() ); // array::ArrayView flags ( nodes.flags() ); // - define cells (only quadrilaterals for now) with - // mesh.cells().add( new mesh::temporary::Quadrilateral(), nquads ); - // mesh.cells().add( new mesh::temporary::Pentagon(), npents ); + // mesh.cells().add( mesh::ElementType::create("Quadrilateral"), nquads ); + // mesh.cells().add( mesh::ElementType::create("Pentagon"), npents ); // further define cells with // array::ArrayView cells_glb_idx( mesh.cells().global_index() // ); @@ -803,11 +803,17 @@ void HealpixMeshGenerator::generate_mesh(const StructuredGrid& grid, const grid: auto halo = array::make_view(nodes.halo()); auto flags = array::make_view(nodes.flags()); + int quad_begin; + int pent_begin; + // define cells and associated properties - mesh.cells().add(new mesh::temporary::Quadrilateral(), nquads); - mesh.cells().add(new mesh::temporary::Pentagon(), npents); - int quad_begin = mesh.cells().elements(0).begin(); - int pent_begin = mesh.cells().elements(1).begin(); + mesh.cells().add(mesh::ElementType::create("Quadrilateral"), nquads); + quad_begin = mesh.cells().elements(0).begin(); + pent_begin = quad_begin + mesh.cells().elements(0).size(); + if (pole_elements == "pentagons") { + mesh.cells().add(mesh::ElementType::create("Pentagon"), npents); + pent_begin = mesh.cells().elements(1).begin(); + } auto cells_part = array::make_view(mesh.cells().partition()); auto cells_glb_idx = array::make_view(mesh.cells().global_index()); auto& node_connectivity = mesh.cells().node_connectivity(); diff --git a/src/atlas/meshgenerator/detail/NodalCubedSphereMeshGenerator.cc b/src/atlas/meshgenerator/detail/NodalCubedSphereMeshGenerator.cc index 2ed0290d2..913e9ea37 100644 --- a/src/atlas/meshgenerator/detail/NodalCubedSphereMeshGenerator.cc +++ b/src/atlas/meshgenerator/detail/NodalCubedSphereMeshGenerator.cc @@ -327,7 +327,7 @@ void NodalCubedSphereMeshGenerator::generate(const Grid& grid, const grid::Distr } // Cells in mesh - mesh.cells().add(new mesh::temporary::Quadrilateral(), nTiles * N * N); + mesh.cells().add(mesh::ElementType::create("Quadrilateral"), nTiles * N * N); //int quad_begin = mesh.cells().elements( 0 ).begin(); auto cells_part = array::make_view(mesh.cells().partition()); auto cells_gidx = array::make_view(mesh.cells().global_index()); diff --git a/src/atlas/meshgenerator/detail/RegularMeshGenerator.cc b/src/atlas/meshgenerator/detail/RegularMeshGenerator.cc index 07f36849b..6035bc64a 100644 --- a/src/atlas/meshgenerator/detail/RegularMeshGenerator.cc +++ b/src/atlas/meshgenerator/detail/RegularMeshGenerator.cc @@ -189,7 +189,7 @@ void RegularMeshGenerator::generate_mesh(const RegularGrid& rg, const grid::Dist // array::ArrayView ghost ( nodes.ghost() ); // array::ArrayView flags ( nodes.flags() ); // - define cells (only quadrilaterals for now) with - // mesh.cells().add( new mesh::temporary::Quadrilateral(), nquads ); + // mesh.cells().add(mesh::ElementType::create("Quadrilateral"), nquads ); // further define cells with // array::ArrayView cells_glb_idx( mesh.cells().global_index() // ); @@ -396,7 +396,7 @@ void RegularMeshGenerator::generate_mesh(const RegularGrid& rg, const grid::Dist auto flags = array::make_view(nodes.flags()); // define cells and associated properties - mesh.cells().add(new mesh::temporary::Quadrilateral(), ncells); + mesh.cells().add(mesh::ElementType::create("Quadrilateral"), ncells); int quad_begin = mesh.cells().elements(0).begin(); auto cells_part = array::make_view(mesh.cells().partition()); mesh::HybridElements::Connectivity& node_connectivity = mesh.cells().node_connectivity(); diff --git a/src/atlas/meshgenerator/detail/StructuredMeshGenerator.cc b/src/atlas/meshgenerator/detail/StructuredMeshGenerator.cc index e8d7a5b0b..0c57e83e5 100644 --- a/src/atlas/meshgenerator/detail/StructuredMeshGenerator.cc +++ b/src/atlas/meshgenerator/detail/StructuredMeshGenerator.cc @@ -1248,8 +1248,8 @@ void StructuredMeshGenerator::generate_mesh(const StructuredGrid& rg, const grid // Now handle elements // ------------------- - mesh.cells().add(new mesh::temporary::Quadrilateral(), nquads); - mesh.cells().add(new mesh::temporary::Triangle(), ntriags); + mesh.cells().add(mesh::ElementType::create("Quadrilateral"), nquads); + mesh.cells().add(mesh::ElementType::create("Triangle"), ntriags); mesh::HybridElements::Connectivity& node_connectivity = mesh.cells().node_connectivity(); array::ArrayView cells_glb_idx = array::make_view(mesh.cells().global_index()); diff --git a/src/atlas/output/detail/GmshIO.cc b/src/atlas/output/detail/GmshIO.cc index 0540ced8f..b69f6f5e5 100644 --- a/src/atlas/output/detail/GmshIO.cc +++ b/src/atlas/output/detail/GmshIO.cc @@ -661,13 +661,13 @@ Mesh GmshIO::read(const PathName& file_path) const { namespace { mesh::ElementType* make_element_type(int type) { if (type == QUAD) { - return new mesh::temporary::Quadrilateral(); + return mesh::ElementType::create("Quadrilateral"); } if (type == TRIAG) { - return new mesh::temporary::Triangle(); + return mesh::ElementType::create("Triangle"); } if (type == LINE) { - return new mesh::temporary::Line(); + return mesh::ElementType::create("Line"); } throw_Exception("Element type not supported", Here()); } 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/src/atlas_f/atlas_module.F90 b/src/atlas_f/atlas_module.F90 index daff7976e..4e184c8a6 100644 --- a/src/atlas_f/atlas_module.F90 +++ b/src/atlas_f/atlas_module.F90 @@ -121,7 +121,7 @@ module atlas_module & atlas_MatchingMeshPartitioner ! Deprecated ! use atlas_MatchingPartitioner instead use atlas_MeshGenerator_module, only: & & atlas_MeshGenerator -#ifdef ATLAS_HAVE_NUMERICS +#if ATLAS_HAVE_NUMERICS use atlas_Method_module, only: & & atlas_Method use atlas_fvm_module, only: & diff --git a/src/atlas_f/mesh/atlas_ElementType_module.F90 b/src/atlas_f/mesh/atlas_ElementType_module.F90 index 98dac057a..a3187212b 100644 --- a/src/atlas_f/mesh/atlas_ElementType_module.F90 +++ b/src/atlas_f/mesh/atlas_ElementType_module.F90 @@ -71,22 +71,25 @@ function atlas_ElementType__cptr(cptr) result(this) function atlas_Quadrilateral__constructor() result(this) use atlas_elementtype_c_binding + use fckit_c_interop_module, only : c_str type(atlas_ElementType) :: this - call this%reset_c_ptr( atlas__mesh__Quadrilateral__create() ) + call this%reset_c_ptr( atlas__mesh__ElementType__create(c_str("Quadrilateral")) ) call this%return() end function function atlas_Triangle__constructor() result(this) use atlas_elementtype_c_binding + use fckit_c_interop_module, only : c_str type(atlas_ElementType) :: this - call this%reset_c_ptr( atlas__mesh__Triangle__create() ) + call this%reset_c_ptr( atlas__mesh__ElementType__create(c_str("Triangle")) ) call this%return() end function function atlas_Line__constructor() result(this) use atlas_elementtype_c_binding + use fckit_c_interop_module, only : c_str type(atlas_ElementType) :: this - call this%reset_c_ptr( atlas__mesh__Line__create() ) + call this%reset_c_ptr( atlas__mesh__ElementType__create(c_str("Line")) ) call this%return() end function diff --git a/src/tests/mesh/test_elements.cc b/src/tests/mesh/test_elements.cc index aa29519ad..f31b82cfc 100644 --- a/src/tests/mesh/test_elements.cc +++ b/src/tests/mesh/test_elements.cc @@ -26,7 +26,6 @@ #include "tests/AtlasTestEnvironment.h" using namespace atlas::mesh; -using namespace atlas::mesh::temporary; namespace atlas { namespace test { @@ -38,7 +37,7 @@ CASE("hybrid_elements") { idx_t triangle_nodes[] = {1, 5, 3, 1, 5, 2}; - idx_t triags_type_idx = hybrid_elements.add(new Triangle(), 2, triangle_nodes); + idx_t triags_type_idx = hybrid_elements.add(ElementType::create("Triangle"), 2, triangle_nodes); EXPECT(triags_type_idx == 0); @@ -50,7 +49,7 @@ CASE("hybrid_elements") { quad_nodes[2] = 2; quad_nodes[3] = 3; - idx_t quads_type_idx = hybrid_elements.add(new Quadrilateral(), 1, quad_nodes); + idx_t quads_type_idx = hybrid_elements.add(ElementType::create("Quadrilateral"), 1, quad_nodes); EXPECT(quads_type_idx == 1); @@ -151,7 +150,7 @@ CASE("elements") { idx_t triangle_nodes[] = {1, 5, 3, 1, 5, 2}; idx_t triag1[3] = {9, 8, 7}; - Elements elements(new Triangle(), 2, triangle_nodes); + Elements elements(ElementType::create("Triangle"), 2, triangle_nodes); EXPECT(elements.begin() == 0); EXPECT(elements.end() == 2); @@ -247,8 +246,8 @@ CASE("zero_elements") { HybridElements hybrid_elements; idx_t* nodes = nullptr; - hybrid_elements.add(new Triangle(), 0, nodes); - hybrid_elements.add(new Quadrilateral(), 0, nodes); + hybrid_elements.add(ElementType::create("Triangle"), 0, nodes); + hybrid_elements.add(ElementType::create("Quadrilateral"), 0, nodes); EXPECT(hybrid_elements.size() == 0); EXPECT(hybrid_elements.nb_types() == 2); @@ -335,9 +334,9 @@ CASE("cells_insert") { Log::info() << "\n\n\ncells_insert \n============ \n\n" << std::endl; HybridElements cells; idx_t c1[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}; - cells.add(new Quadrilateral(), 3, c1); + cells.add(ElementType::create("Quadrilateral"), 3, c1); idx_t c2[] = {13, 14, 15, 16, 17, 18}; - cells.add(new Triangle(), 2, c2); + cells.add(ElementType::create("Triangle"), 2, c2); EXPECT(cells.elements(0).size() == 3); EXPECT(cells.elements(1).size() == 2); @@ -383,8 +382,8 @@ CASE("cells_insert") { CASE("cells_add_add") { HybridElements cells; - cells.add(new Quadrilateral(), 3); - cells.add(new Triangle(), 2); + cells.add(ElementType::create("Quadrilateral"), 3); + cells.add(ElementType::create("Triangle"), 2); HybridElements::Connectivity& conn = cells.node_connectivity(); diff --git a/src/tests/mesh/test_pentagon_element.cc b/src/tests/mesh/test_pentagon_element.cc index 454aaab0d..2a2ca050d 100644 --- a/src/tests/mesh/test_pentagon_element.cc +++ b/src/tests/mesh/test_pentagon_element.cc @@ -68,7 +68,7 @@ CASE("test_pentagon_like_healpix") { // elements { - mesh.cells().add(new mesh::temporary::Pentagon(), 4); + mesh.cells().add(mesh::ElementType::create("Pentagon"), 4); auto connect = [&](idx_t cell, std::array nodes) { mesh.cells().node_connectivity().block(0).set(cell, nodes.data()); }; diff --git a/src/tests/trans/test_transgeneral.cc b/src/tests/trans/test_transgeneral.cc index c95522519..ea8051e91 100644 --- a/src/tests/trans/test_transgeneral.cc +++ b/src/tests/trans/test_transgeneral.cc @@ -1582,6 +1582,7 @@ CASE("test_trans_levels") { #endif +#if ATLAS_HAVE_TRANS #if ATLAS_HAVE_ECTRANS || defined(TRANS_HAVE_INVTRANS_ADJ) CASE("test_2level_adjoint_test_with_powerspectrum_convolution") { std::string grid_uid("F64"); // Regular Gaussian F ( 8 N^2) @@ -1812,7 +1813,7 @@ CASE("test_2level_adjoint_test_with_vortdiv") { EXPECT(std::abs(adj_value - adj_value2) / adj_value < 1e-12); } #endif - +#endif #if 0 CASE( "test_trans_fourier_truncation" ) { diff --git a/src/tests/util/fctest_kdtree.F90 b/src/tests/util/fctest_kdtree.F90 index acbb8537c..1c7aa142f 100644 --- a/src/tests/util/fctest_kdtree.F90 +++ b/src/tests/util/fctest_kdtree.F90 @@ -229,26 +229,26 @@ end module fcta_KDTree_fixture call kdtree%closestPoints(plonlat(1), plonlat(2), k, indices, distances, lonlats(:,1), lonlats(:,2)) do i = 1, k FCTEST_CHECK_EQUAL( indices(i) , result_indices(i) ) - FCTEST_CHECK_CLOSE( lonlats(i,1) , result_lonlats(i,1) , 1.e-12_c_double ) - FCTEST_CHECK_CLOSE( lonlats(i,2) , result_lonlats(i,2) , 1.e-12_c_double ) - FCTEST_CHECK_CLOSE( distances(i) , result_distances(i) , 1.e-12_c_double ) + FCTEST_CHECK_CLOSE( lonlats(i,1) , result_lonlats(i,1) , 1.e-10_c_double ) + FCTEST_CHECK_CLOSE( lonlats(i,2) , result_lonlats(i,2) , 1.e-10_c_double ) + FCTEST_CHECK_CLOSE( distances(i) , result_distances(i) , 1.e-10_c_double ) end do ! Check closestPoint call kdtree%closestPoint(plonlat(1), plonlat(2), indices(1), distances(1), lonlats(1,1), lonlats(1,2)) FCTEST_CHECK_EQUAL( indices(1) , result_indices(1) ) - FCTEST_CHECK_CLOSE( lonlats(1,1) , result_lonlats(1,1) , 1.e-12_c_double ) - FCTEST_CHECK_CLOSE( lonlats(1,2) , result_lonlats(1,2) , 1.e-12_c_double ) - FCTEST_CHECK_CLOSE( distances(1) , result_distances(1) , 1.e-12_c_double ) + FCTEST_CHECK_CLOSE( lonlats(1,1) , result_lonlats(1,1) , 1.e-10_c_double ) + FCTEST_CHECK_CLOSE( lonlats(1,2) , result_lonlats(1,2) , 1.e-10_c_double ) + FCTEST_CHECK_CLOSE( distances(1) , result_distances(1) , 1.e-10_c_double ) ! Check closestPoints call kdtree%closestPointsWithinRadius(plonlat(1), plonlat(2), 3.5e5_c_double, kk, indices_rad, distances_rad, lons_rad, lats_rad) FCTEST_CHECK_EQUAL( kk , 4 ) do i = 1, kk FCTEST_CHECK_EQUAL( indices_rad(i) , result_indices(i) ) - FCTEST_CHECK_CLOSE( lons_rad(i) , result_lonlats(i,1) , 1.e-12_c_double ) - FCTEST_CHECK_CLOSE( lats_rad(i) , result_lonlats(i,2) , 1.e-12_c_double ) - FCTEST_CHECK_CLOSE( distances_rad(i) , result_distances(i) , 1.e-12_c_double ) + FCTEST_CHECK_CLOSE( lons_rad(i) , result_lonlats(i,1) , 1.e-10_c_double ) + FCTEST_CHECK_CLOSE( lats_rad(i) , result_lonlats(i,2) , 1.e-10_c_double ) + FCTEST_CHECK_CLOSE( distances_rad(i) , result_distances(i) , 1.e-10_c_double ) end do ! Finalization @@ -268,26 +268,26 @@ end module fcta_KDTree_fixture call kdtree%closestPoints(plonlat, k, indices, distances, lonlats) do i = 1, k FCTEST_CHECK_EQUAL( indices(i) , result_indices(i) ) - FCTEST_CHECK_CLOSE( lonlats(i,1) , result_lonlats(i,1) , 1.e-12_c_double ) - FCTEST_CHECK_CLOSE( lonlats(i,2) , result_lonlats(i,2) , 1.e-12_c_double ) - FCTEST_CHECK_CLOSE( distances(i) , result_distances(i) , 1.e-12_c_double ) + FCTEST_CHECK_CLOSE( lonlats(i,1) , result_lonlats(i,1) , 1.e-10_c_double ) + FCTEST_CHECK_CLOSE( lonlats(i,2) , result_lonlats(i,2) , 1.e-10_c_double ) + FCTEST_CHECK_CLOSE( distances(i) , result_distances(i) , 1.e-10_c_double ) end do ! Check closestPoint call kdtree%closestPoint(plonlat, indices(1), distances(1), lonlats(1,:)) FCTEST_CHECK_EQUAL( indices(1) , result_indices(1) ) - FCTEST_CHECK_CLOSE( lonlats(1,1) , result_lonlats(1,1) , 1.e-12_c_double ) - FCTEST_CHECK_CLOSE( lonlats(1,2) , result_lonlats(1,2) , 1.e-12_c_double ) - FCTEST_CHECK_CLOSE( distances(1) , result_distances(1) , 1.e-12_c_double ) + FCTEST_CHECK_CLOSE( lonlats(1,1) , result_lonlats(1,1) , 1.e-10_c_double ) + FCTEST_CHECK_CLOSE( lonlats(1,2) , result_lonlats(1,2) , 1.e-10_c_double ) + FCTEST_CHECK_CLOSE( distances(1) , result_distances(1) , 1.e-10_c_double ) ! Check closestPoints call kdtree%closestPointsWithinRadius(plonlat, 3.5e5_c_double, kk, indices_rad, distances_rad, lonlats_rad) FCTEST_CHECK_EQUAL( kk , 4 ) do i = 1, kk FCTEST_CHECK_EQUAL( indices_rad(i) , result_indices(i) ) - FCTEST_CHECK_CLOSE( lonlats_rad(i,1) , result_lonlats(i,1) , 1.e-12_c_double ) - FCTEST_CHECK_CLOSE( lonlats_rad(i,2) , result_lonlats(i,2) , 1.e-12_c_double ) - FCTEST_CHECK_CLOSE( distances_rad(i) , result_distances(i) , 1.e-12_c_double ) + FCTEST_CHECK_CLOSE( lonlats_rad(i,1) , result_lonlats(i,1) , 1.e-10_c_double ) + FCTEST_CHECK_CLOSE( lonlats_rad(i,2) , result_lonlats(i,2) , 1.e-10_c_double ) + FCTEST_CHECK_CLOSE( distances_rad(i) , result_distances(i) , 1.e-10_c_double ) end do ! Finalization 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}