Add function to build mesh from imported connectivity data (#135)

* Add function to build mesh from imported connectivity data

* Review: set remote_index for nodes, but not cells

* Review: pass connectivities via std library containers

* Update copyright to UCAR

* Move shared test code into helper file

* Clean output filename

* Clean documentation

* Clean inconsistent argument name

* Update src/atlas/mesh/BuildMeshFromConnectivities.cc

Co-authored-by: Willem Deconinck <wdeconinck@me.com>

* Code review: correct use of remote_index

* Code review: change API to class with op()

* Code review: update filenames to better match new class name

* Code review: move new code to mesh namespace

src/atlas/CMakeLists.txt

@@ -323,6 +323,8 @@ mesh/HybridElements.cc
 mesh/HybridElements.h
 mesh/Mesh.cc
 mesh/Mesh.h
+mesh/MeshBuilder.cc
+mesh/MeshBuilder.h
 mesh/Nodes.cc
 mesh/Nodes.h
 mesh/PartitionPolygon.cc

src/atlas/mesh/MeshBuilder.cc

@@ -0,0 +1,144 @@
+/*
+ * (C) Copyright 2023 UCAR.
+ *
+ * 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.
+ */
+
+#include "atlas/mesh/MeshBuilder.h"
+
+#include "atlas/array/MakeView.h"
+#include "atlas/mesh/ElementType.h"
+#include "atlas/mesh/HybridElements.h"
+#include "atlas/mesh/Mesh.h"
+#include "atlas/mesh/Nodes.h"
+#include "atlas/parallel/mpi/mpi.h"
+#include "atlas/util/CoordinateEnums.h"
+
+namespace atlas {
+namespace mesh {
+
+//----------------------------------------------------------------------------------------------------------------------
+
+Mesh MeshBuilder::operator()(const std::vector<double>& lons, const std::vector<double>& lats,
+                             const std::vector<int>& ghosts, const std::vector<gidx_t>& global_indices,
+                             const std::vector<idx_t>& remote_indices, const idx_t remote_index_base,
+                             const std::vector<int>& partitions,
+                             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 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();
+
+    ATLAS_ASSERT(nb_nodes == lons.size());
+    ATLAS_ASSERT(nb_nodes == lats.size());
+    ATLAS_ASSERT(nb_nodes == ghosts.size());
+    ATLAS_ASSERT(nb_nodes == remote_indices.size());
+    ATLAS_ASSERT(nb_nodes == partitions.size());
+    ATLAS_ASSERT(nb_tris == tri_boundary_nodes.size());
+    ATLAS_ASSERT(nb_quads == quad_boundary_nodes.size());
+
+    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());
+}
+
+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 {
+    Mesh mesh{};
+
+    // Populate node data
+
+    mesh.nodes().resize(nb_nodes);
+    auto xy        = array::make_view<double, 2>(mesh.nodes().xy());
+    auto lonlat    = array::make_view<double, 2>(mesh.nodes().lonlat());
+    auto ghost     = array::make_view<int, 1>(mesh.nodes().ghost());
+    auto gidx      = array::make_view<gidx_t, 1>(mesh.nodes().global_index());
+    auto ridx      = array::make_indexview<idx_t, 1>(mesh.nodes().remote_index());
+    auto partition = array::make_view<int, 1>(mesh.nodes().partition());
+    auto halo      = array::make_view<int, 1>(mesh.nodes().halo());
+
+    for (size_t i = 0; i < nb_nodes; ++i) {
+        xy(i, size_t(XX)) = lons[i];
+        xy(i, size_t(YY)) = lats[i];
+        // Identity projection, therefore (lon,lat) = (x,y)
+        lonlat(i, size_t(LON)) = lons[i];
+        lonlat(i, size_t(LAT)) = lats[i];
+        ghost(i)               = ghosts[i];
+        gidx(i)                = global_indices[i];
+        ridx(i)                = remote_indices[i] - remote_index_base;
+        partition(i)           = partitions[i];
+    }
+    halo.assign(0);
+
+    // Populate cell/element data
+
+    // 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());
+    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());
+    const bool add_quads = (sum_nb_quads > 0);
+
+    if (add_tris) {
+        mesh.cells().add(new mesh::temporary::Triangle(), nb_tris);
+    }
+    if (add_quads) {
+        mesh.cells().add(new mesh::temporary::Quadrilateral(), nb_quads);
+    }
+
+    atlas::mesh::HybridElements::Connectivity& node_connectivity = mesh.cells().node_connectivity();
+    auto cells_part = array::make_view<int, 1>(mesh.cells().partition());
+    auto cells_gidx = array::make_view<gidx_t, 1>(mesh.cells().global_index());
+
+    // Find position of idx inside global_indices
+    const auto position_of = [&nb_nodes, &global_indices](const gidx_t idx) {
+        const auto& it = std::find(global_indices, global_indices + nb_nodes, idx);
+        ATLAS_ASSERT(it != global_indices + nb_nodes);
+        return std::distance(global_indices, it);
+    };
+
+    size_t idx = 0;
+    if (add_tris) {
+        idx_t buffer[3];
+        for (size_t tri = 0; tri < nb_tris; ++tri) {
+            for (size_t i = 0; i < 3; ++i) {
+                buffer[i] = position_of(tri_boundary_nodes[3 * tri + i]);
+            }
+            node_connectivity.set(idx, buffer);
+            cells_gidx(idx) = tri_global_indices[tri];
+            idx++;
+        }
+    }
+    if (add_quads) {
+        idx_t buffer[4];
+        for (size_t quad = 0; quad < nb_quads; ++quad) {
+            for (size_t i = 0; i < 4; ++i) {
+                buffer[i] = position_of(quad_boundary_nodes[4 * quad + i]);
+            }
+            node_connectivity.set(idx, buffer);
+            cells_gidx(idx) = quad_global_indices[quad];
+            idx++;
+        }
+    }
+
+    ATLAS_ASSERT(idx == nb_tris + nb_quads);
+
+    cells_part.assign(atlas::mpi::comm().rank());
+
+    return mesh;
+}
+
+//----------------------------------------------------------------------------------------------------------------------
+
+}  // namespace mesh
+}  // namespace atlas

