Merge pull request #1180 from LLNL/feature/whitlock/hilbert_ordering

Added utility functions helpful in partitioning and reordering topologies.

CHANGELOG.md

@@ -20,10 +20,16 @@ and this project aspires to adhere to [Semantic Versioning](https://semver.org/s
 - Added a `conduit::blueprint::mpi::mesh::utils::adjset::compare_pointwise()` function that can compare adjsets for multi-domain meshes in parallel. The function is used to diagnose adjsets with points that are out of order on either side of the boundary. The comparison is done point by point within each group and it checks to ensure that the points reference the same spatial location.
 - Added a `conduit::blueprint::mesh::utils::reorder()` function that can accept a vector of element ids and create a reordered topology. The points and coordinates are re-ordered according to their first use in the new element ordering.
 - Added a `conduit::blueprint::mesh::utils::topology::spatial_ordering()` function that takes a topology and computes centroids for each element, passes them through a kdtree, and returns the new element ordering. The new ordering can be used with the `reorder()` function.
+- Added a `conduit::blueprint::mesh::utils::topology::hilbert_ordering()` function that computes a new order for a topology's elements based on their centroids and a Hilbert curve. The new ordering can be used with the `reorder()` function.
 - Added a `conduit::blueprint::mesh::utils::slice_array()` function that can slice Conduit nodes that contain arrays. A new node with the same type is created but it contains only the selected indices.
 - Added a `conduit::blueprint::mesh::utils::slice_field()` function. It is like `slice_array()` but it can handle the mcarray protocol. This functionality was generalized from the partitioner.
+- Added a `conduit::blueprint::mesh::utils::topology::unstructured::rewrite_connectivity()` function that will rewrite a topology's connectivity in terms of a different coordset. The PointQuery is used internally to search for equivalent coordinates in the new coordset.
+- Added a `conduit::blueprint::mesh::utils::copy_fields()` function that helps copy fields from one fields node to another.
+- Added a `conduit::blueprint::mesh::utils::convert()` function that converts a list of nodes to a desired data type.
+- Added a `conduit::blueprint::mesh::generate_boundary_partition_field()` function that can take a topology and a partition field and generate a field for a related boundary topology. This is helpful when partitioning a boundary topology in the same manner as its parent topology.
 - Added `blueprint.mesh.examples.strided_structured` to the blueprint python module.
 
+
 ### Changed
 
 #### General

src/libs/blueprint/conduit_blueprint_mesh.cpp

@@ -7507,6 +7507,7 @@ mesh::partition(const conduit::Node &n_mesh,
     }
 }
 
+//-------------------------------------------------------------------------
 void mesh::partition_map_back(const Node& repart_mesh,
                               const Node& options,
                               Node& orig_mesh)
@@ -7515,6 +7516,78 @@ void mesh::partition_map_back(const Node& repart_mesh,
     p.map_back_fields(repart_mesh, options, orig_mesh);
 }
 
