Add atlas_TriangularMeshBuilder with Fortran interface for serial meshes

src/atlas/CMakeLists.txt

@@ -364,6 +364,8 @@ mesh/detail/MeshImpl.cc
 mesh/detail/MeshImpl.h
 mesh/detail/MeshIntf.cc
 mesh/detail/MeshIntf.h
+mesh/detail/MeshBuilderIntf.cc
+mesh/detail/MeshBuilderIntf.h
 mesh/detail/PartitionGraph.cc
 mesh/detail/PartitionGraph.h
 mesh/detail/AccumulateFacets.h

src/atlas/mesh/MeshBuilder.cc

@@ -160,6 +160,33 @@ Mesh MeshBuilder::operator()(size_t nb_nodes, const double lons[], const double
                              const gidx_t tri_global_indices[], size_t nb_quads, const gidx_t quad_boundary_nodes[],
                              const gidx_t quad_global_indices[],
                              const eckit::Configuration& config) const {
+    return this->operator()(
+        nb_nodes, global_indices,
+        lons, lats, 1, 1,
+        lons, lats, 1, 1,
+        ghosts, partitions, remote_indices, remote_index_base,
+        nb_tris, tri_global_indices, tri_boundary_nodes,
+        nb_quads, quad_global_indices, quad_boundary_nodes,
+        config
+    );
+}
+
+Mesh MeshBuilder::operator()(size_t nb_nodes, const gidx_t global_indices[],
+                             const double x[], const double y[], size_t xstride, size_t ystride,
+                             const double lon[], const double lat[], size_t lonstride, size_t latstride,
+                             const int ghosts[], const int partitions[], const idx_t remote_index[], const idx_t remote_index_base,
+                             size_t nb_triags, const gidx_t triag_global_index[], const gidx_t triag_nodes_global[],
+                             size_t nb_quads,  const gidx_t quad_global_index[],  const gidx_t quad_nodes_global[],
+                             const eckit::Configuration& config) const {
+    auto* lons = lon;
+    auto* lats = lat;
+    auto* remote_indices = remote_index;
+    auto nb_tris = nb_triags;
+    auto* tri_boundary_nodes = triag_nodes_global;
+    auto* tri_global_indices = triag_global_index;
+    auto* quad_boundary_nodes = quad_nodes_global;
+    auto* quad_global_indices = quad_global_index;
+
     // 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());
@@ -199,9 +226,8 @@ Mesh MeshBuilder::operator()(size_t nb_nodes, const double lons[], const double
     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)
+        xy(i, size_t(XX))      = x[i];
+        xy(i, size_t(YY))      = y[i];
         lonlat(i, size_t(LON)) = lons[i];
         lonlat(i, size_t(LAT)) = lats[i];
         ghost(i)               = ghosts[i];
@@ -274,5 +300,27 @@ Mesh MeshBuilder::operator()(size_t nb_nodes, const double lons[], const double
 
 //----------------------------------------------------------------------------------------------------------------------
 
+Mesh TriangularMeshBuilder::operator()(size_t nb_nodes,  const gidx_t nodes_global_index[], const double x[], const double y[], const double lon[], const double lat[],
+                                       size_t nb_triags, const gidx_t triangle_global_index[], const gidx_t triangle_nodes_global_index[]) const {
+    std::vector<int> ghost(nb_nodes,0);
+    std::vector<int> partition(nb_nodes,0);
+    std::vector<idx_t> remote_index(nb_nodes);
+    idx_t remote_index_base = 0;
+    std::iota(remote_index.begin(), remote_index.end(), remote_index_base);
+
+    size_t nb_quads = 0;
+    std::vector<gidx_t> quad_nodes_global_index;
+    std::vector<gidx_t> quad_global_index;
+
+    return meshbuilder_(
+        nb_nodes, nodes_global_index,
+        x, y, 1, 1, lon, lat, 1, 1,
+        ghost.data(), partition.data(), remote_index.data(), remote_index_base,
+        nb_triags, triangle_global_index, triangle_nodes_global_index,
+        nb_quads, quad_global_index.data(), quad_nodes_global_index.data());
+}
+
+//----------------------------------------------------------------------------------------------------------------------
+
 }  // namespace mesh
 }  // namespace atlas

