From 803fe2a8ab58b1ee294e22635bfb5664ef7ade67 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Wed, 8 Nov 2023 17:08:32 +0100 Subject: [PATCH] Add atlas_TriangularMeshBuilder with Fortran interface for serial meshes --- src/atlas/CMakeLists.txt | 2 + src/atlas/mesh/MeshBuilder.cc | 54 +- src/atlas/mesh/MeshBuilder.h | 40 +- src/atlas/mesh/actions/BuildHaloDistance.cc | 1497 +++++++++++++++++ src/atlas/mesh/actions/BuildHaloDistance.h | 69 + src/atlas/mesh/detail/MeshBuilderIntf.cc | 49 + src/atlas/mesh/detail/MeshBuilderIntf.h | 33 + src/atlas_f/CMakeLists.txt | 4 + src/atlas_f/atlas_module.F90 | 2 + src/atlas_f/mesh/atlas_MeshBuilder_module.F90 | 120 ++ src/tests/mesh/CMakeLists.txt | 14 + .../fctest_mesh_triangular_mesh_builder.F90 | 132 ++ .../mesh/test_mesh_triangular_mesh_builder.cc | 73 + 13 files changed, 2085 insertions(+), 4 deletions(-) create mode 100644 src/atlas/mesh/actions/BuildHaloDistance.cc create mode 100644 src/atlas/mesh/actions/BuildHaloDistance.h create mode 100644 src/atlas/mesh/detail/MeshBuilderIntf.cc create mode 100644 src/atlas/mesh/detail/MeshBuilderIntf.h create mode 100644 src/atlas_f/mesh/atlas_MeshBuilder_module.F90 create mode 100644 src/tests/mesh/fctest_mesh_triangular_mesh_builder.F90 create mode 100644 src/tests/mesh/test_mesh_triangular_mesh_builder.cc diff --git a/src/atlas/CMakeLists.txt b/src/atlas/CMakeLists.txt index 2b92a933b..3655519be 100644 --- a/src/atlas/CMakeLists.txt +++ b/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 diff --git a/src/atlas/mesh/MeshBuilder.cc b/src/atlas/mesh/MeshBuilder.cc index 32020d997..40155399f 100644 --- a/src/atlas/mesh/MeshBuilder.cc +++ b/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(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 ghost(nb_nodes,0); + std::vector partition(nb_nodes,0); + std::vector 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 quad_nodes_global_index; + std::vector 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 diff --git a/src/atlas/mesh/MeshBuilder.h b/src/atlas/mesh/MeshBuilder.h index b7ab42547..d5c1b2304 100644 --- a/src/atlas/mesh/MeshBuilder.h +++ b/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 diff --git a/src/atlas/mesh/actions/BuildHaloDistance.cc b/src/atlas/mesh/actions/BuildHaloDistance.cc new file mode 100644 index 000000000..6193002d7 --- /dev/null +++ b/src/atlas/mesh/actions/BuildHaloDistance.cc @@ -0,0 +1,1497 @@ +/* + * (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. + */ + + +#include "BuildHaloDistance.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "atlas/array.h" +#include "atlas/array/IndexView.h" +#include "atlas/field/Field.h" +#include "atlas/library/config.h" +#include "atlas/mesh/Elements.h" +#include "atlas/mesh/HybridElements.h" +#include "atlas/mesh/Mesh.h" +#include "atlas/mesh/Nodes.h" +#include "atlas/mesh/actions/BuildParallelFields.h" +#include "atlas/mesh/detail/AccumulateFacets.h" +#include "atlas/parallel/mpi/Buffer.h" +#include "atlas/parallel/mpi/mpi.h" +#include "atlas/parallel/omp/omp.h" +#include "atlas/parallel/omp/sort.h" +#include "atlas/runtime/Exception.h" +#include "atlas/runtime/Log.h" +#include "atlas/runtime/Trace.h" +#include "atlas/util/CoordinateEnums.h" +#include "atlas/util/LonLatMicroDeg.h" +#include "atlas/util/MicroDeg.h" +#include "atlas/util/PeriodicTransform.h" +#include "atlas/util/Unique.h" + +//#define DEBUG_OUTPUT +#ifdef DEBUG_OUTPUT +#include "atlas/mesh/actions/BuildParallelFields.h" +#include "atlas/mesh/actions/BuildXYZField.h" +#include "atlas/output/Gmsh.h" +#endif + +// #define ATLAS_103 +// #define ATLAS_103_SORT + +using atlas::mesh::detail::accumulate_facets; +using atlas::util::LonLatMicroDeg; +using atlas::util::microdeg; +using atlas::util::PeriodicTransform; +using atlas::util::UniqueLonLat; +using Topology = atlas::mesh::Nodes::Topology; + +namespace atlas { +namespace mesh { +namespace actions { + +namespace { +struct Entity { + Entity(gidx_t gid, idx_t idx) { + g = gid; + i = idx; + } + gidx_t g; + idx_t i; + bool operator<(const Entity& other) const { return (g < other.g); } +}; + +void make_nodes_global_index_human_readable(const mesh::actions::BuildHaloDistance& build_halo, mesh::Nodes& nodes, + bool do_all) { + ATLAS_TRACE(); + // TODO: ATLAS-14: fix renumbering of EAST periodic boundary points + // --> Those specific periodic points at the EAST boundary are not checked for + // uid, + // and could receive different gidx for different tasks + + // unused // int mypart = mpi::rank(); + int nparts = mpi::size(); + size_t root = 0; + + array::ArrayView nodes_glb_idx = array::make_view(nodes.global_index()); + // nodes_glb_idx.dump( Log::info() ); + // ATLAS_DEBUG( "min = " << nodes.global_index().metadata().getLong("min") ); + // ATLAS_DEBUG( "max = " << nodes.global_index().metadata().getLong("max") ); + // ATLAS_DEBUG( "human_readable = " << + // nodes.global_index().metadata().getBool("human_readable") ); + gidx_t glb_idx_max = 0; + + std::vector points_to_edit; + + if (do_all) { + points_to_edit.resize(nodes_glb_idx.size()); + for (idx_t i = 0; i < nodes_glb_idx.size(); ++i) { + points_to_edit[i] = i; + } + } + else { + glb_idx_max = nodes.global_index().metadata().getLong("max", 0); + points_to_edit.resize(build_halo.periodic_points_local_index_.size()); + for (size_t i = 0; i < points_to_edit.size(); ++i) { + points_to_edit[i] = build_halo.periodic_points_local_index_[i]; + } + } + + std::vector glb_idx(points_to_edit.size()); + int nb_nodes = static_cast(glb_idx.size()); + for (idx_t i = 0; i < nb_nodes; ++i) { + glb_idx[i] = nodes_glb_idx(points_to_edit[i]); + } + + // ATLAS_DEBUG_VAR( points_to_edit ); + // ATLAS_DEBUG_VAR( points_to_edit.size() ); + + /* + * Sorting following gidx will define global order of + * gathered fields. Special care needs to be taken for + * pole edges, as their centroid might coincide with + * other edges + */ + // Try to recover + // { + // UniqueLonLat compute_uid(nodes); + // for( int jnode=0; jnode recvcounts(mpi::size()); + std::vector recvdispls(mpi::size()); + + ATLAS_TRACE_MPI(GATHER) { mpi::comm().gather(nb_nodes, recvcounts, root); } + idx_t glb_nb_nodes = std::accumulate(recvcounts.begin(), recvcounts.end(), 0); + + recvdispls[0] = 0; + for (int jpart = 1; jpart < nparts; ++jpart) { // start at 1 + recvdispls[jpart] = recvcounts[jpart - 1] + recvdispls[jpart - 1]; + } + + std::vector glb_idx_gathered(glb_nb_nodes); + ATLAS_TRACE_MPI(GATHER) { + mpi::comm().gatherv(glb_idx.data(), glb_idx.size(), glb_idx_gathered.data(), recvcounts.data(), + recvdispls.data(), root); + } + + // 2) Sort all global indices, and renumber from 1 to glb_nb_edges + std::vector node_sort; + node_sort.reserve(glb_nb_nodes); + const idx_t nb_glb_idx_gathered = static_cast(glb_idx_gathered.size()); + for (idx_t jnode = 0; jnode < nb_glb_idx_gathered; ++jnode) { + node_sort.emplace_back(glb_idx_gathered[jnode], jnode); + } + + ATLAS_TRACE_SCOPE("sort on rank 0") { std::sort(node_sort.begin(), node_sort.end()); } + + gidx_t gid = glb_idx_max + 1; + const idx_t nb_node_sort = static_cast(node_sort.size()); + for (idx_t jnode = 0; jnode < nb_node_sort; ++jnode) { + if (jnode > 0 && node_sort[jnode].g != node_sort[jnode - 1].g) { + ++gid; + } + int inode = node_sort[jnode].i; + glb_idx_gathered[inode] = gid; + } + + // 3) Scatter renumbered back + ATLAS_TRACE_MPI(SCATTER) { + mpi::comm().scatterv(glb_idx_gathered.data(), recvcounts.data(), recvdispls.data(), glb_idx.data(), + static_cast(glb_idx.size()), root); + } + + for (int jnode = 0; jnode < nb_nodes; ++jnode) { + nodes_glb_idx(points_to_edit[jnode]) = glb_idx[jnode]; + } + + // nodes_glb_idx.dump( Log::info() ); + // Log::info() << std::endl; + nodes.global_index().metadata().set("human_readable", true); +} + +// ------------------------------------------------------------------------------------- + +void make_cells_global_index_human_readable(const mesh::actions::BuildHaloDistance& build_halo, mesh::Cells& cells, + bool do_all) { + ATLAS_TRACE(); + + int nparts = static_cast(mpi::size()); + idx_t root = 0; + + array::ArrayView cells_glb_idx = array::make_view(cells.global_index()); + // ATLAS_DEBUG( "min = " << cells.global_index().metadata().getLong("min") ); + // ATLAS_DEBUG( "max = " << cells.global_index().metadata().getLong("max") ); + // ATLAS_DEBUG( "human_readable = " << + // cells.global_index().metadata().getBool("human_readable") ); + gidx_t glb_idx_max = 0; + + std::vector cells_to_edit; + + if (do_all) { + cells_to_edit.resize(cells_glb_idx.size()); + for (idx_t i = 0; i < cells_glb_idx.size(); ++i) { + cells_to_edit[i] = i; + } + } + else { + idx_t nb_cells_to_edit(0); + for (const auto& new_cells : build_halo.periodic_cells_local_index_) { + nb_cells_to_edit += new_cells.size(); + } + cells_to_edit.resize(nb_cells_to_edit); + int c{0}; + int i{0}; + for (int t = 0; t < cells.nb_types(); ++t) { + for (idx_t p : build_halo.periodic_cells_local_index_[t]) { + cells_to_edit[i++] = c + p; + } + c += cells.elements(t).size(); + } + glb_idx_max = cells.global_index().metadata().getLong("max", 0); + } + + std::vector glb_idx(cells_to_edit.size()); + const int nb_cells = static_cast(glb_idx.size()); + for (int i = 0; i < nb_cells; ++i) { + glb_idx[i] = cells_glb_idx(cells_to_edit[i]); + } + + // 1) Gather all global indices, together with location + + std::vector recvcounts(mpi::size()); + std::vector recvdispls(mpi::size()); + + ATLAS_TRACE_MPI(GATHER) { mpi::comm().gather(nb_cells, recvcounts, root); } + idx_t glb_nb_cells = std::accumulate(recvcounts.begin(), recvcounts.end(), 0); + + recvdispls[0] = 0; + for (int jpart = 1; jpart < nparts; ++jpart) { // start at 1 + recvdispls[jpart] = recvcounts[jpart - 1] + recvdispls[jpart - 1]; + } + + std::vector glb_idx_gathered(glb_nb_cells); + ATLAS_TRACE_MPI(GATHER) { + mpi::comm().gatherv(glb_idx.data(), glb_idx.size(), glb_idx_gathered.data(), recvcounts.data(), + recvdispls.data(), root); + } + + // 2) Sort all global indices, and renumber from 1 to glb_nb_edges + std::vector cell_sort; + cell_sort.reserve(glb_nb_cells); + for (idx_t jcell = 0; jcell < glb_nb_cells; ++jcell) { + cell_sort.emplace_back(glb_idx_gathered[jcell], jcell); + } + + ATLAS_TRACE_SCOPE("sort on rank 0") { std::sort(cell_sort.begin(), cell_sort.end()); } + + gidx_t gid = glb_idx_max + 1; + for (idx_t jcell = 0; jcell < glb_nb_cells; ++jcell) { + if (jcell > 0 && cell_sort[jcell].g != cell_sort[jcell - 1].g) { + ++gid; + } + idx_t icell = cell_sort[jcell].i; + glb_idx_gathered[icell] = gid; + } + + // 3) Scatter renumbered back + ATLAS_TRACE_MPI(SCATTER) { + mpi::comm().scatterv(glb_idx_gathered.data(), recvcounts.data(), recvdispls.data(), glb_idx.data(), + static_cast(glb_idx.size()), root); + } + + for (int jcell = 0; jcell < nb_cells; ++jcell) { + cells_glb_idx(cells_to_edit[jcell]) = glb_idx[jcell]; + } + + // nodes_glb_idx.dump( Log::info() ); + // Log::info() << std::endl; + cells.global_index().metadata().set("human_readable", true); +} + +// ------------------------------------------------------------------------------------- + +using uid_t = gidx_t; + +// ------------------------------------------------------------------ +class BuildHaloDistanceHelper; + +void increase_halo(Mesh& mesh); +void increase_halo_interior(BuildHaloDistanceHelper&); + +class EastWest : public util::PeriodicTransform { +public: + EastWest() { x_translation_ = -360.; } +}; + +class WestEast : public util::PeriodicTransform { +public: + WestEast() { x_translation_ = 360.; } +}; + +using Node2Elem = std::vector>; + +void build_lookup_node2elem(const Mesh& mesh, Node2Elem& node2elem) { + ATLAS_TRACE(); + + const mesh::Nodes& nodes = mesh.nodes(); + + node2elem.resize(nodes.size()); + for (idx_t jnode = 0; jnode < nodes.size(); ++jnode) { + node2elem[jnode].clear(); + node2elem[jnode].reserve(12); + } + + const mesh::HybridElements::Connectivity& elem_nodes = mesh.cells().node_connectivity(); + auto field_flags = array::make_view(mesh.cells().flags()); + auto patched = [&field_flags](idx_t e) { + using Topology = atlas::mesh::Nodes::Topology; + return Topology::check(field_flags(e), Topology::PATCH); + }; + + idx_t nb_elems = mesh.cells().size(); + for (idx_t elem = 0; elem < nb_elems; ++elem) { + if (not patched(elem)) { + for (idx_t n = 0; n < elem_nodes.cols(elem); ++n) { + int node = elem_nodes(elem, n); + node2elem[node].push_back(elem); + } + } + } +} + +void accumulate_partition_bdry_nodes_old(Mesh& mesh, std::vector& bdry_nodes) { + ATLAS_TRACE(); + + std::set bdry_nodes_set; + + std::vector facet_nodes; + std::vector connectivity_facet_to_elem; + + facet_nodes.reserve(mesh.nodes().size() * 4); + connectivity_facet_to_elem.reserve(facet_nodes.capacity() * 2); + + idx_t nb_facets(0); + idx_t nb_inner_facets(0); + idx_t missing_value; + accumulate_facets( + /*in*/ mesh.cells(), + /*in*/ mesh.nodes(), + /*out*/ facet_nodes, // shape(nb_facets,nb_nodes_per_facet) + /*out*/ connectivity_facet_to_elem, + /*out*/ nb_facets, + /*out*/ nb_inner_facets, + /*out*/ missing_value); + + for (idx_t jface = 0; jface < nb_facets; ++jface) { + if (connectivity_facet_to_elem[jface * 2 + 1] == missing_value) { + for (idx_t jnode = 0; jnode < 2; ++jnode) // 2 nodes per face + { + bdry_nodes_set.insert(facet_nodes[jface * 2 + jnode]); + } + } + } + bdry_nodes = std::vector(bdry_nodes_set.begin(), bdry_nodes_set.end()); +} + +void accumulate_partition_bdry_nodes(Mesh& mesh, idx_t halo, std::vector& bdry_nodes) { +#ifndef ATLAS_103 + /* deprecated */ + accumulate_partition_bdry_nodes_old(mesh, bdry_nodes); +#else + ATLAS_TRACE(); + const Mesh::Polygon& polygon = mesh.polygon(halo); + bdry_nodes = std::vector(polygon.begin(), polygon.end()); +#endif + +#ifdef ATLAS_103_SORT + /* not required */ + std::sort(bdry_nodes.begin(), bdry_nodes.end()); +#endif +} + +template +std::vector filter_nodes(std::vector nodes, const Predicate& predicate) { + std::vector filtered; + filtered.reserve(nodes.size()); + for (idx_t inode : nodes) { + if (predicate(inode)) { + filtered.push_back(inode); + } + } + return filtered; +} + +class Notification { +public: + void add_error(const std::string& note, const eckit::CodeLocation& loc) { + notes.push_back(note + " @ " + std::string(loc)); + } + void add_error(const std::string& note) { notes.push_back(note); } + + bool error() const { return notes.size() > 0; } + void reset() { notes.clear(); } + + std::string str() const { + std::stringstream stream; + for (size_t jnote = 0; jnote < notes.size(); ++jnote) { + if (jnote > 0) { + stream << "\n"; + } + stream << notes[jnote]; + } + return stream.str(); + } + + operator std::string() const { return str(); } + +private: + friend std::ostream& operator<<(std::ostream& s, const Notification& notes) { + s << notes.str(); + return s; + } + +private: + std::vector notes; +}; + +using Uid2Node = std::unordered_map; + + +void build_lookup_uid2node(Mesh& mesh, Uid2Node& uid2node) { + ATLAS_TRACE(); + Notification notes; + mesh::Nodes& nodes = mesh.nodes(); + array::ArrayView xy = array::make_view(nodes.xy()); + array::ArrayView glb_idx = array::make_view(nodes.global_index()); + idx_t nb_nodes = nodes.size(); + + UniqueLonLat compute_uid(mesh); + + uid2node.clear(); + for (idx_t jnode = 0; jnode < nb_nodes; ++jnode) { + uid_t uid = compute_uid(jnode); + bool inserted = uid2node.insert(std::make_pair(uid, jnode)).second; + if (not inserted) { + int other = uid2node[uid]; + std::stringstream msg; + msg << std::setprecision(10) << std::fixed << "Node uid: " << uid << " " << glb_idx(jnode) << " xy(" + << xy(jnode, XX) << "," << xy(jnode, YY) << ")"; + if (nodes.has_field("ij")) { + auto ij = array::make_view(nodes.field("ij")); + msg << " ij(" << ij(jnode, XX) << "," << ij(jnode, YY) << ")"; + } + msg << " has already been added as node " << glb_idx(other) << " (" << xy(other, XX) << "," << xy(other, YY) + << ")"; + if (nodes.has_field("ij")) { + auto ij = array::make_view(nodes.field("ij")); + msg << " ij(" << ij(other, XX) << "," << ij(other, YY) << ")"; + } + notes.add_error(msg.str()); + } + } + if (notes.error()) { + throw_Exception(notes.str(), Here()); + } +} + +/// For given nodes, find new connected elements +void accumulate_elements(const Mesh& mesh, const mpi::BufferView& request_node_uid, const Uid2Node& uid2node, + const Node2Elem& node2elem, std::vector& found_elements, + std::set& new_nodes_uid) { + + ATLAS_TRACE(); + const mesh::HybridElements::Connectivity& elem_nodes = mesh.cells().node_connectivity(); + const auto elem_part = array::make_view(mesh.cells().partition()); + + const idx_t nb_nodes = mesh.nodes().size(); + const idx_t nb_request_nodes = static_cast(request_node_uid.size()); + const int mpi_rank = static_cast(mpi::rank()); + + std::unordered_set found_elements_set; + found_elements_set.reserve(nb_request_nodes*2); + + for (idx_t jnode = 0; jnode < nb_request_nodes; ++jnode) { + uid_t uid = request_node_uid(jnode); + + idx_t inode = -1; + // search and get node index for uid + auto found = uid2node.find(uid); + if (found != uid2node.end()) { + inode = found->second; + } + if (inode != -1 && inode < nb_nodes) { + for (const idx_t e : node2elem[inode]) { + if (elem_part(e) == mpi_rank) { + found_elements_set.insert(e); + } + } + } + } + + // found_bdry_elements_set now contains elements for the nodes + found_elements = std::vector(found_elements_set.begin(), found_elements_set.end()); + + UniqueLonLat compute_uid(mesh); + + // Collect all nodes + new_nodes_uid.clear(); + for (const idx_t e : found_elements) { + idx_t nb_elem_nodes = elem_nodes.cols(e); + for (idx_t n = 0; n < nb_elem_nodes; ++n) { + new_nodes_uid.insert(compute_uid(elem_nodes(e, n))); + } + } + + // Remove nodes we already have in the request-buffer + for (idx_t jnode = 0; jnode < nb_request_nodes; ++jnode) { + new_nodes_uid.erase(request_node_uid(jnode)); + } +} + +class BuildHaloDistanceHelper { +public: + struct Buffers { + std::vector> node_part; + + std::vector> node_ridx; + + std::vector> node_flags; + + std::vector> node_glb_idx; + + std::vector> node_xy; + + std::vector> elem_glb_idx; + + std::vector> elem_nodes_id; + + std::vector> elem_nodes_displs; + + std::vector> elem_part; + + std::vector> elem_ridx; + + std::vector> elem_flags; + + std::vector> elem_type; + + Buffers(Mesh&) { + const idx_t mpi_size = static_cast(mpi::size()); + + node_part.resize(mpi_size); + node_ridx.resize(mpi_size); + node_flags.resize(mpi_size); + node_glb_idx.resize(mpi_size); + node_xy.resize(mpi_size); + elem_glb_idx.resize(mpi_size); + elem_nodes_id.resize(mpi_size); + elem_nodes_displs.resize(mpi_size); + elem_part.resize(mpi_size); + elem_ridx.resize(mpi_size); + elem_flags.resize(mpi_size); + elem_type.resize(mpi_size); + } + + void print(std::ostream& os) const { + const idx_t mpi_size = static_cast(mpi::size()); + os << "Nodes\n" + << "-----\n"; + idx_t n(0); + for (idx_t jpart = 0; jpart < mpi_size; ++jpart) { + for (idx_t jnode = 0; jnode < node_glb_idx[jpart].size(); ++jnode ) { + auto g = node_glb_idx[jpart][jnode]; + auto p = node_part[jpart][jnode]; + os << std::setw(4) << n++ << " : [ p" << p << " ] " << g << "\n"; + } + } + os << std::flush; + os << "Cells\n" + << "-----\n"; + idx_t e(0); + for (idx_t jpart = 0; jpart < mpi_size; ++jpart) { + const idx_t nb_elem = static_cast(elem_glb_idx[jpart].size()); + for (idx_t jelem = 0; jelem < nb_elem; ++jelem) { + os << std::setw(4) << e++ << " : [ t" << elem_type[jpart][jelem] << " -- p" + << elem_part[jpart][jelem] << "] " << elem_glb_idx[jpart][jelem] << "\n"; + } + } + os << std::flush; + } + friend std::ostream& operator<<(std::ostream& out, const Buffers& b) { + b.print(out); + return out; + } + }; + + static void all_to_all(Buffers& send, Buffers& recv) { + ATLAS_TRACE(); + const eckit::mpi::Comm& comm = mpi::comm(); + + ATLAS_TRACE_MPI(ALLTOALL) { + comm.allToAll(send.node_glb_idx, recv.node_glb_idx); + comm.allToAll(send.node_part, recv.node_part); + comm.allToAll(send.node_ridx, recv.node_ridx); + comm.allToAll(send.node_flags, recv.node_flags); + comm.allToAll(send.node_xy, recv.node_xy); + comm.allToAll(send.elem_glb_idx, recv.elem_glb_idx); + comm.allToAll(send.elem_nodes_id, recv.elem_nodes_id); + comm.allToAll(send.elem_part, recv.elem_part); + comm.allToAll(send.elem_ridx, recv.elem_ridx); + comm.allToAll(send.elem_type, recv.elem_type); + comm.allToAll(send.elem_flags, recv.elem_flags); + comm.allToAll(send.elem_nodes_displs, recv.elem_nodes_displs); + } + } + + struct Status { + std::vector new_periodic_ghost_points; + std::vector> new_periodic_ghost_cells; + } status; + +public: + BuildHaloDistance& builder_; + Mesh& mesh; + array::ArrayView xy; + array::ArrayView lonlat; + array::ArrayView glb_idx; + array::ArrayView part; + array::IndexView ridx; + array::ArrayView flags; + array::ArrayView halo; + array::ArrayView ghost; + mesh::HybridElements::Connectivity* elem_nodes; + array::ArrayView elem_part; + array::IndexView elem_ridx; + array::ArrayView elem_flags; + array::ArrayView elem_glb_idx; + + std::vector bdry_nodes; + Node2Elem node_to_elem; + Uid2Node uid2node; + UniqueLonLat compute_uid; + idx_t halosize; + +public: + BuildHaloDistanceHelper(BuildHaloDistance& builder, Mesh& _mesh): + builder_(builder), + mesh(_mesh), + xy(array::make_view(mesh.nodes().xy())), + lonlat(array::make_view(mesh.nodes().lonlat())), + glb_idx(array::make_view(mesh.nodes().global_index())), + part(array::make_view(mesh.nodes().partition())), + ridx(array::make_indexview(mesh.nodes().remote_index())), + flags(array::make_view(mesh.nodes().flags())), + halo(array::make_view(mesh.nodes().halo())), + ghost(array::make_view(mesh.nodes().ghost())), + elem_nodes(&mesh.cells().node_connectivity()), + elem_part(array::make_view(mesh.cells().partition())), + elem_ridx(array::make_indexview(mesh.cells().remote_index())), + elem_flags(array::make_view(mesh.cells().flags())), + elem_glb_idx(array::make_view(mesh.cells().global_index())), + compute_uid(mesh) { + halosize = 0; + mesh.metadata().get("halo", halosize); + // update(); + } + + void update() { + compute_uid.update(); + mesh::Nodes& nodes = mesh.nodes(); + xy = array::make_view(nodes.xy()); + lonlat = array::make_view(nodes.lonlat()); + glb_idx = array::make_view(nodes.global_index()); + part = array::make_view(nodes.partition()); + ridx = array::make_indexview(nodes.remote_index()); + flags = array::make_view(nodes.flags()); + ghost = array::make_view(nodes.ghost()); + halo = array::make_view(nodes.halo()); + + elem_nodes = &mesh.cells().node_connectivity(); + elem_part = array::make_view(mesh.cells().partition()); + elem_ridx = array::make_indexview(mesh.cells().remote_index()); + elem_flags = array::make_view(mesh.cells().flags()); + elem_glb_idx = array::make_view(mesh.cells().global_index()); + } + + template + void fill_sendbuffer(Buffers& buf, const NodeContainer& nodes_uid, const ElementContainer& elems, const int p) { + ATLAS_TRACE(); + + idx_t nb_nodes = static_cast(nodes_uid.size()); + buf.node_glb_idx[p].resize(nb_nodes); + buf.node_part[p].resize(nb_nodes); + buf.node_ridx[p].resize(nb_nodes); + buf.node_flags[p].resize(nb_nodes, Topology::NONE); + buf.node_xy[p].resize(2 * nb_nodes); + + idx_t jnode = 0; + typename NodeContainer::const_iterator it; + for (it = nodes_uid.begin(); it != nodes_uid.end(); ++it, ++jnode) { + uid_t uid = *it; + auto found = uid2node.find(uid); + if (found != uid2node.end()) // Point exists inside domain + { + idx_t node = found->second; + buf.node_glb_idx[p][jnode] = glb_idx(node); + buf.node_part[p][jnode] = part(node); + buf.node_ridx[p][jnode] = ridx(node); + buf.node_xy[p][jnode * 2 + XX] = xy(node, XX); + buf.node_xy[p][jnode * 2 + YY] = xy(node, YY); + Topology::set(buf.node_flags[p][jnode], flags(node) | Topology::GHOST); + } + else { + Log::warning() << "Node with uid " << uid << " needed by [" << p << "] was not found in [" + << mpi::rank() << "]." << std::endl; + ATLAS_ASSERT(false); + } + } + + idx_t nb_elems = static_cast(elems.size()); + + idx_t nb_elem_nodes(0); + for (idx_t jelem = 0; jelem < nb_elems; ++jelem) { + idx_t ielem = elems[jelem]; + nb_elem_nodes += elem_nodes->cols(ielem); + } + + buf.elem_glb_idx[p].resize(nb_elems); + buf.elem_part[p].resize(nb_elems); + buf.elem_ridx[p].resize(nb_elems); + buf.elem_flags[p].resize(nb_elems, Topology::NONE); + buf.elem_type[p].resize(nb_elems); + buf.elem_nodes_id[p].resize(nb_elem_nodes); + buf.elem_nodes_displs[p].resize(nb_elems); + idx_t jelemnode(0); + for (idx_t jelem = 0; jelem < nb_elems; ++jelem) { + buf.elem_nodes_displs[p][jelem] = jelemnode; + idx_t ielem = elems[jelem]; + + buf.elem_glb_idx[p][jelem] = elem_glb_idx(ielem); + buf.elem_part[p][jelem] = elem_part(ielem); + buf.elem_ridx[p][jelem] = elem_ridx(ielem); + Topology::set(buf.elem_flags[p][jelem], elem_flags(ielem)); + buf.elem_type[p][jelem] = mesh.cells().type_idx(ielem); + for (idx_t jnode = 0; jnode < elem_nodes->cols(ielem); ++jnode) { + buf.elem_nodes_id[p][jelemnode++] = compute_uid((*elem_nodes)(ielem, jnode)); + } + Topology::set(buf.elem_flags[p][jelem], Topology::GHOST); + } + } + + template + void fill_sendbuffer(Buffers& buf, const NodeContainer& nodes_uid, const ElementContainer& elems, + const util::PeriodicTransform& transform, int newflags, const int p) { + // ATLAS_TRACE(); + + idx_t nb_nodes = static_cast(nodes_uid.size()); + buf.node_glb_idx[p].resize(nb_nodes); + buf.node_part[p].resize(nb_nodes); + buf.node_ridx[p].resize(nb_nodes); + buf.node_flags[p].resize(nb_nodes, Topology::NONE); + buf.node_xy[p].resize(2 * nb_nodes); + + int jnode = 0; + typename NodeContainer::const_iterator it; + for (it = nodes_uid.begin(); it != nodes_uid.end(); ++it, ++jnode) { + uid_t uid = *it; + + auto found = uid2node.find(uid); + if (found != uid2node.end()) // Point exists inside domain + { + int node = found->second; + buf.node_part[p][jnode] = part(node); + buf.node_ridx[p][jnode] = ridx(node); + buf.node_xy[p][jnode * 2 + XX] = xy(node, XX); + buf.node_xy[p][jnode * 2 + YY] = xy(node, YY); + transform(&buf.node_xy[p][jnode * 2], -1); + // Global index of node is based on UID of destination + buf.node_glb_idx[p][jnode] = util::unique_lonlat(&buf.node_xy[p][jnode * 2]); + Topology::set(buf.node_flags[p][jnode], newflags); + } + else { + Log::warning() << "Node with uid " << uid << " needed by [" << p << "] was not found in [" + << mpi::rank() << "]." << std::endl; + ATLAS_ASSERT(false); + } + } + + idx_t nb_elems = static_cast(elems.size()); + + idx_t nb_elem_nodes(0); + for (idx_t jelem = 0; jelem < nb_elems; ++jelem) { + idx_t ielem = elems[jelem]; + nb_elem_nodes += elem_nodes->cols(ielem); + } + + buf.elem_glb_idx[p].resize(nb_elems); + buf.elem_part[p].resize(nb_elems); + buf.elem_ridx[p].resize(nb_elems); + buf.elem_flags[p].resize(nb_elems, Topology::NONE); + buf.elem_type[p].resize(nb_elems); + buf.elem_nodes_id[p].resize(nb_elem_nodes); + buf.elem_nodes_displs[p].resize(nb_elems); + idx_t jelemnode(0); + for (idx_t jelem = 0; jelem < nb_elems; ++jelem) { + buf.elem_nodes_displs[p][jelem] = jelemnode; + idx_t ielem = elems[jelem]; + buf.elem_part[p][jelem] = elem_part(ielem); + buf.elem_ridx[p][jelem] = elem_ridx(ielem); + Topology::set(buf.elem_flags[p][jelem], elem_flags(ielem) | newflags); + buf.elem_type[p][jelem] = mesh.cells().type_idx(ielem); + std::vector crds(elem_nodes->cols(ielem) * 2); + for (idx_t jnode = 0; jnode < elem_nodes->cols(ielem); ++jnode) { + double crd[] = {xy((*elem_nodes)(ielem, jnode), XX), xy((*elem_nodes)(ielem, jnode), YY)}; + transform(crd, -1); + buf.elem_nodes_id[p][jelemnode++] = util::unique_lonlat(crd); + crds[jnode * 2 + XX] = crd[XX]; + crds[jnode * 2 + YY] = crd[YY]; + } + // Global index of element is based on UID of destination + + buf.elem_glb_idx[p][jelem] = -util::unique_lonlat(crds.data(), elem_nodes->cols(ielem)); + } + } + + void add_nodes(Buffers& buf) { + ATLAS_TRACE(); + + const idx_t mpi_size = static_cast(mpi::size()); + + mesh::Nodes& nodes = mesh.nodes(); + int nb_nodes = nodes.size(); + + // Nodes might be duplicated from different Tasks. We need to identify + // unique entries + std::vector node_uid(nb_nodes); + std::set new_node_uid; + { + ATLAS_TRACE("compute node_uid"); + for (int jnode = 0; jnode < nb_nodes; ++jnode) { + node_uid[jnode] = compute_uid(jnode); + } + std::sort(node_uid.begin(), node_uid.end()); + } + auto node_already_exists = [&node_uid, &new_node_uid](uid_t uid) { + std::vector::iterator it = std::lower_bound(node_uid.begin(), node_uid.end(), uid); + bool not_found = (it == node_uid.end() || uid < *it); + if (not_found) { + bool inserted = new_node_uid.insert(uid).second; + return not inserted; + } + else { + return true; + } + }; + + std::vector> rfn_idx(mpi_size); + for (idx_t jpart = 0; jpart < mpi_size; ++jpart) { + rfn_idx[jpart].reserve(buf.node_glb_idx[jpart].size()); + } + + int nb_new_nodes = 0; + for (idx_t jpart = 0; jpart < mpi_size; ++jpart) { + const idx_t nb_nodes_on_part = static_cast(buf.node_glb_idx[jpart].size()); + for (idx_t n = 0; n < nb_nodes_on_part; ++n) { + std::array crd{buf.node_xy[jpart][n * 2 + XX], buf.node_xy[jpart][n * 2 + YY]}; + mesh.projection().xy2lonlat(crd.data()); + + if (not node_already_exists(util::unique_lonlat(crd))) { + rfn_idx[jpart].push_back(n); + } + } + nb_new_nodes += rfn_idx[jpart].size(); + } + + // Resize nodes + // ------------ + nodes.resize(nb_nodes + nb_new_nodes); + flags = array::make_view(nodes.flags()); + halo = array::make_view(nodes.halo()); + glb_idx = array::make_view(nodes.global_index()); + part = array::make_view(nodes.partition()); + ridx = array::make_indexview(nodes.remote_index()); + xy = array::make_view(nodes.xy()); + lonlat = array::make_view(nodes.lonlat()); + ghost = array::make_view(nodes.ghost()); + + compute_uid.update(); + + // Add new nodes + // ------------- + int new_node = 0; + for (idx_t jpart = 0; jpart < mpi_size; ++jpart) { + for (size_t n = 0; n < rfn_idx[jpart].size(); ++n) { + int loc_idx = nb_nodes + new_node; + halo(loc_idx) = halosize + 1; + Topology::reset(flags(loc_idx), buf.node_flags[jpart][rfn_idx[jpart][n]]); + ghost(loc_idx) = Topology::check(flags(loc_idx), Topology::GHOST); + glb_idx(loc_idx) = buf.node_glb_idx[jpart][rfn_idx[jpart][n]]; + part(loc_idx) = buf.node_part[jpart][rfn_idx[jpart][n]]; + ridx(loc_idx) = buf.node_ridx[jpart][rfn_idx[jpart][n]]; + PointXY pxy(&buf.node_xy[jpart][rfn_idx[jpart][n] * 2]); + xy(loc_idx, XX) = pxy.x(); + xy(loc_idx, YY) = pxy.y(); + PointLonLat pll = mesh.projection().lonlat(pxy); + lonlat(loc_idx, XX) = pll.lon(); + lonlat(loc_idx, YY) = pll.lat(); + + if (Topology::check(flags(loc_idx), Topology::PERIODIC) and + not Topology::check(flags(loc_idx), Topology::BC)) { + status.new_periodic_ghost_points.push_back(loc_idx); + } + + // make sure new node was not already there + { + uid_t uid = compute_uid(loc_idx); + auto found = uid2node.find(uid); + if (found != uid2node.end()) { + int other = found->second; + std::stringstream msg; + msg << "New node loc " << loc_idx << " with uid " << uid << ":\n" + << glb_idx(loc_idx) << "(" << xy(loc_idx, XX) << "," << xy(loc_idx, YY) << ")\n"; + msg << "Existing already loc " << other << " : " << glb_idx(other) << "(" << xy(other, XX) + << "," << xy(other, YY) << ")\n"; + throw_Exception(msg.str(), Here()); + } + uid2node[uid] = nb_nodes + new_node; + } + ++new_node; + } + } + } + + void add_elements(Buffers& buf) { + ATLAS_TRACE(); + + const idx_t mpi_size = static_cast(mpi::size()); + auto cell_gidx = array::make_view(mesh.cells().global_index()); + // Elements might be duplicated from different Tasks. We need to identify + // unique entries + int nb_elems = mesh.cells().size(); + // std::set elem_uid; + std::vector elem_uid(2 * nb_elems); + std::set new_elem_uid; + { + ATLAS_TRACE("compute elem_uid"); + atlas_omp_parallel_for (int jelem = 0; jelem < nb_elems; ++jelem) { + elem_uid[jelem * 2 + 0] = -compute_uid(elem_nodes->row(jelem)); + elem_uid[jelem * 2 + 1] = cell_gidx(jelem); + } + omp::sort(elem_uid.begin(), elem_uid.end()); + } + auto element_already_exists = [&elem_uid, &new_elem_uid](uid_t uid) -> bool { + std::vector::iterator it = std::lower_bound(elem_uid.begin(), elem_uid.end(), uid); + bool not_found = (it == elem_uid.end() || uid < *it); + if (not_found) { + bool inserted = new_elem_uid.insert(uid).second; + return not inserted; + } + else { + return true; + } + }; + + if (not status.new_periodic_ghost_cells.size()) { + status.new_periodic_ghost_cells.resize(mesh.cells().nb_types()); + } + + std::vector> received_new_elems(mpi_size); + for (idx_t jpart = 0; jpart < mpi_size; ++jpart) { + received_new_elems[jpart].reserve(buf.elem_glb_idx[jpart].size()); + } + + idx_t nb_new_elems(0); + for (idx_t jpart = 0; jpart < mpi_size; ++jpart) { + const idx_t nb_elems_from_part = static_cast(buf.elem_glb_idx[jpart].size()); + for (idx_t e = 0; e < nb_elems_from_part; ++e) { + if (element_already_exists(buf.elem_glb_idx[jpart][e]) == false) { + received_new_elems[jpart].emplace_back(e); + } + } + nb_new_elems += received_new_elems[jpart].size(); + } + + std::vector>> elements_of_type(mesh.cells().nb_types(), + std::vector>(mpi_size)); + std::vector nb_elements_of_type(mesh.cells().nb_types(), 0); + + for (idx_t jpart = 0; jpart < mpi_size; ++jpart) { + for (const idx_t ielem : received_new_elems[jpart]) { + elements_of_type[buf.elem_type[jpart][ielem]][jpart].emplace_back(ielem); + ++nb_elements_of_type[buf.elem_type[jpart][ielem]]; + } + } + + for (idx_t t = 0; t < mesh.cells().nb_types(); ++t) { + const std::vector>& elems = elements_of_type[t]; + mesh::Elements& elements = mesh.cells().elements(t); + + // Add new elements + BlockConnectivity& node_connectivity = elements.node_connectivity(); + if (nb_elements_of_type[t] == 0) { + continue; + } + + idx_t old_size = elements.size(); + idx_t new_elems_pos = elements.add(nb_elements_of_type[t]); + + auto elem_type_glb_idx = elements.view(mesh.cells().global_index()); + auto elem_type_part = elements.view(mesh.cells().partition()); + auto elem_type_ridx = elements.indexview(mesh.cells().remote_index()); + auto elem_type_halo = elements.view(mesh.cells().halo()); + auto elem_type_flags = elements.view(mesh.cells().flags()); + + // Copy information in new elements + idx_t new_elem(0); + for (idx_t jpart = 0; jpart < mpi_size; ++jpart) { + for (const idx_t jelem : elems[jpart]) { + int loc_idx = new_elems_pos + new_elem; + elem_type_glb_idx(loc_idx) = std::abs(buf.elem_glb_idx[jpart][jelem]); + elem_type_part(loc_idx) = buf.elem_part[jpart][jelem]; + elem_type_ridx(loc_idx) = buf.elem_ridx[jpart][jelem]; + elem_type_halo(loc_idx) = halosize + 1; + elem_type_flags(loc_idx) = buf.elem_flags[jpart][jelem]; + for (idx_t n = 0; n < node_connectivity.cols(); ++n) { + node_connectivity.set( + loc_idx, n, uid2node[buf.elem_nodes_id[jpart][buf.elem_nodes_displs[jpart][jelem] + n]]); + } + + if (Topology::check(elem_type_flags(loc_idx), Topology::PERIODIC)) { + status.new_periodic_ghost_cells[t].push_back(old_size + new_elem); + } + ++new_elem; + } + } + } + } + + std::vector> gather_nb_elements() { + idx_t nb_types = mesh.cells().nb_types(); + idx_t mpi_size = mpi::comm().size(); + std::vector elems_per_type(nb_types); + for (idx_t t = 0; t < nb_types; ++t) { + elems_per_type[t] = mesh.cells().elements(t).size(); + } + atlas::mpi::Buffer recv(mpi_size); + ATLAS_TRACE_MPI(ALLGATHER) { mpi::comm().allGatherv(elems_per_type.begin(), elems_per_type.end(), recv); } + + std::vector> result(nb_types, std::vector(mpi_size)); + for (idx_t p = 0; p < mpi_size; ++p) { + auto recv_p = recv[p]; + for (idx_t t = 0; t < nb_types; ++t) { + result[t][p] = recv_p[t]; + } + } + return result; + }; + + + void add_buffers(Buffers& buf) { + add_nodes(buf); + + auto nb_elements_pre = gather_nb_elements(); + + add_elements(buf); + + auto nb_elements_post = gather_nb_elements(); + + // Renumber remote_index as the number of elements per type might have grown + { + auto nb_partitions = mpi::comm().size(); + auto nb_types = mesh.cells().nb_types(); + + std::vector> accumulated_diff(nb_types, std::vector(nb_partitions)); + for (idx_t t = 1; t < nb_types; ++t) { + for (idx_t p = 0; p < nb_partitions; ++p) { + accumulated_diff[t][p] = + accumulated_diff[t - 1][p] + nb_elements_post[t - 1][p] - nb_elements_pre[t - 1][p]; + } + } + + for (idx_t t = 0; t < nb_types; ++t) { + auto& elements = mesh.cells().elements(t); + auto partition = elements.view(mesh.cells().partition()); + auto ridx = elements.indexview(mesh.cells().remote_index()); + + auto& accumulated_diff_t = accumulated_diff[t]; + + for (idx_t elem = 0; elem < elements.size(); ++elem) { + ridx(elem) += accumulated_diff_t[partition(elem)]; + } + } + } + + update(); + } +}; +} + +namespace { +void gather_bdry_nodes(const BuildHaloDistanceHelper& helper, const std::vector& send, + atlas::mpi::Buffer& recv, bool periodic = false) { + auto& comm = mpi::comm(); +#ifndef ATLAS_103 + /* deprecated */ + ATLAS_TRACE("gather_bdry_nodes old way"); + { + ATLAS_TRACE_MPI(ALLGATHER) { comm.allGatherv(send.begin(), send.end(), recv); } + } +#else + ATLAS_TRACE(); + + Mesh::PartitionGraph::Neighbours neighbours = helper.mesh.nearestNeighbourPartitions(); + if (periodic) { + // add own rank to neighbours to allow periodicity with self (pole caps) + idx_t rank = comm.rank(); + neighbours.insert(std::upper_bound(neighbours.begin(), neighbours.end(), rank), rank); + } + + const idx_t mpi_size = comm.size(); + const int counts_tag = 0; + const int buffer_tag = 1; + + std::vector counts_requests; + counts_requests.reserve(neighbours.size()); + std::vector buffer_requests; + buffer_requests.reserve(neighbours.size()); + + int sendcnt = send.size(); + ATLAS_TRACE_MPI(ISEND) { + for (idx_t to : neighbours) { + counts_requests.push_back(comm.iSend(sendcnt, to, counts_tag)); + } + } + + recv.counts.assign(0, mpi_size); + + ATLAS_TRACE_MPI(IRECEIVE) { + for (idx_t from : neighbours) { + counts_requests.push_back(comm.iReceive(recv.counts[from], from, counts_tag)); + } + } + + ATLAS_TRACE_MPI(ISEND) { + for (idx_t to : neighbours) { + buffer_requests.push_back(comm.iSend(send.data(), send.size(), to, buffer_tag)); + } + } + + ATLAS_TRACE_MPI(WAIT) { + for (auto request : counts_requests) { + comm.wait(request); + } + } + + recv.displs[0] = 0; + recv.cnt = recv.counts[0]; + for (idx_t jpart = 1; jpart < mpi_size; ++jpart) { + recv.displs[jpart] = recv.displs[jpart - 1] + recv.counts[jpart - 1]; + recv.cnt += recv.counts[jpart]; + } + recv.buffer.resize(recv.cnt); + + ATLAS_TRACE_MPI(IRECEIVE) { + for (idx_t from : neighbours) { + buffer_requests.push_back( + comm.iReceive(recv.buffer.data() + recv.displs[from], recv.counts[from], from, buffer_tag)); + } + } + + ATLAS_TRACE_MPI(WAIT) { + for (auto request : buffer_requests) { + comm.wait(request); + } + } +#endif +} +} // namespace + +namespace { +void increase_halo_interior(BuildHaloDistanceHelper& helper) { + helper.update(); + if (helper.node_to_elem.size() == 0) { + build_lookup_node2elem(helper.mesh, helper.node_to_elem); + } + + if (helper.uid2node.size() == 0) { + build_lookup_uid2node(helper.mesh, helper.uid2node); + } + + // All buffers needed to move elements and nodes + BuildHaloDistanceHelper::Buffers sendmesh(helper.mesh); + BuildHaloDistanceHelper::Buffers recvmesh(helper.mesh); + + // 1) Find boundary nodes of this partition: + + accumulate_partition_bdry_nodes(helper.mesh, helper.halosize, helper.bdry_nodes); + const std::vector& bdry_nodes = helper.bdry_nodes; + const idx_t nb_bdry_nodes = static_cast(bdry_nodes.size()); + + // 2) Communicate uid of these boundary nodes to other partitions + + std::vector send_bdry_nodes_uid(bdry_nodes.size()); + for (idx_t jnode = 0; jnode < nb_bdry_nodes; ++jnode) { + send_bdry_nodes_uid[jnode] = helper.compute_uid(bdry_nodes[jnode]); + } + + idx_t mpi_size = mpi::size(); + atlas::mpi::Buffer recv_bdry_nodes_uid_from_parts(mpi_size); + + gather_bdry_nodes(helper, send_bdry_nodes_uid, recv_bdry_nodes_uid_from_parts); + + { + runtime::trace::Barriers set_barriers(false); + runtime::trace::Logging set_logging(false); +#ifndef ATLAS_103 + /* deprecated */ + atlas_omp_parallel_for (idx_t jpart = 0; jpart < mpi_size; ++jpart) +#else + const Mesh::PartitionGraph::Neighbours neighbours = helper.mesh.nearestNeighbourPartitions(); + for (idx_t jpart : neighbours) +#endif + { + // 3) Find elements and nodes completing these elements in + // other tasks that have my nodes through its UID + + mpi::BufferView recv_bdry_nodes_uid = recv_bdry_nodes_uid_from_parts[jpart]; + + std::vector found_bdry_elems; + std::set found_bdry_nodes_uid; + + accumulate_elements(helper.mesh, recv_bdry_nodes_uid, helper.uid2node, helper.node_to_elem, found_bdry_elems, + found_bdry_nodes_uid); + + // 4) Fill node and element buffers to send back + helper.fill_sendbuffer(sendmesh, found_bdry_nodes_uid, found_bdry_elems, jpart); + } + } + + // 5) Now communicate all buffers + helper.all_to_all(sendmesh, recvmesh); + +// 6) Adapt mesh +#ifdef DEBUG_OUTPUT + Log::debug() << "increase_halo_interior -- recv: \n" << recvmesh << std::endl; +#endif + helper.add_buffers(recvmesh); +} + +class PeriodicPoints { +public: + PeriodicPoints(Mesh& mesh, int flag, idx_t N): flags_(array::make_view(mesh.nodes().flags())) { + flag_ = flag; + N_ = N; + } + + bool operator()(int j) const { + if (j >= N_) { + return false; + } + if (Topology::check(flags_(j), flag_)) { + return true; + } + return false; + } + +private: + int N_; + int flag_; + array::ArrayView flags_; + + friend std::ostream& operator<<(std::ostream& os, const PeriodicPoints& periodic_points) { + os << "["; + for (idx_t j = 0; j < periodic_points.flags_.shape(0); ++j) { + if (periodic_points(j)) { + os << " " << j + 1; + } + } + os << " ]"; + return os; + } +}; + +void increase_halo_periodic(BuildHaloDistanceHelper& helper, const PeriodicPoints& periodic_points, + const util::PeriodicTransform& transform, int newflags) { + helper.update(); + // if (helper.node_to_elem.size() == 0 ) !!! NOT ALLOWED !!! (atlas_test_halo + // will fail) + build_lookup_node2elem(helper.mesh, helper.node_to_elem); + + // if( helper.uid2node.size() == 0 ) !!! NOT ALLOWED !!! (atlas_test_halo will + // fail) + build_lookup_uid2node(helper.mesh, helper.uid2node); + + // All buffers needed to move elements and nodes + BuildHaloDistanceHelper::Buffers sendmesh(helper.mesh); + BuildHaloDistanceHelper::Buffers recvmesh(helper.mesh); + + // 1) Find boundary nodes of this partition: + + if (!helper.bdry_nodes.size()) { + accumulate_partition_bdry_nodes(helper.mesh, helper.halosize, helper.bdry_nodes); + } + + std::vector bdry_nodes = filter_nodes(helper.bdry_nodes, periodic_points); + const idx_t nb_bdry_nodes = static_cast(bdry_nodes.size()); + + // 2) Compute transformed uid of these boundary nodes and send to other + // partitions + + std::vector send_bdry_nodes_uid(nb_bdry_nodes); + for (idx_t jnode = 0; jnode < nb_bdry_nodes; ++jnode) { + double crd[] = {helper.xy(bdry_nodes[jnode], XX), helper.xy(bdry_nodes[jnode], YY)}; + transform(crd, +1); + // Log::info() << " crd " << crd[0] << " " << crd[1] << " uid " << + // util::unique_lonlat(crd) << std::endl; + send_bdry_nodes_uid[jnode] = util::unique_lonlat(crd); + } + + idx_t mpi_size = mpi::size(); + atlas::mpi::Buffer recv_bdry_nodes_uid_from_parts(mpi_size); + + gather_bdry_nodes(helper, send_bdry_nodes_uid, recv_bdry_nodes_uid_from_parts, + /* periodic = */ true); + + { + runtime::trace::Barriers set_barriers(false); + runtime::trace::Logging set_logging(false); +#ifndef ATLAS_103 + /* deprecated */ + atlas_omp_parallel_for (idx_t jpart = 0; jpart < mpi_size; ++jpart) +#else + Mesh::PartitionGraph::Neighbours neighbours = helper.mesh.nearestNeighbourPartitions(); + // add own rank to neighbours to allow periodicity with self (pole caps) + idx_t rank = mpi::rank(); + neighbours.insert(std::upper_bound(neighbours.begin(), neighbours.end(), rank), rank); + for (idx_t jpart : neighbours) +#endif + { + + // 3) Find elements and nodes completing these elements in + // other tasks that have my nodes through its UID + + atlas::mpi::BufferView recv_bdry_nodes_uid = recv_bdry_nodes_uid_from_parts[jpart]; + + std::vector found_bdry_elems; + std::set found_bdry_nodes_uid; + + accumulate_elements(helper.mesh, recv_bdry_nodes_uid, helper.uid2node, helper.node_to_elem, found_bdry_elems, + found_bdry_nodes_uid); + + // 4) Fill node and element buffers to send back + helper.fill_sendbuffer(sendmesh, found_bdry_nodes_uid, found_bdry_elems, transform, newflags, jpart); + } + } + + // 5) Now communicate all buffers + helper.all_to_all(sendmesh, recvmesh); + +// 6) Adapt mesh +#ifdef DEBUG_OUTPUT + Log::debug() << "increase_halo_periodic -- recv: \n" << recvmesh << std::endl; +#endif + helper.add_buffers(recvmesh); +} + +} + +BuildHaloDistance::BuildHaloDistance(Mesh& mesh, const eckit::Configuration& config): mesh_(mesh), periodic_cells_local_index_(mesh.cells().nb_types()) {} + +void BuildHaloDistance::operator()(double nb_elems) { + ATLAS_TRACE("BuildHaloDistance"); + + mpi::Scope mpi_scope(mesh_.mpi_comm()); + + int halo = 0; + mesh_.metadata().get("halo", halo); + + // Locked halos must be set at mesh generation. + bool haloLocked = false; + if (mesh_.metadata().get("halo_locked", haloLocked)) { + if (haloLocked && halo < nb_elems) { + const auto errMsg = + "Error: Halo size locked. Please set MeshGenerator " + "parameter \"halo\" to at least " + + std::to_string(nb_elems) + + ".\n" + "Current halo size is " + + std::to_string(halo) + ".\n"; + throw_Exception(errMsg, Here()); + } + } + + if (halo == nb_elems) { + return; + } + + ATLAS_TRACE("Increasing mesh halo"); + + for (int jhalo = halo; jhalo < nb_elems; ++jhalo) { + Log::debug() << "Increase halo " << jhalo + 1 << std::endl; + idx_t nb_nodes_before_halo_increase = mesh_.nodes().size(); + + BuildHaloDistanceHelper helper(*this, mesh_); + + ATLAS_TRACE_SCOPE("increase_halo_interior") { increase_halo_interior(helper); } + + PeriodicPoints westpts(mesh_, Topology::PERIODIC | Topology::WEST, nb_nodes_before_halo_increase); + +#ifdef DEBUG_OUTPUT + Log::debug() << " periodic west : " << westpts << std::endl; +#endif + ATLAS_TRACE_SCOPE("increase_halo_periodic West") { + increase_halo_periodic(helper, westpts, WestEast(), Topology::PERIODIC | Topology::WEST | Topology::GHOST); + } + + PeriodicPoints eastpts(mesh_, Topology::PERIODIC | Topology::EAST, nb_nodes_before_halo_increase); + +#ifdef DEBUG_OUTPUT + Log::debug() << " periodic east : " << eastpts << std::endl; +#endif + ATLAS_TRACE_SCOPE("increase_halo_periodic East") { + increase_halo_periodic(helper, eastpts, EastWest(), Topology::PERIODIC | Topology::EAST | Topology::GHOST); + } + + for (idx_t p : helper.status.new_periodic_ghost_points) { + periodic_points_local_index_.push_back(p); + } + for (int t = 0; t < mesh_.cells().nb_types(); ++t) { + for (idx_t p : helper.status.new_periodic_ghost_cells[t]) { + periodic_cells_local_index_[t].push_back(p); + } + } + + { + std::stringstream ss; + ss << "nb_nodes_including_halo[" << jhalo + 1 << "]"; + mesh_.metadata().set(ss.str(), mesh_.nodes().size()); + } + + for (int t = 0; t < mesh_.cells().nb_types(); ++t) { + std::stringstream ss; + ss << "nb_cells_including_halo[" << t << "][" << jhalo + 1 << "]"; + mesh_.metadata().set(ss.str(), mesh_.cells().elements(t).size()); + } + mesh_.metadata().set("halo", jhalo + 1); + mesh_.nodes().global_index().metadata().set("human_readable", false); + mesh_.cells().global_index().metadata().set("human_readable", false); + +#ifdef DEBUG_OUTPUT + output::Gmsh gmsh2d("build-halo-mesh2d.msh", util::Config("ghost", true)("coordinates", "xy")); + output::Gmsh gmsh3d("build-halo-mesh3d.msh", util::Config("ghost", true)("coordinates", "xyz")); + renumber_nodes_glb_idx(mesh_.nodes()); + BuildXYZField("xyz", true)(mesh_.nodes()); + mesh_.metadata().set("halo", jhalo + 1); + gmsh2d.write(mesh_); + gmsh3d.write(mesh_); +#endif + } + + make_nodes_global_index_human_readable(*this, mesh_.nodes(), + /*do_all*/ false); + + make_cells_global_index_human_readable(*this, mesh_.cells(), + /*do_all*/ false); + // renumber_nodes_glb_idx (mesh_.nodes()); +} + +void build_halo_distance(Mesh& mesh, double distance) { + BuildHaloDistance f(mesh, util::NoConfig()); + f(distance); +} + +// ------------------------------------------------------------------ +// C wrapper interfaces to C++ routines + +void atlas__build_halo_distance(Mesh::Implementation* mesh, double radius) { + ATLAS_ASSERT(mesh != nullptr, "Cannot access uninitialised atlas_Mesh"); + Mesh m(mesh); + build_halo_distance(m, radius); +} + +// ------------------------------------------------------------------ + +} // namespace actions +} // namespace mesh +} // namespace atlas diff --git a/src/atlas/mesh/actions/BuildHaloDistance.h b/src/atlas/mesh/actions/BuildHaloDistance.h new file mode 100644 index 000000000..dde24e55f --- /dev/null +++ b/src/atlas/mesh/actions/BuildHaloDistance.h @@ -0,0 +1,69 @@ +/* + * (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. + */ + +#pragma once + +#include +#include + +#include "atlas/library/config.h" +#include "atlas/util/Config.h" + +namespace atlas { +class Mesh; + +namespace mesh { +namespace actions { + +class BuildHaloDistance { +public: + BuildHaloDistance(Mesh& mesh, const eckit::Configuration&); + void operator()(double distance); + +private: + Mesh& mesh_; + +public: + std::vector periodic_points_local_index_; + std::vector> periodic_cells_local_index_; +}; + +/// @brief Enlarge each partition of the mesh with a halo of elements +/// @param [inout] mesh The mesh to enlarge +/// @param [in] radius Size of the halo +/// @author Willem Deconinck +/// @date October 2023 +void build_halo_distance(Mesh& mesh, double distance); + +} // namespace actions +} // namespace mesh +} // namespace atlas + +// ------------------------------------------------------------------ +// C wrapper interfaces to C++ routines + +namespace atlas { +namespace mesh { +namespace detail { +class MeshImpl; +} +namespace actions { + +// ------------------------------------------------------------------ + +extern "C" { +void atlas__build_halo_distance(detail::MeshImpl* mesh, double distance); +} + +// ------------------------------------------------------------------ + +} // namespace actions +} // namespace mesh +} // namespace atlas diff --git a/src/atlas/mesh/detail/MeshBuilderIntf.cc b/src/atlas/mesh/detail/MeshBuilderIntf.cc new file mode 100644 index 000000000..3d04fa75b --- /dev/null +++ b/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 diff --git a/src/atlas/mesh/detail/MeshBuilderIntf.h b/src/atlas/mesh/detail/MeshBuilderIntf.h new file mode 100644 index 000000000..4b5d5eb11 --- /dev/null +++ b/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 diff --git a/src/atlas_f/CMakeLists.txt b/src/atlas_f/CMakeLists.txt index f01da0392..7075be1bb 100644 --- a/src/atlas_f/CMakeLists.txt +++ b/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 diff --git a/src/atlas_f/atlas_module.F90 b/src/atlas_f/atlas_module.F90 index 4e184c8a6..d82553e86 100644 --- a/src/atlas_f/atlas_module.F90 +++ b/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, & diff --git a/src/atlas_f/mesh/atlas_MeshBuilder_module.F90 b/src/atlas_f/mesh/atlas_MeshBuilder_module.F90 new file mode 100644 index 000000000..5b1fce4dd --- /dev/null +++ b/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 diff --git a/src/tests/mesh/CMakeLists.txt b/src/tests/mesh/CMakeLists.txt index c901de5bc..f9d7825e3 100644 --- a/src/tests/mesh/CMakeLists.txt +++ b/src/tests/mesh/CMakeLists.txt @@ -24,6 +24,14 @@ if( HAVE_FCTEST ) ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} ) + add_fctest( TARGET atlas_fctest_mesh_triangular_mesh_builder + LINKER_LANGUAGE Fortran + CONDITION ON + SOURCES fctest_mesh_triangular_mesh_builder.F90 + LIBS atlas_f + ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} + ) + add_fctest( TARGET atlas_fctest_connectivity LINKER_LANGUAGE Fortran SOURCES fctest_connectivity.F90 @@ -135,6 +143,12 @@ ecbuild_add_test( TARGET atlas_test_mesh_builder_parallel ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} ) +ecbuild_add_test( TARGET atlas_test_mesh_triangular_mesh_builder + SOURCES test_mesh_triangular_mesh_builder.cc + LIBS atlas + ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} +) + ecbuild_add_executable( TARGET atlas_test_mesh_reorder SOURCES test_mesh_reorder.cc LIBS atlas diff --git a/src/tests/mesh/fctest_mesh_triangular_mesh_builder.F90 b/src/tests/mesh/fctest_mesh_triangular_mesh_builder.F90 new file mode 100644 index 000000000..63e70e2bf --- /dev/null +++ b/src/tests/mesh/fctest_mesh_triangular_mesh_builder.F90 @@ -0,0 +1,132 @@ +! (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. + +! This File contains Unit Tests for testing the +! C++ / Fortran Interfaces to the Mesh Datastructure +! @author Willem Deconinck + +#include "fckit/fctest.h" + +! ----------------------------------------------------------------------------- + +TESTSUITE(fctest_atlas_TriangularMeshBuilder) + +! ----------------------------------------------------------------------------- + +TESTSUITE_INIT + use atlas_module + call atlas_library%initialise() +END_TESTSUITE_INIT + +! ----------------------------------------------------------------------------- + +TESTSUITE_FINALIZE + use atlas_module + call atlas_library%finalise() +END_TESTSUITE_FINALIZE + +! ----------------------------------------------------------------------------- + +TEST( test_mesbuilder ) + use, intrinsic :: iso_c_binding + use atlas_module + implicit none + type(atlas_Mesh) :: mesh + + integer :: nb_nodes + integer :: nb_triags + integer(ATLAS_KIND_GIDX), allocatable :: node_global_index(:) + integer(ATLAS_KIND_GIDX), allocatable :: triag_global_index(:) + integer(ATLAS_KIND_GIDX), allocatable :: triag_nodes_global_index(:,:) + real(c_double), allocatable :: x(:), y(:), lon(:), lat(:) + + ! small mesh + ! + ! 1 ---- 5 ------- 6 + ! | 3 / 4 \ 1 / 2 | + ! 2 ------ 3 ----- 4 + ! + nb_nodes = 6 + nb_triags = 4 + + allocate(node_global_index(nb_nodes)) + allocate(x(nb_nodes)) + allocate(y(nb_nodes)) + allocate(lon(nb_nodes)) + allocate(lat(nb_nodes)) + allocate(triag_global_index(nb_triags)) + allocate(triag_nodes_global_index(3,nb_triags)) + + node_global_index = [1, 2, 3, 4, 5, 6] + lon = [0.0, 0.0, 10.0, 15.0, 5.0, 15.0] + lat = [5.0, 0.0, 0.0, 0.0, 5.0, 5.0] + x = lon / 10._c_double + y = lat / 10._c_double + + triag_global_index = [1, 2, 3, 4] + triag_nodes_global_index = reshape([3,6,5 , 3,4,6, 2,5,1, 2,3,5], shape(triag_nodes_global_index)) + + + ! Manually build mesh + BLOCK + type(atlas_TriangularMeshBuilder) :: meshbuilder + meshbuilder = atlas_TriangularMeshBuilder() + mesh = meshbuilder%build(nb_nodes, node_global_index, x, y, lon, lat, & + nb_triags, triag_global_index, triag_nodes_global_index) + call meshbuilder%final() + END BLOCK + + ! Mesh should created now. Verify some fields + BLOCK + type(atlas_mesh_Nodes) :: nodes + real(c_double), pointer :: xy(:,:) + real(c_double), pointer :: lonlat(:,:) + integer(ATLAS_KIND_GIDX), pointer :: global_index(:) + + type(atlas_Field) :: field_xy + type(atlas_Field) :: field_lonlat + type(atlas_Field) :: field_global_index + + integer :: jnode + integer :: nb_nodes + nodes = mesh%nodes() + nb_nodes = nodes%size() + FCTEST_CHECK_EQUAL( nb_nodes, 6 ) + field_xy = nodes%xy() + field_lonlat = nodes%lonlat() + field_global_index = nodes%global_index() + call field_xy%data(xy) + call field_lonlat%data(lonlat) + call field_global_index%data(global_index) + do jnode=1,nb_nodes + FCTEST_CHECK_EQUAL( xy(:,jnode), ([x(jnode), y(jnode)])) + FCTEST_CHECK_EQUAL( lonlat(:,jnode), ([lon(jnode), lat(jnode)])) + FCTEST_CHECK_EQUAL( global_index(jnode), jnode) + enddo + call field_xy%final() + call field_lonlat%final() + call field_global_index%final() + call nodes%final() + END BLOCK + + ! Output mesh + BLOCK + type(atlas_Output) :: gmsh + gmsh = atlas_output_Gmsh("out.msh",coordinates="lonlat") + call gmsh%write(mesh) + call gmsh%final() + END BLOCK + + + call mesh%final() + +END_TEST + +! ----------------------------------------------------------------------------- + +END_TESTSUITE + diff --git a/src/tests/mesh/test_mesh_triangular_mesh_builder.cc b/src/tests/mesh/test_mesh_triangular_mesh_builder.cc new file mode 100644 index 000000000..27ed21a50 --- /dev/null +++ b/src/tests/mesh/test_mesh_triangular_mesh_builder.cc @@ -0,0 +1,73 @@ +/* + * (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 +#include +#include + +#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 "atlas/output/Gmsh.h" + +#include "tests/AtlasTestEnvironment.h" + +using namespace atlas::mesh; + +//#include "atlas/output/Gmsh.h" +//using namespace atlas::output; + +namespace atlas { +namespace test { + +//----------------------------------------------------------------------------- + +CASE("test_tiny_mesh") { + // small regional grid whose cell-centers are connected as (global nodes and cells): + // + // 1 ---- 5 ----- 6 + // | 3 / 4 | 1 /2 | + // 2 ----- 3 ---- 4 + // + size_t nb_nodes = 6; + std::vector lon{{0.0, 0.0, 10.0, 15.0, 5.0, 15.0}}; + std::vector lat{{5.0, 0.0, 0.0, 0.0, 5.0, 5.0}}; + std::vector x(nb_nodes); + std::vector y(nb_nodes); + for (size_t j=0; j global_index(6); + std::iota(global_index.begin(), global_index.end(), 1); // 1-based numbering + + // triangles + size_t nb_triags = 4; + std::vector> triag_nodes_global = {{{3, 6, 5}}, {{3, 4, 6}}, {{2, 5, 1}}, {{2, 3, 5}}}; + std::vector triag_global_index = {1, 2, 3, 4}; + + const TriangularMeshBuilder mesh_builder{}; + const Mesh mesh = mesh_builder(nb_nodes, global_index.data(), x.data(), y.data(), lon.data(), lat.data(), + nb_triags, triag_global_index.data(), triag_nodes_global.data()->data()); + + output::Gmsh gmsh("out.msh", util::Config("coordinates", "xy")); + gmsh.write(mesh); +} + +//----------------------------------------------------------------------------- + +} // namespace test +} // namespace atlas + +int main(int argc, char** argv) { + return atlas::test::run(argc, argv); +}