+//-------------------------------------------------------------------------
+void mesh::generate_boundary_partition_field(const conduit::Node &topo,
+                                             const conduit::Node &partField,
+                                             const conduit::Node &btopo,
+                                             conduit::Node &bpartField)
+{
+    // Basic checks.
+    if(topo.has_child("coordset") && btopo.has_child("coordset"))
+    {
+        if(topo.fetch_existing("coordset").as_string() != btopo.fetch_existing("coordset").as_string())
+        {
+            CONDUIT_ERROR("Input topologies must use the same coordset.");
+        }
+    }
+    if(partField.has_child("topology") && (partField.fetch_existing("topology").as_string() != topo.name()))
+    {
+        CONDUIT_ERROR("The partition field must be associated with the " << topo.name() << " topology.");
+    }
+
+    const Node *coordset = bputils::find_reference_node(topo, "coordset");
+
+    // Produce the external "faces" of the domain and for each "face", hash
+    // its node ids and associate that hash with the parent zone for the face.
+    auto d = static_cast<size_t>(conduit::blueprint::mesh::utils::topology::dims(btopo));
+    std::vector<std::pair<size_t, size_t>> desired_maps{{d, d + 1}};
+    bputils::TopologyMetadata md(topo, *coordset, d, desired_maps);
+    const conduit::Node &dtopo = md.get_topology(d);
+    auto nent = md.get_topology_length(d);
+    std::map<conduit::uint64, int> hashToZone;
+    for (conduit::index_t ei = 0; ei < nent; ei++)
+    {
+        // The global associatino for this "face" tells us how many elements
+        // it belongs to. If there is just one parent, it is external and can
+        // belong to the boundary.
+        const auto vv = md.get_global_association(ei, d, d + 1);
+        if (vv.size() == 1)
+        {
+            // Get the ids that make up the entity and hash them.
+            auto ids = conduit::blueprint::mesh::utils::topology::unstructured::points(dtopo, ei);
+            std::sort(ids.begin(), ids.end());
+            conduit::uint64 h = conduit::utils::hash(&ids[0], static_cast<unsigned int>(ids.size()));
+
+            // Save hash to parent zone.
+            hashToZone[h] = vv[0];
+        }
+    }
+
+    // Get the partition field.
+    const auto f = partField.fetch_existing("values").as_int32_accessor();
+
+    // Now, iterate through the boundary topology, hash each entity's ids
+    // and try to look up the parent zone. The hashToZone map should contain
+    // all possible external faces for the domain so the boundary should be
+    // a subset of that.
+    auto blen = conduit::blueprint::mesh::topology::length(btopo);
+    bpartField.reset();
+    bpartField["association"] = "element";
+    bpartField["topology"] = btopo.name();
+    bpartField["values"].set(conduit::DataType::int32(blen));
+    auto bndPartition = bpartField["values"].as_int32_ptr();
+    for (conduit::index_t ei = 0; ei < blen; ei++)
+    {
+        // Get the ids that make up the entity and hash them.
+        auto ids = conduit::blueprint::mesh::utils::topology::unstructured::points(btopo, ei);
+        std::sort(ids.begin(), ids.end());
+        conduit::uint64 h = conduit::utils::hash(&ids[0], static_cast<unsigned int>(ids.size()));
+
+        auto it = hashToZone.find(h);
+        bndPartition[ei] = (it != hashToZone.end()) ? f[it->second] : 0;
+    }
+}
+
 //-------------------------------------------------------------------------
 void
 mesh::flatten(const conduit::Node &mesh,

src/libs/blueprint/conduit_blueprint_mesh.hpp

@@ -162,7 +162,7 @@ void CONDUIT_BLUEPRINT_API paint_adjset(const std::string &adjset_name,
 /**
  @brief Partition an input mesh or set of mesh domains into a different decomposition,
         according to options. This is the serial implementation.
- @param n_mesh  A Conduit node containing a Blueprint mesh or set of mesh domains.
+ @param mesh    A Conduit node containing a Blueprint mesh or set of mesh domains.
  @param options A Conduit node containing options that govern the partitioning.
  @param[out]    A Conduit node to accept the repartitioned mesh(es).
  */
@@ -188,7 +188,24 @@ void CONDUIT_BLUEPRINT_API partition(const conduit::Node &mesh,
 void CONDUIT_BLUEPRINT_API partition_map_back(const conduit::Node& repart_mesh,
                                               const conduit::Node& options,
                                               conduit::Node& orig_mesh);
+//-------------------------------------------------------------------------
+/**
+ @brief Take a partition field on the source mesh and make a corresponding field
+        on the boundary mesh so it can be partitioned in a compatible way in
+        another call to partition.
+
+ @param topo       The source topology.
+ @param partField  The partition field.
+ @param btopo      The boundary topology.
+ @param bpartField The node that will contain the new field.
 
+ */
+void CONDUIT_BLUEPRINT_API generate_boundary_partition_field(const conduit::Node &topo,
+                                                             const conduit::Node &partField,
+                                                             const conduit::Node &btopo,
+                                                             conduit::Node &bpartField);
+
+//-------------------------------------------------------------------------
 /**
  @brief Convert the given blueprint mesh into a blueprint table.
  @param mesh    A Conduit node containing a blueprint mesh or set of mesh domains.

src/libs/blueprint/conduit_blueprint_mesh_examples_tiled.cpp

@@ -632,9 +632,9 @@ Tiler::generate(conduit::index_t nx, conduit::index_t ny, conduit::index_t nz,
     if(options.has_path("tile"))
         initialize(options.fetch_existing("tile"));
 
-    bool reorder = true;
+    std::string reorder;
     if(options.has_path("reorder"))
-        reorder = options.fetch_existing("reorder").as_string() == "kdtree";
+        reorder = options.fetch_existing("reorder").as_string();
 
     if(options.has_path("meshname"))
         meshName = options.fetch_existing("meshname").as_string();
@@ -882,13 +882,18 @@ Tiler::generate(conduit::index_t nx, conduit::index_t ny, conduit::index_t nz,
 
     // Reorder the elements unless it was turned off.
     std::vector<conduit::index_t> old2NewPoint;
-    if(reorder)
+    bool doReorder = (reorder == "kdtree" || reorder == "hilbert");
+    if(doReorder)
     {
         // We need offsets.
         conduit::blueprint::mesh::utils::topology::unstructured::generate_offsets(topo, topo["elements/offsets"]);
 
         // Create a new order for the mesh elements.
-        const auto elemOrder = conduit::blueprint::mesh::utils::topology::spatial_ordering(topo);
+        std::vector<conduit::index_t> elemOrder;
+        if(reorder == "kdtree")
+            elemOrder = conduit::blueprint::mesh::utils::topology::spatial_ordering(topo);
+        else
+            elemOrder = conduit::blueprint::mesh::utils::topology::hilbert_ordering(topo);
 
 #ifdef CONDUIT_USE_PARTITIONER_FOR_REORDER
         // NOTE: This was an idea I had after I made reorder. Reordering is like
@@ -934,7 +939,7 @@ Tiler::generate(conduit::index_t nx, conduit::index_t ny, conduit::index_t nz,
     {
         // 2D
         bshape = "line";
-        if(reorder)
+        if(doReorder)
         {
             iterateBoundary2D(tiles, nx, ny, flags,
                 [&](const conduit::index_t *ids, conduit::index_t npts, int bnd)
@@ -962,7 +967,7 @@ Tiler::generate(conduit::index_t nx, conduit::index_t ny, conduit::index_t nz,
         // 3D
         bshape = "quad";
         bool anyNonQuads = false;
-        if(reorder)
+        if(doReorder)
         {
             iterateBoundary3D(tiles, nx, ny, nz, ptsPerPlane, flags,
                 [&](const conduit::index_t *ids, conduit::index_t npts, int bnd)
@@ -1011,7 +1016,7 @@ Tiler::generate(conduit::index_t nx, conduit::index_t ny, conduit::index_t nz,
     }
 
     // Build an adjacency set.
-    addAdjset(tiles, nx, ny, nz, ptsPerPlane, reorder, old2NewPoint, options, res);
+    addAdjset(tiles, nx, ny, nz, ptsPerPlane, doReorder, old2NewPoint, options, res);
 
 #if 0
     // Print for debugging.

src/libs/blueprint/conduit_blueprint_mesh_topology_metadata.cpp

@@ -3647,7 +3647,6 @@ TopologyMetadata::TopologyMetadata(const conduit::Node &topology,
     const Node &topo_elem_offsets = dim_topos[topo_shape.dim]["elements/offsets"];
     const auto elem_conn_access = topo_elem_conn.as_index_t_accessor();
     const auto elem_offsets_access = topo_elem_offsets.as_index_t_accessor();
-    index_t topo_conn_len = 0;
     for(index_t ei = 0; ei < topo_num_elems; ei++)
     {
         index_t bi = topo_num_coords + ei;
@@ -3656,7 +3655,6 @@ TopologyMetadata::TopologyMetadata(const conduit::Node &topology,
         index_t entity_end_index = (ei < topo_num_elems - 1) ? elem_offsets_access[ei + 1] : 
                                        topo_elem_conn.dtype().number_of_elements();
         index_t entity_size = entity_end_index - entity_start_index;
-        topo_conn_len += entity_size;
 
         // Get the vector we'll populate.
         std::vector<int64> &elem_indices = entity_index_bag[bi];