src/atlas/mesh/MeshBuilder.h

@@ -0,0 +1,75 @@
+/*
+ * (C) Copyright 2023 UCAR.
+ *
+ * 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.
+ */
+
+#pragma once
+
+#include "atlas/mesh/Mesh.h"
+#include "atlas/util/Config.h"
+
+#include <array>
+#include <vector>
+
+namespace atlas {
+namespace mesh {
+
+//-----------------------------------------------------------------------------
+
+/**
+ * \brief Construct a Mesh by importing external connectivity data
+ *
+ * Given a list of nodes and corresponding cell-to-node connectivity, sets up a Mesh. Not all Mesh
+ * fields are initialized, but enough are to build halos and construct a NodeColumns FunctionSpace.
+ *
+ * Some limitations of the current implementation (but not inherent):
+ * - can only set up triangle and quad cells.
+ * - 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.
+ */
+class MeshBuilder {
+public:
+    MeshBuilder(const eckit::Configuration& = util::NoConfig()) {}
+
+    /**
+     * \brief C-interface to construct a Mesh from external connectivity data
+     *
+     * The inputs lons, lats, ghosts, global_indices, remote_indices, and partitions are vectors of
+     * size nb_nodes, ranging over the nodes locally owned by (or in the ghost nodes of) the MPI
+     * task. The global index is a uniform labeling of the nodes across all MPI tasks; the remote
+     * index is a remote_index_base-based vector index for the node on its owning task.
+     *
+     * The tri/quad connectivities (boundary_nodes and global_indices) are vectors ranging over the
+     * cells owned by the MPI task. Each cell is defined by a list of nodes defining its boundary;
+     * note that each boundary node must be locally known (whether as an owned of ghost node on the
+     * 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.
+     */
+    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;
+
+    /**
+     * \brief C++-interface to construct a Mesh from external connectivity data
+     *
+     * Provides a wrapper to the C-interface using STL containers.
+     */
+    Mesh operator()(const std::vector<double>& lons, const std::vector<double>& lats, const std::vector<int>& ghosts,
+                    const std::vector<gidx_t>& global_indices, const std::vector<idx_t>& remote_indices,
+                    const idx_t remote_index_base, const std::vector<int>& partitions,
+                    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;
+};
+
+//-----------------------------------------------------------------------------
+
+}  // namespace mesh
+}  // namespace atlas

src/tests/mesh/CMakeLists.txt

@@ -114,6 +114,20 @@ if( TEST atlas_test_cubedsphere_meshgen )
     set_tests_properties( atlas_test_cubedsphere_meshgen PROPERTIES TIMEOUT 30 )
 endif()
 
+ecbuild_add_test( TARGET atlas_test_mesh_builder
+  SOURCES test_mesh_builder.cc
+  LIBS atlas
+  ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT}
+)
+
+ecbuild_add_test( TARGET atlas_test_mesh_builder_parallel
+  MPI        6
+  CONDITION  eckit_HAVE_MPI AND MPI_SLOTS GREATER_EQUAL 6
+  SOURCES test_mesh_builder_parallel.cc
+  LIBS atlas
+  ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT}
+)
+
 ecbuild_add_executable( TARGET atlas_test_mesh_reorder
   SOURCES  test_mesh_reorder.cc
   LIBS     atlas

src/tests/mesh/helper_mesh_builder.h

@@ -0,0 +1,83 @@
+/*
+ * (C) Copyright 2023 UCAR.
+ *
+ * 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.
+ */
+
+#include <algorithm>
+#include <array>
+#include <vector>
+
+#include "atlas/array/MakeView.h"
+#include "atlas/mesh/Elements.h"
+#include "atlas/mesh/HybridElements.h"
+#include "atlas/mesh/Mesh.h"
+#include "atlas/mesh/Nodes.h"
+
+namespace atlas {
+namespace test {
+namespace helper {
+
+void check_mesh_nodes_and_cells(const Mesh& mesh, const std::vector<double>& lons, const std::vector<double>& lats,
+                                const std::vector<int>& ghosts, const std::vector<gidx_t>& global_indices,
+                                const std::vector<idx_t>& remote_indices, const idx_t remote_index_base,
+                                const std::vector<int>& partitions,
+                                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 auto mesh_xy        = array::make_view<double, 2>(mesh.nodes().xy());
+    const auto mesh_lonlat    = array::make_view<double, 2>(mesh.nodes().lonlat());
+    const auto mesh_ghost     = array::make_view<int, 1>(mesh.nodes().ghost());
+    const auto mesh_gidx      = array::make_view<gidx_t, 1>(mesh.nodes().global_index());
+    const auto mesh_ridx      = array::make_indexview<idx_t, 1>(mesh.nodes().remote_index());
+    const auto mesh_partition = array::make_view<int, 1>(mesh.nodes().partition());
+    const auto mesh_halo      = array::make_view<int, 1>(mesh.nodes().halo());
+
+    EXPECT(mesh.nodes().size() == lons.size());
+    for (size_t i = 0; i < mesh.nodes().size(); ++i) {
+        EXPECT(mesh_xy(i, 0) == lons[i]);
+        EXPECT(mesh_xy(i, 1) == lats[i]);
+        EXPECT(mesh_lonlat(i, 0) == lons[i]);
+        EXPECT(mesh_lonlat(i, 1) == lats[i]);
+        EXPECT(mesh_ghost(i) == ghosts[i]);
+        EXPECT(mesh_gidx(i) == global_indices[i]);
+        EXPECT(mesh_ridx(i) == remote_indices[i] - remote_index_base);
+        EXPECT(mesh_partition(i) == partitions[i]);
+        EXPECT(mesh_halo(i) == 0.);
+        // Don't expect (or test) any node-to-cell connectivities
+    }
+
+    EXPECT(mesh.cells().nb_types() == 2);
+    EXPECT(mesh.cells().size() == tri_boundary_nodes.size() + quad_boundary_nodes.size());
+
+    const auto position_of = [&global_indices](const gidx_t idx) {
+        const auto& it = std::find(global_indices.begin(), global_indices.end(), idx);
+        ATLAS_ASSERT(it != global_indices.end());
+        return std::distance(global_indices.begin(), it);
+    };
+
+    // Check triangle cell-to-node connectivities
+    EXPECT(mesh.cells().elements(0).size() == tri_boundary_nodes.size());
+    EXPECT(mesh.cells().elements(0).nb_nodes() == 3);
+    for (size_t tri = 0; tri < mesh.cells().elements(0).size(); ++tri) {
+        for (size_t node = 0; node < mesh.cells().elements(0).nb_nodes(); ++node) {
+            EXPECT(mesh.cells().elements(0).node_connectivity()(tri, node) ==
+                   position_of(tri_boundary_nodes[tri][node]));
+        }
+    }
+    // Check quad cell-to-node connectivities
+    EXPECT(mesh.cells().elements(1).size() == quad_boundary_nodes.size());
+    EXPECT(mesh.cells().elements(1).nb_nodes() == 4);
+    for (size_t quad = 0; quad < mesh.cells().elements(1).size(); ++quad) {
+        for (size_t node = 0; node < mesh.cells().elements(1).nb_nodes(); ++node) {
+            EXPECT(mesh.cells().elements(1).node_connectivity()(quad, node) ==
+                   position_of(quad_boundary_nodes[quad][node]));
+        }
+    }
+}
+
+}  // namespace helper
+}  // namespace test
+}  // namespace atlas