Enable MeshBuilder to set up the Grid (#152)

src/atlas/mesh/Mesh.h

@@ -39,6 +39,7 @@ class PartitionPolygons;
 
 namespace atlas {
 namespace mesh {
+class MeshBuilder;
 class Nodes;
 class HybridElements;
 typedef HybridElements Edges;
@@ -140,6 +141,7 @@ class Mesh : DOXYGEN_HIDE(public util::ObjectHandle<mesh::detail::MeshImpl>) {
         return s;
     }
 
+    friend class mesh::MeshBuilder;
     friend class meshgenerator::MeshGeneratorImpl;
     void setProjection(const Projection& p) { get()->setProjection(p); }
     void setGrid(const Grid& p) { get()->setGrid(p); }

src/atlas/mesh/MeshBuilder.cc

@@ -7,7 +7,13 @@
 
 #include "atlas/mesh/MeshBuilder.h"
 
+#include <algorithm>
+
 #include "atlas/array/MakeView.h"
+#include "atlas/grid/Grid.h"
+#include "atlas/grid/Iterator.h"
+#include "atlas/grid/StructuredGrid.h"
+#include "atlas/grid/UnstructuredGrid.h"
 #include "atlas/mesh/ElementType.h"
 #include "atlas/mesh/HybridElements.h"
 #include "atlas/mesh/Mesh.h"
@@ -16,6 +22,106 @@
 #include "atlas/util/CoordinateEnums.h"
 
 namespace atlas {
+
+namespace detail {
+atlas::UnstructuredGrid assemble_unstructured_grid(size_t nb_nodes, const double lons[], const double lats[],
+                                                   const int ghosts[], const eckit::mpi::Comm& comm) {
+
+    // First serialize owned lons and lats into single vector
+    const size_t nb_owned_nodes = std::count(ghosts, ghosts + nb_nodes, 0);
+    std::vector<double> owned_lonlats(2 * nb_owned_nodes);
+    int counter = 0;
+    for (size_t i = 0; i < nb_nodes; ++i) {
+        if (ghosts[i] == 0) {
+            owned_lonlats[counter] = lons[i];
+            counter++;
+            owned_lonlats[counter] = lats[i];
+            counter++;
+        }
+    }
+    ATLAS_ASSERT(counter == 2 * nb_owned_nodes);
+
+    // Gather points across MPI ranks
+    size_t nb_nodes_global = 0;
+    comm.allReduce(nb_owned_nodes, nb_nodes_global, eckit::mpi::sum());
+
+    std::vector<double> global_lonlats(2 * nb_nodes_global);
+    eckit::mpi::Buffer<double> buffer(comm.size());
+    comm.allGatherv(owned_lonlats.begin(), owned_lonlats.end(), buffer);
+    global_lonlats = std::move(buffer.buffer);
+
+    std::vector<atlas::PointXY> points(nb_nodes_global);
+    for (size_t i = 0; i < nb_nodes_global; ++i) {
+        points[i] = atlas::PointXY({global_lonlats[2*i], global_lonlats[2*i + 1]});
+    }
+
+    return atlas::UnstructuredGrid(new std::vector<atlas::PointXY>(points.begin(), points.end()));
+}
+
+// Check if grid points match mesh points -- not obviously true as the MeshBuilder sets the Grid
+// independently from the Mesh. This check can give confidence that the grid and mesh are specified
+// consistently with each other.
+//
+// This check makes a few assumptions
+// - the global_indices passed to the MeshBuilder are a 1-based contiguous index
+// - the Grid is one of StructuredGrid or UnstructedGrid. This restriction is because a
+//   CubedSphereGrid doesn't present a simple interface for the coordinates at the N'th point.
+// - the Mesh grid points are the grid nodes, and not the grid cell-centers (e.g., HEALPix or CubedSphere meshes)
+void validate_grid_vs_mesh(const atlas::Grid& grid, size_t nb_nodes, const double lons[], const double lats[],
+                           const int ghosts[], const gidx_t global_indices[], const eckit::mpi::Comm& comm) {
+    // Check assumption that global_indices look like a 1-based contiguous index
+    const size_t nb_owned_nodes = std::count(ghosts, ghosts + nb_nodes, 0);
+    size_t nb_nodes_global = 0;
+    comm.allReduce(nb_owned_nodes, nb_nodes_global, eckit::mpi::sum());
+    for (size_t i = 0; i < nb_nodes; ++i) {
+        if (ghosts[i] == 0) {
+            // Check global_indices is consistent with a 1-based contiguous index over nodes
+            ATLAS_ASSERT(global_indices[i] >= 1);
+            ATLAS_ASSERT(global_indices[i] <= nb_nodes_global);
+        }
+    }
+
+    double lonlat[2];
+    const auto equal_within_roundoff = [](const double a, const double b) -> bool {
+        return fabs(a - b) <= 360.0 * 1.0e-16;
+    };
+
+    // Check lonlats for each supported grid type
+    if (grid.type() == "unstructured") {
+        const atlas::UnstructuredGrid ugrid(grid);
+        for (size_t i = 0; i < nb_nodes; ++i) {
+            if (ghosts[i] == 0) {
+                ugrid.lonlat(global_indices[i] - 1, lonlat);
+                if (!equal_within_roundoff(lonlat[0], lons[i]) || !equal_within_roundoff(lonlat[1], lats[i])) {
+                    throw_Exception("In MeshBuilder: UnstructuredGrid from config does not match mesh coordinates", Here());
+                }
+            }
+        }
+    } else if (grid.type() == "structured") {
+        const atlas::StructuredGrid sgrid(grid);
+        for (size_t i = 0; i < nb_nodes; ++i) {
+            if (ghosts[i] == 0) {
+                idx_t i, j;
+                sgrid.index2ij(global_indices[i] - 1, i, j);
+                sgrid.lonlat(i, j, lonlat);
+                if (!equal_within_roundoff(lonlat[0], lons[i]) || !equal_within_roundoff(lonlat[1], lats[i])) {
+                    throw_Exception("In MeshBuilder: StructuredGrid from config does not match mesh coordinates", Here());
+                }
+            }
+        }
+    } else {
+        for (size_t i = 0; i < nb_nodes; ++i) {
+            if (ghosts[i] == 0) {
+                auto point = *(grid.lonlat().begin() + static_cast<size_t>(global_indices[i] - 1));
+                if (!equal_within_roundoff(point.lon(), lons[i]) || !equal_within_roundoff(point.lat(), lats[i])) {
+                    throw_Exception("In MeshBuilder: Grid from config does not match mesh coordinates", Here());
+                }
+            }
+        }
+    }
+}
+}  // namespace detail
+
 namespace mesh {
 
 //----------------------------------------------------------------------------------------------------------------------
@@ -27,7 +133,8 @@ Mesh MeshBuilder::operator()(const std::vector<double>& lons, const std::vector<
                              const std::vector<std::array<gidx_t, 3>>& tri_boundary_nodes,
                              const std::vector<gidx_t>& tri_global_indices,
                              const std::vector<std::array<gidx_t, 4>>& quad_boundary_nodes,
-                             const std::vector<gidx_t>& quad_global_indices) const {
+                             const std::vector<gidx_t>& quad_global_indices,
+                             const eckit::Configuration& config) const {
     const size_t nb_nodes = global_indices.size();
     const size_t nb_tris  = tri_global_indices.size();
     const size_t nb_quads = quad_global_indices.size();
@@ -43,16 +150,42 @@ Mesh MeshBuilder::operator()(const std::vector<double>& lons, const std::vector<
     return operator()(nb_nodes, lons.data(), lats.data(), ghosts.data(), global_indices.data(), remote_indices.data(),
                       remote_index_base, partitions.data(), nb_tris,
                       reinterpret_cast<const gidx_t*>(tri_boundary_nodes.data()), tri_global_indices.data(), nb_quads,
-                      reinterpret_cast<const gidx_t*>(quad_boundary_nodes.data()), quad_global_indices.data());
+                      reinterpret_cast<const gidx_t*>(quad_boundary_nodes.data()), quad_global_indices.data(),
+                      config);
 }
 
 Mesh MeshBuilder::operator()(size_t nb_nodes, const double lons[], const double lats[], const int ghosts[],
                              const gidx_t global_indices[], const idx_t remote_indices[], const idx_t remote_index_base,
                              const int partitions[], size_t nb_tris, const gidx_t tri_boundary_nodes[],
                              const gidx_t tri_global_indices[], size_t nb_quads, const gidx_t quad_boundary_nodes[],
-                             const gidx_t quad_global_indices[]) const {
+                             const gidx_t quad_global_indices[],
+                             const eckit::Configuration& config) const {
+    // Get MPI comm from config name or fall back to atlas default comm
+    auto mpi_comm_name = [](const auto& config) {
+        return config.getString("mpi_comm", atlas::mpi::comm().name());
+    };
+    const eckit::mpi::Comm& comm = eckit::mpi::comm(mpi_comm_name(config).c_str());
+
     Mesh mesh{};
 
+    // Setup a grid, if requested via config argument
+    if (config.has("grid")) {
+        atlas::Grid grid;
+        if (config.has("grid.type") && (config.getString("grid.type") == "unstructured")
+            && !config.has("grid.xy")) {
+            // Assemble the unstructured grid by gathering input lons,lats across ranks
+            grid = ::atlas::detail::assemble_unstructured_grid(nb_nodes, lons, lats, ghosts, comm);
+        } else {
+            // Build grid directly from config
+            grid = atlas::Grid(config.getSubConfiguration("grid"));
+            const bool validate = config.getBool("validate", false);
+            if (validate) {
+                ::atlas::detail::validate_grid_vs_mesh(grid, nb_nodes, lons, lats, ghosts, global_indices, comm);
+            }
+        }
+        mesh.setGrid(grid);
+    }
+
     // Populate node data
 
     mesh.nodes().resize(nb_nodes);
@@ -82,11 +215,11 @@ Mesh MeshBuilder::operator()(size_t nb_nodes, const double lons[], const double
     // First, count how many cells of each type are on this processor
     // Then optimize away the element type if globally nb_tris or nb_quads is zero
     size_t sum_nb_tris = 0;
-    atlas::mpi::comm().allReduce(nb_tris, sum_nb_tris, eckit::mpi::sum());
+    comm.allReduce(nb_tris, sum_nb_tris, eckit::mpi::sum());
     const bool add_tris = (sum_nb_tris > 0);
 
     size_t sum_nb_quads = 0;
-    atlas::mpi::comm().allReduce(nb_quads, sum_nb_quads, eckit::mpi::sum());
+    comm.allReduce(nb_quads, sum_nb_quads, eckit::mpi::sum());
     const bool add_quads = (sum_nb_quads > 0);
 
     if (add_tris) {
@@ -133,7 +266,7 @@ Mesh MeshBuilder::operator()(size_t nb_nodes, const double lons[], const double
 
     ATLAS_ASSERT(idx == nb_tris + nb_quads);
 
-    cells_part.assign(atlas::mpi::comm().rank());
+    cells_part.assign(comm.rank());
 
     return mesh;
 }

src/atlas/mesh/MeshBuilder.h

@@ -8,8 +8,11 @@
 #pragma once
 
 #include "atlas/mesh/Mesh.h"
+#include "atlas/parallel/mpi/mpi.h"
 #include "atlas/util/Config.h"
 
+#include "eckit/config/Configuration.h"
+
 #include <array>
 #include <vector>
 
@@ -29,6 +32,8 @@ namespace mesh {
  * - cannot import halos, i.e., cells owned by other MPI tasks; halos can still be subsequently
  *   computed by calling the BuildMesh action.
  * - cannot import node-to-cell connectivity information.
+ *
+ * The Mesh's Grid can be initialized via the call operator's optional config argument.
  */
 class MeshBuilder {
 public:
@@ -48,12 +53,21 @@ class MeshBuilder {
      * MPI task), in other words, must be an element of the node global_indices. The boundary nodes
      * are ordered node-varies-fastest, element-varies-slowest order. The cell global index is,
      * here also, a uniform labeling over the of the cells across all MPI tasks.
+     *
+     * The config argument can be used to
+     * - Request the Mesh's Grid to be constructed, usually from the config. If the Grid is either
+     *   an UnstructuredGrid or a Structured grid, the `validate` bool option can be used to trigger
+     *   a simple check that the grid is consistent with the lons/lats passed in to the MeshBuilder.
+     *   In the special case where `grid.type` is unstructured and the `grid.xy` coordinates are
+     *   _not_ given, then the grid is constructed from the lons/lats passed to the MeshBuilder.
+     * - Select which MPI communicator to use.
      */
     Mesh operator()(size_t nb_nodes, const double lons[], const double lats[], const int ghosts[],
                     const gidx_t global_indices[], const idx_t remote_indices[], const idx_t remote_index_base,
                     const int partitions[], size_t nb_tris, const gidx_t tri_boundary_nodes[],
                     const gidx_t tri_global_indices[], size_t nb_quads, const gidx_t quad_boundary_nodes[],
-                    const gidx_t quad_global_indices[]) const;
+                    const gidx_t quad_global_indices[],
+                    const eckit::Configuration& config = util::NoConfig()) const;
 
     /**
      * \brief C++-interface to construct a Mesh from external connectivity data
@@ -66,7 +80,8 @@ class MeshBuilder {
                     const std::vector<std::array<gidx_t, 3>>& tri_boundary_nodes,
                     const std::vector<gidx_t>& tri_global_indices,
                     const std::vector<std::array<gidx_t, 4>>& quad_boundary_nodes,
-                    const std::vector<gidx_t>& quad_global_indices) const;
+                    const std::vector<gidx_t>& quad_global_indices,
+                    const eckit::Configuration& config = util::NoConfig()) const;
 };
 
 //-----------------------------------------------------------------------------

src/tests/mesh/test_mesh_builder.cc

@@ -10,11 +10,13 @@
 #include <vector>
 
 #include "atlas/array/MakeView.h"
+#include "atlas/grid.h"
 #include "atlas/mesh/Elements.h"
 #include "atlas/mesh/HybridElements.h"
 #include "atlas/mesh/Mesh.h"
 #include "atlas/mesh/MeshBuilder.h"
 #include "atlas/mesh/Nodes.h"
+#include "atlas/util/Config.h"
 
 #include "tests/AtlasTestEnvironment.h"
 #include "tests/mesh/helper_mesh_builder.h"
@@ -129,15 +131,78 @@ CASE("test_cs_c2_mesh_serial") {
     std::iota(quad_global_indices.begin(), quad_global_indices.end(), 9);  // nb_tris + 1
 
     const MeshBuilder mesh_builder{};
-    const Mesh mesh = mesh_builder(lons, lats, ghosts, global_indices, remote_indices, remote_index_base, partitions,
-                                   tri_boundary_nodes, tri_global_indices, quad_boundary_nodes, quad_global_indices);
-
-    helper::check_mesh_nodes_and_cells(mesh, lons, lats, ghosts, global_indices, remote_indices, remote_index_base,
-                                       partitions, tri_boundary_nodes, tri_global_indices, quad_boundary_nodes,
-                                       quad_global_indices);
 
-    //Gmsh gmsh("out.msh", util::Config("coordinates", "xyz"));
-    //gmsh.write(mesh);
+    SECTION("Build Mesh without a Grid") {
+        const Mesh mesh = mesh_builder(lons, lats, ghosts, global_indices, remote_indices, remote_index_base, partitions,
+                                       tri_boundary_nodes, tri_global_indices, quad_boundary_nodes, quad_global_indices);
+
+        helper::check_mesh_nodes_and_cells(mesh, lons, lats, ghosts, global_indices, remote_indices, remote_index_base,
+                                           partitions, tri_boundary_nodes, tri_global_indices, quad_boundary_nodes,
+                                           quad_global_indices);
+        EXPECT(!mesh.grid());
+
+        //Gmsh gmsh("out.msh", util::Config("coordinates", "xyz"));
+        //gmsh.write(mesh);
+    }
+
+    SECTION("Build Mesh with an UnstructuredGrid from config") {
+        std::vector<double> lonlats(2 * lons.size());
+        int counter = 0;
+        for (size_t i = 0; i < lons.size(); ++i) {
+            lonlats[counter] = lons[i];
+            counter++;
+            lonlats[counter] = lats[i];
+            counter++;
+        }
+
+        util::Config config{};
+        config.set("grid.type", "unstructured");
+        config.set("grid.xy", lonlats);
+        config.set("validate", true);
+
+        const Mesh mesh = mesh_builder(lons, lats, ghosts, global_indices, remote_indices, remote_index_base, partitions,
+                                       tri_boundary_nodes, tri_global_indices, quad_boundary_nodes, quad_global_indices,
+                                       config);
+
+        helper::check_mesh_nodes_and_cells(mesh, lons, lats, ghosts, global_indices, remote_indices, remote_index_base,
+                                           partitions, tri_boundary_nodes, tri_global_indices, quad_boundary_nodes,
+                                           quad_global_indices);
+
+        EXPECT(mesh.grid());
+        EXPECT(mesh.grid().type() == "unstructured");
+    }
+
+    SECTION("Build Mesh with an UnstructuredGrid, with Grid assembled from MeshBuilder arguments") {
+        util::Config config{};
+        config.set("grid.type", "unstructured");
+
+        const Mesh mesh = mesh_builder(lons, lats, ghosts, global_indices, remote_indices, remote_index_base, partitions,
+                                       tri_boundary_nodes, tri_global_indices, quad_boundary_nodes, quad_global_indices,
+                                       config);
+
+        helper::check_mesh_nodes_and_cells(mesh, lons, lats, ghosts, global_indices, remote_indices, remote_index_base,
+                                           partitions, tri_boundary_nodes, tri_global_indices, quad_boundary_nodes,
+                                           quad_global_indices);
+
+        EXPECT(mesh.grid());
+        EXPECT(mesh.grid().type() == "unstructured");
+    }
+
+    SECTION("Build Mesh with a CubedSphereGrid from config") {  // but can't validate this
+        util::Config config{};
+        config.set("grid.name", "CS-LFR-2");
+
+        const Mesh mesh = mesh_builder(lons, lats, ghosts, global_indices, remote_indices, remote_index_base, partitions,
+                                       tri_boundary_nodes, tri_global_indices, quad_boundary_nodes, quad_global_indices,
+                                       config);
+
+        helper::check_mesh_nodes_and_cells(mesh, lons, lats, ghosts, global_indices, remote_indices, remote_index_base,
+                                           partitions, tri_boundary_nodes, tri_global_indices, quad_boundary_nodes,
+                                           quad_global_indices);
+
+        EXPECT(mesh.grid());
+        EXPECT(mesh.grid().type() == "cubedsphere");
+    }
 }
 
 //-----------------------------------------------------------------------------