src/atlas/mesh/MeshBuilder.h

@@ -35,7 +35,7 @@ namespace mesh {
  *
  * The Mesh's Grid can be initialized via the call operator's optional config argument.
  */
-class MeshBuilder {
+class MeshBuilder : public eckit::Owned {
 public:
     MeshBuilder(const eckit::Configuration& = util::NoConfig()) {}
 
@@ -69,6 +69,15 @@ class MeshBuilder {
                     const gidx_t quad_global_indices[],
                     const eckit::Configuration& config = util::NoConfig()) const;
 
+    Mesh operator()(size_t nb_nodes, const gidx_t global_index[],
+                    const double x[], const double y[], size_t xstride, size_t ystride,
+                    const double lon[], const double lat[], size_t lonstride, size_t latstride,
+                    const int ghost[], const int partition[], const idx_t remote_index[], const idx_t remote_index_base,
+                    size_t nb_triags, const gidx_t triag_global_index[], const gidx_t triag_nodes_global[],
+                    size_t nb_quads,  const gidx_t quad_global_index[],  const gidx_t quad_nodes_global[],
+                    const eckit::Configuration& config = util::NoConfig()) const;
+
+
     /**
      * \brief C++-interface to construct a Mesh from external connectivity data
      *
@@ -86,5 +95,34 @@ class MeshBuilder {
 
 //-----------------------------------------------------------------------------
 
+
+class TriangularMeshBuilder {
+public:
+    TriangularMeshBuilder(const eckit::Configuration& config = util::NoConfig()) :
+        meshbuilder_(config) {}
+
+    /**
+     * \brief C-interface to construct a Triangular Mesh from external connectivity data
+     *
+     * The inputs x, y, lons, lats, ghost, global_index, remote_index, and partition 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 triangle connectivities (boundary_nodes and global_index) 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 gidx_t node_global_index[], const double x[], const double y[], const double lon[], const double lat[],
+                    size_t nb_triags, const gidx_t triangle_global_index[], const gidx_t triangle_nodes_global_index[]) const;
+private:
+    MeshBuilder meshbuilder_;
+};
+
+//-----------------------------------------------------------------------------
+
 }  // namespace mesh
 }  // namespace atlas

src/atlas/mesh/detail/MeshBuilderIntf.cc

@@ -0,0 +1,49 @@
+/*
+ * (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.
+ */
+
+#include "atlas/mesh/detail/MeshBuilderIntf.h"
+#include "atlas/runtime/Exception.h"
+
+namespace atlas {
+namespace mesh {
+
+//----------------------------------------------------------------------------------------------------------------------
+// C wrapper interfaces to C++ routines
+
+TriangularMeshBuilder* atlas__TriangularMeshBuilder__new() {
+    return new TriangularMeshBuilder();
+}
+
+void atlas__TriangularMeshBuilder__delete(TriangularMeshBuilder* This) {
+    ATLAS_ASSERT(This != nullptr, "Cannot access uninitialisd atlas_TriangularMeshBuilder");
+    delete This;
+}
+
+Mesh::Implementation* atlas__TriangularMeshBuilder__operator(TriangularMeshBuilder* This,
+        size_t nb_nodes, const gidx_t node_global_index[], const double x[], const double y[], const double lon[], const double lat[],
+        size_t nb_triags, const gidx_t triangle_global_index[], const gidx_t triangle_nodes_global_index[]) {
+
+    ATLAS_ASSERT(This != nullptr, "Cannot access uninitialisd atlas_TriangularMeshBuilder");
+
+    Mesh::Implementation* m;
+    {
+        Mesh mesh = This->operator()(nb_nodes, node_global_index, x, y, lon, lat,
+                                     nb_triags, triangle_global_index, triangle_nodes_global_index);
+        mesh.get()->attach();
+        m = mesh.get();
+    }
+    m->detach();
+    return m;
+}
+
+//----------------------------------------------------------------------------------------------------------------------
+
+}  // namespace mesh
+}  // namespace atlas

src/atlas/mesh/detail/MeshBuilderIntf.h

@@ -0,0 +1,33 @@
+/*
+ * (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 "atlas/mesh/MeshBuilder.h"
+
+//----------------------------------------------------------------------------------------------------------------------
+
+namespace atlas {
+namespace mesh {
+
+// C wrapper interfaces to C++ routines
+extern "C" {
+TriangularMeshBuilder* atlas__TriangularMeshBuilder__new();
+void atlas__TriangularMeshBuilder__delete(TriangularMeshBuilder* This);
+
+Mesh::Implementation* atlas__TriangularMeshBuilder__operator(TriangularMeshBuilder* This,
+  size_t nb_nodes, const gidx_t node_global_index[], const double x[], const double y[], const double lon[], const double lat[],
+  size_t nb_triags, const gidx_t triangle_global_index[], const gidx_t triangle_nodes_global_index[]);
+}
+
+//----------------------------------------------------------------------------------------------------------------------
+
+}  // namespace mesh
+}  // namespace atlas

src/atlas_f/CMakeLists.txt

@@ -83,6 +83,9 @@ generate_fortran_bindings(FORTRAN_BINDINGS ../atlas/domain/detail/Domain.h
 generate_fortran_bindings(FORTRAN_BINDINGS ../atlas/mesh/detail/MeshIntf.h
     MODULE    atlas_mesh_c_binding
     OUTPUT    atlas_mesh_c_binding.f90 )
+generate_fortran_bindings(FORTRAN_BINDINGS ../atlas/mesh/detail/MeshBuilderIntf.h
+    MODULE    atlas_meshbuilder_c_binding
+    OUTPUT    atlas_meshbuilder_c_binding.f90 )
 generate_fortran_bindings(FORTRAN_BINDINGS ../atlas/mesh/Nodes.h)
 generate_fortran_bindings(FORTRAN_BINDINGS ../atlas/mesh/Connectivity.h)
 generate_fortran_bindings(FORTRAN_BINDINGS ../atlas/mesh/HybridElements.h)
@@ -219,6 +222,7 @@ ecbuild_add_library( TARGET atlas_f
         grid/atlas_GridDistribution_module.F90
         grid/atlas_Vertical_module.F90
         grid/atlas_Partitioner_module.F90
+        mesh/atlas_MeshBuilder_module.F90
         mesh/atlas_MeshGenerator_module.F90
         mesh/atlas_Mesh_module.F90
         mesh/atlas_mesh_Nodes_module.F90

src/atlas_f/atlas_module.F90

@@ -129,6 +129,8 @@ module atlas_module
 use atlas_Nabla_module, only: &
     & atlas_Nabla
 #endif
+use atlas_meshbuilder_module, only : &
+    & atlas_TriangularMeshBuilder
 use atlas_mesh_actions_module, only: &
     & atlas_build_parallel_fields, &
     & atlas_build_nodes_parallel_fields, &

src/atlas_f/mesh/atlas_MeshBuilder_module.F90

@@ -0,0 +1,120 @@
+! (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.
+
+#include "atlas/atlas_f.h"
+
+module atlas_MeshBuilder_module
+
+
+use fckit_owned_object_module, only: fckit_owned_object
+
+implicit none
+
+private :: fckit_owned_object
+
+public :: atlas_TriangularMeshBuilder
+
+private
+
+!-----------------------------!
+! atlas_TriangularMeshBuilder !
+!-----------------------------!
+
+!------------------------------------------------------------------------------
+TYPE, extends(fckit_owned_object) :: atlas_TriangularMeshBuilder
+contains
+  procedure, public :: build => atlas_TriangularMeshBuilder__build
+
+#if FCKIT_FINAL_NOT_INHERITING
+  final :: atlas_TriangularMeshBuilder__final_auto
+#endif
+
+END TYPE atlas_TriangularMeshBuilder
+
+interface atlas_TriangularMeshBuilder
+  module procedure atlas_TriangularMeshBuilder__cptr
+  module procedure atlas_TriangularMeshBuilder__config
+end interface
+
+!------------------------------------------------------------------------------
+
+
+!========================================================
+contains
+!========================================================
+
+
+function atlas_TriangularMeshBuilder__cptr(cptr) result(this)
+  use, intrinsic :: iso_c_binding, only: c_ptr
+  type(atlas_TriangularMeshBuilder) :: this
+  type(c_ptr), intent(in) :: cptr
+  call this%reset_c_ptr( cptr )
+  call this%return()
+end function
+
+function atlas_TriangularMeshBuilder__config(config) result(this)
+  use fckit_c_interop_module, only: c_str
+  use atlas_MeshBuilder_c_binding
+  use atlas_Config_module, only: atlas_Config
+  type(atlas_TriangularMeshBuilder) :: this
+  type(atlas_Config), intent(in), optional :: config
+  call this%reset_c_ptr( atlas__TriangularMeshBuilder__new() )
+  call this%return()
+end function
+
+
+!    Mesh operator()(size_t nb_nodes,  const gidx_t node_global_index[], const double x[], const double y[], const double lon[], const double lat[],
+!                    size_t nb_triags, const gidx_t triangle_global_index[], const gidx_t triangle_nodes_global_index[]) const;
+
+function atlas_TriangularMeshBuilder__build(this, &
+    nb_nodes, node_global_index, x, y, lon, lat, &
+    nb_triags, triag_global_index, triag_nodes) result(mesh)
+  use, intrinsic :: iso_c_binding, only: c_double, c_size_t
+  use atlas_MeshBuilder_c_binding
+  use atlas_Mesh_module, only: atlas_Mesh
+  use atlas_kinds_module, only : ATLAS_KIND_GIDX
+  use fckit_array_module, only : array_strides, array_view1d
+
+  type(atlas_Mesh) :: mesh
+  class(atlas_TriangularMeshBuilder), intent(in) :: this
+  integer, intent(in) :: nb_nodes
+  integer(ATLAS_KIND_GIDX), intent(in) :: node_global_index(:)
+  real(c_double), contiguous, intent(in) :: x(:), y(:), lon(:), lat(:)
+  integer, intent(in) :: nb_triags
+  integer(ATLAS_KIND_GIDX), intent(in) :: triag_global_index(nb_triags)
+  integer(ATLAS_KIND_GIDX), intent(in) :: triag_nodes(3,nb_triags)
+
+  integer(ATLAS_KIND_GIDX), pointer :: triag_nodes_1d(:)
+  triag_nodes_1d => array_view1d( triag_nodes, int(0,ATLAS_KIND_GIDX) )
+
+  call mesh%reset_c_ptr() ! Somehow needed with PGI/16.7 and build-type "bit"
+  mesh = atlas_Mesh( atlas__TriangularMeshBuilder__operator(this%CPTR_PGIBUG_A, &
+       & int(nb_nodes,c_size_t), node_global_index, x, y, lon, lat, &
+       & int(nb_triags,c_size_t), triag_global_index, triag_nodes_1d) )
+  call mesh%return()
+end function
+
+
+!-------------------------------------------------------------------------------
+
+#if FCKIT_FINAL_NOT_INHERITING
+ATLAS_FINAL subroutine atlas_TriangularMeshBuilder__final_auto(this)
+  type(atlas_TriangularMeshBuilder), intent(inout) :: this
+#if FCKIT_FINAL_DEBUGGING
+  write(0,*) "atlas_MeshBuilder__final_auto"
+#endif
+#if FCKIT_FINAL_NOT_PROPAGATING
+  call this%final()
+#endif
+  FCKIT_SUPPRESS_UNUSED( this )
+end subroutine
+#endif
+
+! ----------------------------------------------------------------------------------------
+
+end module atlas_MeshBuilder_module