eckit::geo triangulation

src/eckit/geo/CMakeLists.txt

@@ -50,6 +50,8 @@ list(APPEND eckit_geo_srcs
     Sphere.cc
     Sphere.h
     SphereT.h
+    Triangulation.cc
+    Triangulation.h
     UnitSphere.h
     area/BoundingBox.cc
     area/BoundingBox.h
@@ -140,3 +142,8 @@ if(HAVE_CODEC)
     target_link_libraries(eckit_geo PUBLIC eckit_codec)
 endif()
 
+if(Qhull_FOUND)
+    target_sources(eckit_geo PRIVATE triangulation/Triangulation3.cc triangulation/Triangulation3.h)
+    target_link_libraries(eckit_geo PUBLIC Qhull::Qhull)
+endif()
+

src/eckit/geo/Triangulation.h

@@ -0,0 +1,70 @@
+/*
+ * (C) Copyright 1996- 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.
+ */
+
+
+#include <array>
+#include <vector>
+
+#include "eckit/memory/Builder.h"
+#include "eckit/memory/Factory.h"
+
+
+namespace eckit {
+class Configuration;
+}  // namespace eckit
+
+
+namespace eckit::geo {
+
+
+using Triangle = std::array<size_t, 3>;
+
+
+class Triangulation : protected std::vector<Triangle> {
+public:
+    // -- Types
+
+    using builder_t = BuilderT1<Triangulation>;
+    using ARG1      = const Configuration&;
+
+    // -- Destructor
+
+    virtual ~Triangulation() = default;
+
+    // -- Operators
+
+    using vector::operator[];
+
+    // -- Methods
+
+    static std::string className() { return "Triangulation"; }
+
+    using vector::at;
+    using vector::size;
+
+    using vector::begin;
+    using vector::cbegin;
+    using vector::cend;
+    using vector::end;
+
+protected:
+    // -- Constructors
+
+    Triangulation() /*noexcept */ = default;
+};
+
+
+using TriangulationFactory = Factory<Triangulation>;
+
+template <typename T>
+using TriangulationBuilder = ConcreteBuilderT1<Triangulation, T>;
+
+
+}  // namespace eckit::geo

src/eckit/geo/triangulation/Triangulation3.cc

@@ -0,0 +1,57 @@
+/*
+ * (C) Copyright 1996- 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.
+ */
+
+
+#include "eckit/geo/triangulation/Triangulation3.h"
+
+#include <memory>
+
+#include "eckit/exception/Exceptions.h"
+#include "eckit/geo/Grid.h"
+#include "eckit/geo/projection/LonLatToXYZ.h"
+
+#include "libqhullcpp/Qhull.h"
+#include "libqhullcpp/QhullFacetList.h"
+#include "libqhullcpp/QhullVertexSet.h"
+
+
+namespace eckit::geo::triangulation {
+
+
+Triangulation3::Triangulation3(const Grid& grid) {
+    auto [lat, lon] = grid.to_latlon();
+    ASSERT(lat.size() == lon.size());
+
+    auto N = static_cast<int>(lon.size());
+    std::vector<double> coord(N * 3);
+
+    geo::projection::LonLatToXYZ lonlat2xyz(1.);
+    for (size_t i = 0; i < N; ++i) {
+        auto p = lonlat2xyz.fwd({lon[i], lat[i]});
+
+        coord[i * 3 + 0] = p.X;
+        coord[i * 3 + 1] = p.Y;
+        coord[i * 3 + 2] = p.Z;
+    }
+
+    auto qh = std::make_unique<orgQhull::Qhull>("", 3, N, coord.data(), "Qt");
+
+    reserve(qh->facetCount());
+    for (const auto& facet : qh->facetList()) {
+        auto vertices = facet.vertices();
+        ASSERT(vertices.size() == 3);
+        emplace_back(Triangle{static_cast<size_t>(vertices[0].id()),
+                              static_cast<size_t>(vertices[1].id()),
+                              static_cast<size_t>(vertices[2].id())});
+    }
+}
+
+
+}  // namespace eckit::geo::triangulation

src/eckit/geo/triangulation/Triangulation3.h

@@ -0,0 +1,31 @@
+/*
+ * (C) Copyright 1996- 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 "eckit/geo/Triangulation.h"
+
+
+namespace eckit::geo {
+class Grid;
+}  // namespace eckit::geo
+
+
+namespace eckit::geo::triangulation {
+
+
+class Triangulation3 : public Triangulation {
+public:
+    explicit Triangulation3(const Grid&);
+};
+
+
+}  // namespace eckit::geo::triangulation

src/tools/eckit-grid-triangulation.cc

@@ -9,81 +9,23 @@
  */
 
 
-#include <array>
 #include <memory>
 #include <ostream>
-#include <string>
-#include <vector>
 
 #include "eckit/config/MappedConfiguration.h"
 #include "eckit/exception/Exceptions.h"
 #include "eckit/geo/Grid.h"
 #include "eckit/geo/projection/LonLatToXYZ.h"
+#include "eckit/geo/triangulation/Triangulation3.h"
 #include "eckit/log/Log.h"
 #include "eckit/option/CmdArgs.h"
 #include "eckit/option/EckitTool.h"
 #include "eckit/option/SimpleOption.h"
 
-#include "libqhullcpp/Qhull.h"
-#include "libqhullcpp/QhullFacet.h"
-#include "libqhullcpp/QhullFacetList.h"
-#include "libqhullcpp/QhullPoints.h"
-#include "libqhullcpp/QhullVertexSet.h"
-
 
 namespace eckit::tools {
 
 
-struct Triangulation {
-    using triangle_type = std::array<countT, 3>;
-
-    size_t n_triangles() const {
-        ASSERT(qhull_);
-        return static_cast<size_t>(qhull_->facetCount());
-    }
-
-    std::vector<triangle_type> triangles() const {
-        ASSERT(qhull_);
-
-        std::vector<triangle_type> tri;
-        tri.reserve(qhull_->facetCount());
-
-        for (const auto& facet : qhull_->facetList()) {
-            auto vertices = facet.vertices();
-            ASSERT(vertices.size() == 3);
-
-            tri.emplace_back(triangle_type{vertices[0].id(), vertices[1].id(), vertices[2].id()});
-        }
-
-        return tri;
-    }
-
-protected:
-    std::unique_ptr<const orgQhull::Qhull> qhull_;
-};
-
-
-struct Triangulation3D : Triangulation {
-    Triangulation3D(const std::vector<double>& lon, const std::vector<double>& lat) {
-        ASSERT(lon.size() == lat.size());
-
-        auto N = static_cast<int>(lon.size());
-        std::vector<double> coord(N * 3);
-
-        geo::projection::LonLatToXYZ lonlat2xyz(1.);
-        for (size_t i = 0; i < N; ++i) {
-            auto p = lonlat2xyz.fwd({lon[i], lat[i]});
-
-            coord[i * 3 + 0] = p.X;
-            coord[i * 3 + 1] = p.Y;
-            coord[i * 3 + 2] = p.Z;
-        }
-
-        qhull_ = std::make_unique<orgQhull::Qhull>("", 3, N, coord.data(), "Qt");
-    }
-};
-
-
 struct EckitGridTriangulation final : EckitTool {
     EckitGridTriangulation(int argc, char** argv) :
         EckitTool(argc, argv) {
@@ -104,13 +46,34 @@ void EckitGridTriangulation::execute(const option::CmdArgs& args) {
     std::unique_ptr<const geo::Grid> grid(
         geo::GridFactory::build(MappedConfiguration{{{"grid", args.getString("grid")}}}));
 
+    auto& out = std::cout;
+
+    geo::triangulation::Triangulation3 tri(*grid);
+
     auto [lat, lon] = grid->to_latlon();
-    ASSERT(lon.size() == lat.size());
 
-    std::unique_ptr<Triangulation> tri(new Triangulation3D(lon, lat));
-    for (const auto& t : tri->triangles()) {
-        Log::info() << t[0] << ' ' << t[1] << ' ' << t[2] << std::endl;
+    auto N = grid->size();
+    ASSERT(N == lat.size());
+    ASSERT(N == lon.size());
+
+    out << "$MeshFormat\n"
+           "4.1 0 8\n"
+           "$EndMeshFormat\n";
+
+    out << "$Nodes\n" << N << '\n';
+    geo::projection::LonLatToXYZ lonlat2xyz(1.);
+    for (size_t i = 0; i < N; ++i) {
+        auto p = lonlat2xyz.fwd(geo::PointLonLat{lon[i], lat[i]});
+        out << (i + 1) << ' ' << p.X << ' ' << p.Y << ' ' << p.Z << '\n';
+    }
+    out << "$EndNodes\n";
+
+    out << "$Elements\n" << tri.size() << '\n';
+    size_t i = 0;
+    for (const auto& t : tri) {
+        out << (i + 1) << " 2 4 1 1 1 0 " << t[0] << ' ' << t[1] << ' ' << t[2] << '\n';
     }
+    out << "$EndElements\n";
 };