diff --git a/CHANGELOG.md b/CHANGELOG.md index 7575c9a1b..39f6a88a2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,10 +10,18 @@ and this project aspires to adhere to [Semantic Versioning](https://semver.org/s #### General - Added `conduit_json_external` protocol. Creates a json schema representation of a node that includes all addresses that the node is pointing to. Parsing this schema will create a node equivalent to `set_external()`. +- Added a `conduit_generate_data` executable that can generate datasets using the `tiled()` and `braid()` functions and save the datasets to files. #### Relay - Added ability to read N-dimensional hyperslabs from HDF5 leaf arrays into linear memory arrays. +#### Blueprint +- Added a `conduit::blueprint::mesh::examples::tiled()` function that can generate meshes by repeating a tiled pattern. +- 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::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. ### Changed diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index a24e64517..0f9ea7f70 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -15,7 +15,7 @@ endif() # Conduit ################################ -project(conduit VERSION "0.8.8") +project(conduit VERSION "0.8.9") ################################ # Build Options diff --git a/src/docs/sphinx/blueprint_mesh.rst b/src/docs/sphinx/blueprint_mesh.rst index 9324fe0b5..0f14bd6b8 100644 --- a/src/docs/sphinx/blueprint_mesh.rst +++ b/src/docs/sphinx/blueprint_mesh.rst @@ -1687,6 +1687,91 @@ equal to half the number of prisms. The resulting data is placed the Node ``res``, which is passed in via reference. + +tiled ++++++++++ + +The ``tiled()`` function repeats a tile (given as a 2D Blueprint topology) into a larger mesh composed of +a regular square tiling of the supplied tile. If no topology is given, a default pattern consisting of quads +is used. The output mesh can be either 2D or 3D. For 3D, supply a ``nz`` parameter greater than zero. Note +that the input tile must consist of a homogeneous set of triangles or quads to extrude the tile into 3D since +polyhedral output is not yet supported. The ``tiled()`` function produces a single domain comprised of a +main mesh, a boundary mesh, and adjacency sets if the output mesh is to be part of a multi-domain dataset. + +.. code:: cpp + + conduit::blueprint::mesh::examples::tiled(index_t nx, // number of tiles along x + index_t ny, // number of tiles along y + index_t nz, // number of elements along z (0 for 2D) + conduit::Node &res, // result container + const conduit::Node &options);// options node + +The ``tiled()`` function accepts a Conduit node containing options that influence how the mesh is generated. +If the options contain a ``tile`` node that contains a 2D blueprint topology, the first supplied topology will +be used to override the default tile pattern. The ``reorder`` option indicates the type of point and element +reordering that will be done. Reordering can improve cache-friendliness. The default is to reorder points +and elements, using "kdtree" method. Passing "none" or an empty string will prevent reordering. +The name of the mesh can be given by passing a ``meshname`` option string. Likewise, the name of the boundary +mesh can be supplied using the ``boundarymeshname`` option. The optional ``translate/x`` and ``translate/y`` +options determine the tile spacing. If the translation values are not given, they will be determined from +the coordset extents. The output mesh topology will store its integer connectivity information as index_t +by default. The precision of the integer output can turned to int32 by passing a ``datatype`` option containing +the "int", "int32", "integer" strings. + +An important set of options define the left, right, bottom, and top sets of points within the supplied tile +pattern. The values in the ``left`` option identify the list of points that define the left edge of the tile. +These are indices into the coordset and the values should be in consecutive order along the edge. Opposite point +sets must match. In other words, the left and right point sets must contain the same number of points and +they need to proceed along their edges in the same order. The same is true of the bottom and top point sets. + +The ``tiled()`` function also accepts options that simplify the task of generating +mesh domains for a multi-domain dataset. The coordinate extents of the current mesh domain are given using +the ``extents`` option, which contains 6 double values: {xmin, xmax, ymin, ymax, zmin, zmax}. The ``domains`` +option contains a triple of {domainsI, domainsJ, domainsK} values that indicate how many divisions there are +of the extents in the I,J,K dimensions. The product of these numbers determines the total number of domains. +The ``domain`` option specifies a triple indicating the I,J,K domain id within the overall set of domains. +This is used to help construct adjacency sets. + +.. code:: yaml + + # Define the tile topology + tile: + coordsets: + coords: + type: explicit + values: + x: [0., 1., 2., 0., 1., 2., 0., 1., 2.] + y: [0., 0.5, 0., 1., 1.5, 1., 2., 2.5, 2.] + topologies: + tile: + type: unstructured + coordset: coords + elements: + shape: tri + connectivity: [0,1,4, 0,4,3, 1,2,5, 1,5,4, 3,4,7, 3,7,6, 4,5,8, 4,8,7] + # Set some options that aid tiling. + reorder: 1 + left: [0,3,6] + right: [2,5,8] + bottom: [0,1,2] + top: [6,7,8] + translate: + x: 2. + y: 2. + +.. figure:: tiled_single.png + :width: 400px + :align: center + + Pseudocolor plot of zoneid for default tile mesh that has been reordered. + +.. figure:: tiled.png + :width: 600px + :align: center + + Subset plots of multi-domain datasets created using the ``tiled()`` function. + + miscellaneous ++++++++++++++ diff --git a/src/docs/sphinx/index.rst b/src/docs/sphinx/index.rst index 49e104037..a931196fb 100644 --- a/src/docs/sphinx/index.rst +++ b/src/docs/sphinx/index.rst @@ -115,7 +115,7 @@ Contributors - Mark Miller (LLNL) - Todd Gamblin (LLNL) - Kevin Huynh (LLNL) -- Brad Whitlock (Intelligent Light) +- Brad Whitlock (LLNL) - Chris Laganella (Intelligent Light) - George Aspesi (Harvey Mudd) - Justin Bai (Harvey Mudd) diff --git a/src/docs/sphinx/tiled.png b/src/docs/sphinx/tiled.png new file mode 100644 index 000000000..9229398b4 Binary files /dev/null and b/src/docs/sphinx/tiled.png differ diff --git a/src/docs/sphinx/tiled_single.png b/src/docs/sphinx/tiled_single.png new file mode 100644 index 000000000..e008b05de Binary files /dev/null and b/src/docs/sphinx/tiled_single.png differ diff --git a/src/executables/CMakeLists.txt b/src/executables/CMakeLists.txt index bd9ac22ed..ba8b2d069 100644 --- a/src/executables/CMakeLists.txt +++ b/src/executables/CMakeLists.txt @@ -3,3 +3,4 @@ # other details. No copyright assignment is required to contribute to Conduit. add_subdirectory(adjset_validate) +add_subdirectory(generate_data) diff --git a/src/executables/generate_data/CMakeLists.txt b/src/executables/generate_data/CMakeLists.txt new file mode 100644 index 000000000..cfa521edb --- /dev/null +++ b/src/executables/generate_data/CMakeLists.txt @@ -0,0 +1,36 @@ +# Copyright (c) Lawrence Livermore National Security, LLC and other Conduit +# Project developers. See top-level LICENSE AND COPYRIGHT files for dates and +# other details. No copyright assignment is required to contribute to Conduit. + +if(ENABLE_UTILS) + blt_add_executable( + NAME conduit_generate_data + SOURCES conduit_generate_data.cpp + OUTPUT_DIR ${CMAKE_CURRENT_BINARY_DIR} + DEPENDS_ON conduit conduit_blueprint conduit_relay + FOLDER utils + ) + + # add install target + install(TARGETS conduit_generate_data + RUNTIME DESTINATION bin) + + ################################################################ + # If we have mpi, add a parallel version. + ################################################################ + + if(MPI_FOUND) + blt_add_executable( + NAME conduit_generate_data_mpi + SOURCES conduit_generate_data.cpp + OUTPUT_DIR ${CMAKE_CURRENT_BINARY_DIR} + DEFINES CONDUIT_PARALLEL + DEPENDS_ON conduit conduit_blueprint conduit_relay conduit_relay_mpi_io ${conduit_blt_mpi_deps} + FOLDER utils + ) + + # add install target + install(TARGETS conduit_generate_data_mpi + RUNTIME DESTINATION bin) + endif() +endif() diff --git a/src/executables/generate_data/conduit_generate_data.cpp b/src/executables/generate_data/conduit_generate_data.cpp new file mode 100644 index 000000000..284b77b54 --- /dev/null +++ b/src/executables/generate_data/conduit_generate_data.cpp @@ -0,0 +1,290 @@ +// Copyright (c) Lawrence Livermore National Security, LLC and other Conduit +// Project developers. See top-level LICENSE AND COPYRIGHT files for dates and +// other details. No copyright assignment is required to contribute to Conduit. + +//----------------------------------------------------------------------------- +/// +/// file: conduit_generate_data.hpp +/// +//----------------------------------------------------------------------------- +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#ifdef CONDUIT_PARALLEL +#include +#include +#else +#include +#endif + +//#define MAKE_UNIQUE(T) std::make_unique() +#define MAKE_UNIQUE(T) std::unique_ptr(new T()) + +//----------------------------------------------------------------------------- +/// Generate data for a domain +class DomainGenerator +{ +public: + virtual void generate(int domain[3], conduit::Node &n, conduit::Node &opts) = 0; + + void setDims(const int d[3]) + { + for(int i = 0; i < 3; i++) + dims[i] = d[i]; + } + void setDomains(const int d[3]) + { + for(int i = 0; i < 3; i++) + domains[i] = d[i]; + } + void setExtents(const double d[3]) + { + for(int i = 0; i < 6; i++) + extents[i] = d[i]; + } + void setMeshType(const std::string &m) + { + meshType = m; + } +protected: + int dims[3]{0,0,0}; + int domains[3]{1,1,1}; + double extents[6]{0., 1., 0., 1., 0., 1.}; + std::string meshType{}; +}; + +//----------------------------------------------------------------------------- +/// Generate data using tiled data generator. +class TiledDomainGenerator : public DomainGenerator +{ +public: + virtual void generate(int domain[3], conduit::Node &n, conduit::Node &opts) override + { + // Determine the size and location of this domain in the whole. + double sideX = (extents[1] - extents[0]) / static_cast(domains[0]); + double sideY = (extents[3] - extents[2]) / static_cast(domains[1]); + double sideZ = (extents[5] - extents[4]) / static_cast(domains[2]); + double domainExt[] = {extents[0] + domain[0] * sideX, + extents[0] + (domain[0]+1) * sideX, + extents[2] + domain[1] * sideY, + extents[2] + (domain[1]+1) * sideY, + extents[4] + domain[2] * sideZ, + extents[4] + (domain[2]+1) * sideZ}; + opts["extents"].set(domainExt, 6); + opts["domain"].set(domain, 3); + opts["domains"].set(domains, 3); + + conduit::blueprint::mesh::examples::tiled(dims[0], dims[1], dims[2], n, opts); + } +}; + +//----------------------------------------------------------------------------- +/// Generate data using braid data generator. +class BraidDomainGenerator : public DomainGenerator +{ +public: + virtual void generate(int /*domain*/[3], conduit::Node &n, conduit::Node &) override + { + conduit::blueprint::mesh::examples::braid(meshType, dims[0] + 1, dims[1] + 1, dims[2] + 1, n); + + // TODO: Use domain,domains to adjust coordinates to get them to line up nicely. + } +}; + +//----------------------------------------------------------------------------- +void +printUsage(const char *exeName) +{ + std::cout << "Usage: " << exeName << "[-dims x,y,z] [-domains x,y,z] [-tile]\n" + << " [-braid] [-output fileroot] [-protocol name] [-meshtype type]\n" + << " [-tiledef filename] [-extents x0,x1,y0,y1[,z0,z1]] [-help]\n"; + std::cout << "\n"; + std::cout << "Argument Description\n"; + std::cout << "=================== ==========================================================\n"; + std::cout << "-dims x,y,z The number of mesh zones in each dimension. For 2D meshes,\n"; + std::cout << " supply 0 for z.\n"; + std::cout << "\n"; + std::cout << "-domains x,y,z The number of domains in each dimension.\n"; + std::cout << "\n"; + std::cout << "-tile Generate a mesh using the tiled data generator.\n"; + std::cout << "\n"; + std::cout << "-braid Generate a mesh using the braid data generator.\n"; + std::cout << "\n"; + std::cout << "-output fileroot Specify the root used in filenames that are created.\n"; + std::cout << "\n"; + std::cout << "-protocol name Specify the protocol used in writing the data. The default\n"; + std::cout << " is \"hdf5\".\n"; + std::cout << "\n"; + std::cout << "-meshtype type The mesh type used when generating data using braid.\n"; + std::cout << "\n"; + std::cout << "-tiledef filename A file containing a tile definition.\n"; + std::cout << "\n"; + std::cout << "-extents ext A list of 4 or 6 comma-separated values indicating extents\n"; + std::cout << " as pairs of min,max values for each dimension.\n"; + std::cout << "\n"; + std::cout << "-help Print the usage and exit.\n"; +} + +//----------------------------------------------------------------------------- +int +main(int argc, char *argv[]) +{ + int rank = 0; +#ifdef CONDUIT_PARALLEL + int size = 1; + MPI_Init(&argc, &argv); + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + MPI_Comm_size(MPI_COMM_WORLD, &size); +#endif + + // Some basic arg parsing. + int dims[3] = {10, 10, 10}; + int domains[3] = {1, 1, 1}; + double extents[6] = {0., 1., 0., 1., 0., 1.}; + conduit::Node n, opts; + std::unique_ptr g = MAKE_UNIQUE(TiledDomainGenerator); + std::string meshType("hexs"),meshTypeDefault("hexs"), output("output"), protocol("hdf5"); + for(int i = 1; i < argc; i++) + { + if(strcmp(argv[i], "-dims") == 0 && (i+1) < argc) + { + int d[3] = {1,1,1}; + if(sscanf(argv[i+1], "%d,%d,%d", &d[0], &d[1], &d[2]) == 3) + { + dims[0] = std::max(d[0], 1); + dims[1] = std::max(d[1], 1); + dims[2] = std::max(d[2], 0); // Allow 0 for 2D + + // If we have not set the mesh type, set it according to dimension. + if(meshType == meshTypeDefault) + { + meshType = (dims[2] > 0) ? "hexs" : "quads"; + } + } + i++; + } + else if(strcmp(argv[i], "-domains") == 0 && (i+1) < argc) + { + int d[3] = {1,1,1}; + if(sscanf(argv[i+1], "%d,%d,%d", &d[0], &d[1], &d[2]) == 3) + { + domains[0] = std::max(d[0], 1); + domains[1] = std::max(d[1], 1); + domains[2] = std::max(d[2], 1); + } + i++; + } + else if(strcmp(argv[i], "-tile") == 0) + { + g = MAKE_UNIQUE(TiledDomainGenerator); + } + else if(strcmp(argv[i], "-braid") == 0) + { + g = MAKE_UNIQUE(BraidDomainGenerator); + } + else if(strcmp(argv[i], "-output") == 0 && (i+1) < argc) + { + output = argv[i+1]; + i++; + } + else if(strcmp(argv[i], "-protocol") == 0 && (i+1) < argc) + { + protocol = argv[i+1]; + i++; + } + else if(strcmp(argv[i], "-meshtype") == 0 && (i+1) < argc) + { + meshType = argv[i+1]; + i++; + } + else if(strcmp(argv[i], "-tiledef") == 0 && (i+1) < argc) + { + // Load a tile definition file into the options tile node. + conduit::relay::io::load(argv[i+1], opts["tile"]); + i++; + } + else if(strcmp(argv[i], "-extents") == 0 && (i+1) < argc) + { + double e[6] = {0., 1., 0., 1., 0., 1.}; + if(sscanf(argv[i + 1], "%lg,%lg,%lg,%lg,%lg,%lg", &e[0], &e[1], &e[2], &e[3], &e[4], &e[5]) == 6) + { + memcpy(extents, e, 6 * sizeof(double)); + } + else if(sscanf(argv[i + 1], "%lg,%lg,%lg,%lg", &e[0], &e[1], &e[2], &e[3]) == 4) + { + memcpy(extents, e, 4 * sizeof(double)); + } + i++; + } + else if(strcmp(argv[i], "-help") == 0 || + strcmp(argv[i], "--help") == 0) + { + if(rank == 0) + printUsage(argv[0]); +#ifdef CONDUIT_PARALLEL + MPI_Finalize(); +#endif + return 0; + } + } + + g->setDims(dims); + g->setDomains(domains); + g->setMeshType(meshType); + g->setExtents(extents); + + int ndoms = domains[0] * domains[1] * domains[2]; + if(ndoms == 1 && rank == 0) + { + // Single domain. + int domain[] = {0, 0, 0}; + g->generate(domain, n, opts); + } + else + { + int domainid = 0; + for(int k = 0; k < domains[2]; k++) + { + for(int j = 0; j < domains[1]; j++) + { + for(int i = 0; i < domains[0]; i++, domainid++) + { + int domain[] = {i, j, k}; + +#ifdef CONDUIT_PARALLEL + if(domainid % size == rank) + { +#endif + // Make the new domain. + char domainName[32]; + sprintf(domainName, "domain_%07d", domainid); + conduit::Node &d = n[domainName]; + g->generate(domain, d, opts); +#ifdef CONDUIT_PARALLEL + } +#endif + } + } + } + } + + // Save the output domains. +#ifdef CONDUIT_PARALLEL + conduit::relay::mpi::io::blueprint::save_mesh(n, output, protocol, MPI_COMM_WORLD); + + MPI_Finalize(); +#else + conduit::relay::io::save(n, output + "-inspect.yaml", "yaml"); + conduit::relay::io::blueprint::save_mesh(n, output, protocol); +#endif + + return 0; +} diff --git a/src/libs/blueprint/CMakeLists.txt b/src/libs/blueprint/CMakeLists.txt index f7efc3a71..d457699b6 100644 --- a/src/libs/blueprint/CMakeLists.txt +++ b/src/libs/blueprint/CMakeLists.txt @@ -39,6 +39,7 @@ set(blueprint_headers conduit_blueprint_mesh_examples_related_boundary.hpp conduit_blueprint_mesh_examples_polystar.hpp conduit_blueprint_mesh_examples_rz_cylinder.hpp + conduit_blueprint_mesh_examples_tiled.hpp conduit_blueprint_mcarray.hpp conduit_blueprint_mcarray_examples.hpp conduit_blueprint_o2mrelation.hpp @@ -81,6 +82,7 @@ set(blueprint_sources conduit_blueprint_mesh_examples_related_boundary.cpp conduit_blueprint_mesh_examples_polystar.cpp conduit_blueprint_mesh_examples_rz_cylinder.cpp + conduit_blueprint_mesh_examples_tiled.cpp conduit_blueprint_mesh_flatten.cpp conduit_blueprint_mcarray.cpp conduit_blueprint_mcarray_examples.cpp diff --git a/src/libs/blueprint/conduit_blueprint_mesh.cpp b/src/libs/blueprint/conduit_blueprint_mesh.cpp index cbabce8f3..fe32538bf 100644 --- a/src/libs/blueprint/conduit_blueprint_mesh.cpp +++ b/src/libs/blueprint/conduit_blueprint_mesh.cpp @@ -150,7 +150,7 @@ std::vector intersect_sets(const Container1 &v1, } } - return std::vector(std::move(res)); + return res; } //----------------------------------------------------------------------------- @@ -182,7 +182,7 @@ std::vector subtract_sets(const std::vector &v1, res.push_back(v1[i1]); } } - return std::vector(std::move(res)); + return res; } //----------------------------------------------------------------------------- @@ -2008,7 +2008,7 @@ mesh::domains(conduit::Node &n) } } - return std::vector(std::move(doms)); + return doms; } @@ -2035,7 +2035,7 @@ mesh::domains(const conduit::Node &mesh) } } - return std::vector(std::move(doms)); + return doms; } //------------------------------------------------------------------------- @@ -2744,7 +2744,7 @@ group_domains_and_maps(conduit::Node &mesh, conduit::Node &s2dmap, conduit::Node } } - return std::vector(std::move(doms_and_maps)); + return doms_and_maps; } //----------------------------------------------------------------------------- @@ -3596,7 +3596,7 @@ mesh::generate_sides(conduit::Node& mesh, side_dims.push_back(2); } - return std::vector(std::move(side_dims)); + return side_dims; }; verify_generate_mesh(mesh, src_adjset_name); @@ -3644,7 +3644,7 @@ mesh::generate_corners(conduit::Node& mesh, corner_dims.push_back(di); } - return std::vector(std::move(corner_dims)); + return corner_dims; }; verify_generate_mesh(mesh, src_adjset_name); @@ -3966,8 +3966,22 @@ mesh::coordset::generate_strip(const Node& coordset, } } +//----------------------------------------------------------------------------- +void +mesh::coordset::to_explicit(const conduit::Node& coordset, + conduit::Node& coordset_dest) +{ + std::string type = coordset.fetch_existing("type").as_string(); -//------------------------------------------------------------------------- + if(type == "uniform") + mesh::coordset::uniform::to_explicit(coordset, coordset_dest); + else if(type == "rectilinear") + mesh::coordset::rectilinear::to_explicit(coordset, coordset_dest); + else if(type == "explicit") + coordset_dest.set_external(coordset); +} + +//----------------------------------------------------------------------------- void mesh::coordset::uniform::to_rectilinear(const conduit::Node &coordset, conduit::Node &dest) @@ -3976,7 +3990,7 @@ mesh::coordset::uniform::to_rectilinear(const conduit::Node &coordset, } -//------------------------------------------------------------------------- +//----------------------------------------------------------------------------- void mesh::coordset::uniform::to_explicit(const conduit::Node &coordset, conduit::Node &dest) @@ -3985,7 +3999,7 @@ mesh::coordset::uniform::to_explicit(const conduit::Node &coordset, } -//------------------------------------------------------------------------- +//----------------------------------------------------------------------------- void mesh::coordset::rectilinear::to_explicit(const conduit::Node &coordset, conduit::Node &dest) diff --git a/src/libs/blueprint/conduit_blueprint_mesh.hpp b/src/libs/blueprint/conduit_blueprint_mesh.hpp index c5c2ca5cf..c338d22ed 100644 --- a/src/libs/blueprint/conduit_blueprint_mesh.hpp +++ b/src/libs/blueprint/conduit_blueprint_mesh.hpp @@ -355,6 +355,16 @@ namespace coordset void CONDUIT_BLUEPRINT_API generate_strip(const conduit::Node& coordset, conduit::Node& coordset_dest); + //------------------------------------------------------------------------- + /** + @brief Convert the coordset, no matter its type, to explicit. + + @param coordset The input coordset. + @param[out] coordset_dest The output coordset. + */ + void CONDUIT_BLUEPRINT_API to_explicit(const conduit::Node& coordset, + conduit::Node& coordset_dest); + //------------------------------------------------------------------------- // blueprint::mesh::coordset::uniform protocol interface //------------------------------------------------------------------------- diff --git a/src/libs/blueprint/conduit_blueprint_mesh_examples.hpp b/src/libs/blueprint/conduit_blueprint_mesh_examples.hpp index 07a04dd89..6a24d08c1 100644 --- a/src/libs/blueprint/conduit_blueprint_mesh_examples.hpp +++ b/src/libs/blueprint/conduit_blueprint_mesh_examples.hpp @@ -22,6 +22,7 @@ #include "conduit_blueprint_mesh_examples_related_boundary.hpp" #include "conduit_blueprint_mesh_examples_polystar.hpp" #include "conduit_blueprint_mesh_examples_rz_cylinder.hpp" +#include "conduit_blueprint_mesh_examples_tiled.hpp" //----------------------------------------------------------------------------- // -- begin conduit::-- diff --git a/src/libs/blueprint/conduit_blueprint_mesh_examples_tiled.cpp b/src/libs/blueprint/conduit_blueprint_mesh_examples_tiled.cpp new file mode 100644 index 000000000..b54100ede --- /dev/null +++ b/src/libs/blueprint/conduit_blueprint_mesh_examples_tiled.cpp @@ -0,0 +1,1563 @@ +// Copyright (c) Lawrence Livermore National Security, LLC and other Conduit +// Project developers. See top-level LICENSE AND COPYRIGHT files for dates and +// other details. No copyright assignment is required to contribute to Conduit. + +//----------------------------------------------------------------------------- +/// +/// file: conduit_blueprint_mesh_examples_tiled.hpp +/// +//----------------------------------------------------------------------------- + +//----------------------------------------------------------------------------- +// conduit lib includes +//----------------------------------------------------------------------------- +#include "conduit_blueprint_mesh_examples_tiled.hpp" +#include "conduit.hpp" +#include "conduit_blueprint.hpp" +#include "conduit_blueprint_exports.h" +#include "conduit_blueprint_mesh_utils.hpp" + +#include +#include + +// Uncomment this to add some fields on the filed mesh prior to reordering. +// #define CONDUIT_TILER_DEBUG_FIELDS + +// Uncomment this to try an experimental mode that uses the partitioner to +// do reordering. +// #define CONDUIT_USE_PARTITIONER_FOR_REORDER + +// Uncomment this to use a simpler tiled pattern for debugging. +// #define CONDUIT_SIMPLE_TILED_PATTERN + +//----------------------------------------------------------------------------- +// -- begin conduit::-- +//----------------------------------------------------------------------------- +namespace conduit +{ + +//----------------------------------------------------------------------------- +// -- begin conduit::blueprint -- +//----------------------------------------------------------------------------- +namespace blueprint +{ + +//----------------------------------------------------------------------------- +// -- begin conduit::blueprint::mesh -- +//----------------------------------------------------------------------------- +namespace mesh +{ + +//----------------------------------------------------------------------------- +/// Methods that generate example meshes. +//----------------------------------------------------------------------------- +namespace examples +{ + +//----------------------------------------------------------------------------- +/// Detail +//----------------------------------------------------------------------------- +namespace detail +{ + +/** + \brief Keep track of some tile information. + */ +class Tile +{ +public: + static const conduit::index_t INVALID_POINT; + + Tile() : ptids() + { + } + + /// Reset the tile. + void reset(size_t npts) + { + ptids = std::vector(npts, -1); + } + + /// Return the point ids. + std::vector &getPointIds() { return ptids; } + const std::vector &getPointIds() const { return ptids; } + + /// Get the specified point ids for this tile using the supplied indices. + std::vector getPointIds(const std::vector &indices) const + { + std::vector ids; + ids.reserve(indices.size()); + for(const auto &idx : indices) + ids.push_back(ptids[idx]); + return ids; + } + + // Set the point ids + void setPointIds(const std::vector &indices, const std::vector &ids) + { + for(size_t i = 0; i < indices.size(); i++) + { + ptids[indices[i]] = ids[i]; + } + } + +private: + std::vector ptids; //!< This tile's point ids. +}; + +const conduit::index_t Tile::INVALID_POINT = -1; + +/** + \brief Build a mesh from tiles. There is a default tile pattern, although it can + be replaced using an options Node containing new tile information. + */ +class Tiler +{ +public: + static const int BoundaryLeft = 0; + static const int BoundaryRight = 1; + static const int BoundaryBottom = 2; + static const int BoundaryTop = 3; + static const int BoundaryBack = 4; + static const int BoundaryFront = 5; + + static const conduit::index_t InvalidDomain = -1; + + Tiler(); + + /// Generate the tiled mesh. + void generate(conduit::index_t nx, conduit::index_t ny, conduit::index_t nz, + conduit::Node &res, + const conduit::Node &options); +protected: + /// Fill a default tile pattern into the filter. + void initialize(); + + /// Fill in the tile pattern from a Node. + void initialize(const conduit::Node &t); + + /// Return the topology (the first one) + const conduit::Node &getTopology() const + { + const conduit::Node &t = m_options.fetch_existing("topologies"); + return t[0]; + } + + /// Return the coordset + const conduit::Node &getCoordset() const + { + const conduit::Node &t = getTopology(); + std::string coordsetName(t.fetch_existing("coordset").as_string()); + return m_options.fetch_existing("coordsets/" + coordsetName); + } + + /// Return point indices of points along left edge. + const std::vector &left() const { return m_left; } + + /// Return point indices of points along right edge. + const std::vector &right() const { return m_right; } + + /// Return point indices of points along bottom edge. + const std::vector &bottom() const { return m_bottom; } + + /// Return point indices of points along top edge. + const std::vector &top() const { return m_top; } + + /// Return tile width + double width() const { return m_width; } + + /// Return tile height + double height() const { return m_height; } + + /// Creates the points for the tile (if they need to be created). + void addPoints(const double M[3][3], + std::vector &ptids, + std::vector &x, + std::vector &y) + { + // Iterate through points in the template and add them if they have + // not been created yet. + const auto &xpts = getCoordset().fetch_existing("values/x").as_double_array(); + const auto &ypts = getCoordset().fetch_existing("values/y").as_double_array(); + for(conduit::index_t i = 0; i < xpts.number_of_elements(); i++) + { + if(ptids[i] == Tile::INVALID_POINT) + { + ptids[i] = static_cast(x.size()); + + // (x,y,1) * M + double xc = xpts[i] * M[0][0] + ypts[i] * M[1][0] + M[2][0]; + double yc = xpts[i] * M[0][1] + ypts[i] * M[1][1] + M[2][1]; + double h = xpts[i] * M[0][2] + ypts[i] * M[1][2] + M[2][2]; + xc /= h; + yc /= h; + x.push_back(xc); + y.push_back(yc); + } + } + } + + /// Iterate over faces + template + void iterateFaces(const Connectivity &conn, + conduit::index_t nelem, + conduit::index_t sides, + conduit::index_t *idlist, + const std::vector &ptids, + conduit::index_t offset, + bool reverse, + int stype, + Body &&body) const + { + if(reverse) + { + auto s1 = sides - 1; + for(conduit::index_t i = 0; i < nelem; i++) + { + auto start = i * sides; + for(conduit::index_t s = 0; s < sides; s++) + idlist[s1 - s] = offset + ptids[conn[start + s]]; + body(idlist, sides, stype); + } + } + else + { + for(conduit::index_t i = 0; i < nelem; i++) + { + auto start = i * sides; + for(conduit::index_t s = 0; s < sides; s++) + idlist[s] = offset + ptids[conn[start + s]]; + body(idlist, sides, stype); + } + } + } + + /// Iterate over the tile's elements and apply a lambda. + template + void iterateFaces(const std::vector &ptids, + conduit::index_t offset, + bool reverse, + int stype, + Body &&body) const + { + const conduit::Node &topo = getTopology(); + std::string shape = topo.fetch_existing("elements/shape").as_string(); + conduit::index_t sides = 0; + if(shape == "tri") + sides = 3; + else if(shape == "quad") + sides = 4; + + // Handle triangles and quads. + const conduit::Node &n_conn = topo.fetch_existing("elements/connectivity"); + if(sides == 3 || sides == 4) + { + conduit::index_t idlist[4]; + bool handled = false; + if(n_conn.dtype().spanned_bytes() == n_conn.dtype().strided_bytes()) + { + if(n_conn.dtype().is_index_t()) + { + iterateFaces(n_conn.as_index_t_ptr(), + n_conn.dtype().number_of_elements() / sides, + sides, idlist, ptids, + offset, reverse, stype, body); + handled = true; + } + else if(n_conn.dtype().is_int32()) + { + iterateFaces(n_conn.as_int32_ptr(), + n_conn.dtype().number_of_elements() / sides, + sides, idlist, ptids, + offset, reverse, stype, body); + handled = true; + } + } + if(!handled) + { + iterateFaces(n_conn.as_index_t_accessor(), + n_conn.dtype().number_of_elements() / sides, + sides, idlist, ptids, + offset, reverse, stype, body); + } + } + else if(shape == "polygonal") + { + // handle polygons + const auto conn = n_conn.as_index_t_accessor(); + const auto sizes = topo.fetch_existing("elements/sizes").as_index_t_accessor(); + const conduit::index_t nelem = sizes.number_of_elements(); + conduit::index_t start = 0; + std::vector idlist(10); + if(reverse) + { + for(conduit::index_t i = 0; i < nelem; i++) + { + auto esides = sizes[i]; + idlist.reserve(esides); + for(conduit::index_t s = 0; s < esides; s++) + idlist[s] = offset + ptids[conn[start + esides - s]]; + body(&idlist[0], esides, stype); + start += esides; + } + } + else + { + for(conduit::index_t i = 0; i < nelem; i++) + { + auto esides = sizes[i]; + idlist.reserve(esides); + for(conduit::index_t s = 0; s < esides; s++) + idlist[s] = offset + ptids[conn[start + s]]; + body(&idlist[0], esides, stype); + start += esides; + } + } + } + } + + /// Emit the hex cells using this tile's point ids. + void addVolumeElements(const std::vector &ptids, + conduit::index_t plane1Offset, + conduit::index_t plane2Offset, + std::vector &conn, + std::vector &sizes) const + { + const conduit::Node &topo = getTopology(); + std::string shape = topo.fetch_existing("elements/shape").as_string(); + conduit::index_t sides = 0; + if(shape == "tri") + sides = 3; + else if(shape == "quad") + sides = 4; + + if(sides == 3 || sides == 4) + { + const conduit::Node &n_conn = topo.fetch_existing("elements/connectivity"); + const auto tileconn = n_conn.as_index_t_accessor(); + const conduit::index_t nelem = tileconn.number_of_elements() / sides; + for(conduit::index_t i = 0; i < nelem; i++) + { + conduit::index_t start = i * sides; + for(conduit::index_t s = 0; s < sides; s++) + conn.push_back(plane1Offset + ptids[tileconn[start + s]]); + + for(conduit::index_t s = 0; s < sides; s++) + conn.push_back(plane2Offset + ptids[tileconn[start + s]]); + + sizes.push_back(2 * sides); + } + } + else + { + CONDUIT_ERROR("Tiling polygonal shapes into 3D polyhedra is not yet supported."); + } + } + + /// Turn a node into an int vector. + std::vector toIndexVector(const conduit::Node &n) const + { + auto acc = n.as_index_t_accessor(); + std::vector vec; + vec.reserve(acc.number_of_elements()); + for(conduit::index_t i = 0; i < acc.number_of_elements(); i++) + vec.push_back(acc[i]); + return vec; + } + + /// Determine which boundaries are needed. + void boundaryFlags(const conduit::Node &options, bool flags[6]) const; + + // Iterate over 2D boundaries + template + void iterateBoundary2D(const std::vector &tiles, + conduit::index_t nx, + conduit::index_t ny, + const bool flags[6], + Body &&body) const; + + /// Iterate over 3D boundaries. + template + void iterateBoundary3D(const std::vector &tiles, + conduit::index_t nx, + conduit::index_t ny, + conduit::index_t nz, + conduit::index_t nPtsPerPlane, + const bool flags[6], + Body &&body) const; + + /// Compute domain id, if it exists, or return InvalidDomain. + conduit::index_t domainIndex(conduit::index_t d0, conduit::index_t d1, conduit::index_t d2, + const std::vector &domain, + const std::vector &domains) const; + + /// Compute domain id, if it exists, or return InvalidDomain. + conduit::index_t domainIndex(const conduit::index_t delta[3], + const std::vector &domain, + const std::vector &domains) const; + + /// Add ptids to a values node for an adjset, possibly transforming the point ids. + void addAdjsetValues(const std::vector &ptids, + conduit::index_t planeOffset, + bool reorder, + const std::vector &old2NewPoint, + conduit::Node &values) const; + + /// Add adjacency set + void addAdjset(const std::vector &tiles, + conduit::index_t nx, + conduit::index_t ny, + conduit::index_t nz, + conduit::index_t ptsPerPlane, + bool reorder, + const std::vector &old2NewPoint, + const conduit::Node &options, + conduit::Node &out) const; +private: + conduit::Node m_options; + double m_width, m_height; + std::vector m_left, m_right, m_bottom, m_top; + std::string meshName, boundaryMeshName; +}; + +//--------------------------------------------------------------------------- +Tiler::Tiler() : m_options(), m_width(0.), m_height(0.), + m_left(), m_right(), m_bottom(), m_top(), + meshName("mesh"), boundaryMeshName("boundary") +{ + initialize(); +} + +//--------------------------------------------------------------------------- +void +Tiler::initialize() +{ +#ifdef CONDUIT_SIMPLE_TILED_PATTERN + // Simpler case for debugging. + const double x[] = {0., 1., 0., 1.}; + const double y[] = {0., 0., 1., 1.}; + const conduit::index_t conn[] = {0,1,3,2}; + + const conduit::index_t left[] = {0,2}; + const conduit::index_t right[] = {1,3}; + const conduit::index_t bottom[] = {0,1}; + const conduit::index_t top[] = {2,3}; +#else + // Default pattern + const double x[] = { + 0., 3., 10., 17., 20., + 0., 3., 17., 20., + 5., 15., + 7., 10., 13., + 0., 7., 10., 13., 20., + 7., 10., 13., + 5., 15., + 0., 3., 17., 20., + 0., 3., 10., 17., 20., + }; + + const double y[] = { + 0., 0., 0., 0., 0., + 3., 3., 3., 3., + 5., 5., + 7., 7., 7., + 10., 10., 10., 10., 10., + 13., 13., 13., + 15., 15., + 17., 17., 17., 17., + 20., 20., 20., 20., 20., + }; + + const conduit::index_t conn[] = { + // lower-left quadrant + 0,1,6,5, + 1,2,9,6, + 2,12,11,9, + 5,6,9,14, + 9,11,15,14, + 11,12,16,15, + // lower-right quadrant + 2,3,7,10, + 3,4,8,7, + 7,8,18,10, + 2,10,13,12, + 12,13,17,16, + 10,18,17,13, + // upper-left quadrant + 14,22,25,24, + 14,15,19,22, + 15,16,20,19, + 24,25,29,28, + 22,30,29,25, + 19,20,30,22, + // upper-right quadrant + 16,17,21,20, + 17,18,23,21, + 18,27,26,23, + 20,21,23,30, + 23,26,31,30, + 26,27,32,31 + }; + + const conduit::index_t left[] = {0,5,14,24,28}; + const conduit::index_t right[] = {4,8,18,27,32}; + const conduit::index_t bottom[] = {0,1,2,3,4}; + const conduit::index_t top[] = {28,29,30,31,32}; +#endif + + // Define the tile as a topology. + conduit::Node opts; + opts["coordsets/coords/type"] = "explicit"; + opts["coordsets/coords/values/x"].set(x, sizeof(x) / sizeof(double)); + opts["coordsets/coords/values/y"].set(y, sizeof(y) / sizeof(double)); + opts["topologies/tile/type"] = "unstructured"; + opts["topologies/tile/coordset"] = "coords"; + opts["topologies/tile/elements/shape"] = "quad"; + constexpr auto nelem = (sizeof(conn) / sizeof(conduit::index_t)) / 4; + opts["topologies/tile/elements/connectivity"].set(conn, nelem * 4); + std::vector size(nelem, 4); + opts["topologies/tile/elements/sizes"].set(size.data(), size.size()); + + // Define tile boundary indices. + opts["left"].set(left, sizeof(left) / sizeof(conduit::index_t)); + opts["right"].set(right, sizeof(right) / sizeof(conduit::index_t)); + opts["bottom"].set(bottom, sizeof(bottom) / sizeof(conduit::index_t)); + opts["top"].set(top, sizeof(top) / sizeof(conduit::index_t)); + + initialize(opts); +} + +//--------------------------------------------------------------------------- +void +Tiler::initialize(const conduit::Node &t) +{ + std::vector required{"coordsets", "topologies", "left", "right", "bottom", "top"}; + for(const auto &name : required) + { + if(!t.has_child(name)) + { + CONDUIT_ERROR("Node does not contain key: " << name); + } + } + + // Get the tile boundaries and convert them. + m_left = toIndexVector(t.fetch_existing("left")); + m_right = toIndexVector(t.fetch_existing("right")); + m_bottom = toIndexVector(t.fetch_existing("bottom")); + m_top = toIndexVector(t.fetch_existing("top")); + if(m_left.size() != m_right.size()) + { + CONDUIT_ERROR("left/right vectors have different lengths."); + } + if(m_bottom.size() != m_top.size()) + { + CONDUIT_ERROR("bottom/top vectors have different lengths."); + } + + // Save the tile definition options. + m_options.set(t); + + // Make sure the coordset is 2D, explicit. + if(conduit::blueprint::mesh::coordset::dims(getCoordset()) != 2) + { + CONDUIT_ERROR("The tile coordset must be 2D."); + } + if(getCoordset()["type"].as_string() != "explicit") + { + CONDUIT_ERROR("The tile coordset must be explicit."); + } + // Make sure the topology is 2D, unstructured + if(conduit::blueprint::mesh::topology::dims(getTopology()) != 2) + { + CONDUIT_ERROR("The tile topology must be 2D."); + } + if(getTopology()["type"].as_string() != "unstructured") + { + CONDUIT_ERROR("The tile topology must be 2D."); + } + + // Compute the tile extents. + if(t.has_path("translate/x")) + m_width = t["translate/x"].to_double(); + else + { + const auto &xc = getCoordset().fetch_existing("values/x").as_double_array(); + m_width = xc.max() - xc.min(); + } + + if(t.has_path("translate/y")) + m_height = t["translate/y"].to_double(); + else + { + const auto &yc = getCoordset().fetch_existing("values/y").as_double_array(); + m_height = yc.max() - yc.min(); + } +} + +//--------------------------------------------------------------------------- +/** + \brief Generate coordinate and connectivity arrays using a tiled mesh pattern, + given by the Tile class. + + \param nx The number of tiles in the X dimension. + \param ny The number of tiles in the Y dimension. + \param nz The number of tiles in the Z dimension. + \param[out] res The output node. + \param options A node that may contain additional control options. + */ +void +Tiler::generate(conduit::index_t nx, conduit::index_t ny, conduit::index_t nz, + conduit::Node &res, + const conduit::Node &options) +{ + std::vector x, y, z; + std::vector conn, sizes, bconn, bsizes; + std::vector btype; + + // Process any options. + if(options.has_path("tile")) + initialize(options.fetch_existing("tile")); + + bool reorder = true; + if(options.has_path("reorder")) + reorder = options.fetch_existing("reorder").as_string() == "kdtree"; + + if(options.has_path("meshname")) + meshName = options.fetch_existing("meshname").as_string(); + if(options.has_path("boundarymeshname")) + boundaryMeshName = options.fetch_existing("boundarymeshname").as_string(); + + conduit::DataType indexDT(conduit::DataType::index_t()); + if(options.has_child("datatype")) + { + auto s = options.fetch_existing("datatype").as_string(); + if((s == "int") || (s == "int32") || (s == "integer")) + { + indexDT = conduit::DataType::int32(); + } + } + + // Make a transformation matrix for the tile points. + double origin[] = {0., 0., 0.}; + double tx = width(), ty = height(); + double z1 = std::max(width(), height()) * nz; + double M[3][3] = {{1., 0., 0.}, {0., 1., 0.}, {0., 0., 1.}}; + if(options.has_path("extents")) + { + // Extents for the domain were given. Fit the domain into it. + auto extents = options.fetch_existing("extents").as_double_accessor(); + tx = (extents[1] - extents[0]) / nx; + ty = (extents[3] - extents[2]) / ny; + origin[0] = extents[0]; + origin[1] = extents[2]; + origin[2] = extents[4]; + z1 = extents[5]; + } + else + { + // There are no extents so figure out some based on the domains, if present. + if(options.has_path("domain") && options.has_path("domains")) + { + auto domain = options.fetch_existing("domain").as_int_accessor(); + auto domains = options.fetch_existing("domains").as_int_accessor(); + if(domain.number_of_elements() == 3 && + domain.number_of_elements() == domains.number_of_elements()) + { + origin[0] = domain[0] * nx * width(); + origin[1] = domain[1] * ny * height(); + origin[2] = domain[2] * z1; + z1 = origin[2] + z1; + } + } + } + // Scaling + M[0][0] = tx / width(); + M[1][1] = ty / height(); + // Translation + M[2][0] = origin[0]; + M[2][1] = origin[1]; + + // Number of tile points. + const auto nTilePts = getCoordset().fetch_existing("values/x").dtype().number_of_elements(); + + // Make a pass where we make nx*ny tiles so we can generate their points. + std::vector tiles(nx * ny); + for(conduit::index_t j = 0; j < ny; j++) + { + M[2][0] = origin[0]; + for(conduit::index_t i = 0; i < nx; i++) + { + Tile ¤t = tiles[(j*nx + i)]; + + // The first time we've used the tile, set its size. + current.reset(nTilePts); + + // Copy some previous points over so they can be shared. + if(i > 0) + { + Tile &prevX = tiles[(j*nx + i - 1)]; + current.setPointIds(left(), prevX.getPointIds(right())); + } + if(j > 0) + { + Tile &prevY = tiles[((j-1)*nx + i)]; + current.setPointIds(bottom(), prevY.getPointIds(top())); + } + + addPoints(M, current.getPointIds(), x, y); + M[2][0] += tx; + } + M[2][1] += ty; + } + + conduit::index_t ptsPerPlane = 0; + std::string shape2, shape3; + shape2 = getTopology().fetch_existing("elements/shape").as_string(); + if(nz < 1) + { + // Iterate over the tiles and add their quads. + // TODO: reserve size for conn, sizes + for(conduit::index_t j = 0; j < ny; j++) + { + for(conduit::index_t i = 0; i < nx; i++) + { + Tile ¤t = tiles[(j*nx + i)]; + iterateFaces(current.getPointIds(), 0, false, BoundaryBack, + [&](const conduit::index_t *ids, conduit::index_t npts, int) + { + for(conduit::index_t pi = 0; pi < npts; pi++) + conn.push_back(ids[pi]); + sizes.push_back(npts); + }); + } + } + // NOTE: z coords in output will be empty. + } + else + { + ptsPerPlane = static_cast(x.size()); + + shape3 = (shape2 == "tri") ? "wedge" : "hex"; + + // We have x,y points now. We need to replicate them to make multiple planes. + // We make z coordinates too. + conduit::index_t nplanes = nz + 1; + x.reserve(ptsPerPlane * nplanes); + y.reserve(ptsPerPlane * nplanes); + z.reserve(ptsPerPlane * nplanes); + for(conduit::index_t i = 0; i < ptsPerPlane; i++) + z.push_back(origin[2]); + for(conduit::index_t p = 1; p < nplanes; p++) + { + double t = static_cast(p) / static_cast(nplanes - 1); + double zvalue = (1. - t) * origin[2] + t * z1; + for(conduit::index_t i = 0; i < ptsPerPlane; i++) + { + x.push_back(x[i]); + y.push_back(y[i]); + z.push_back(zvalue); + } + } + + // Iterate over the tiles and add their hexs. + // TODO: reserve size for conn, sizes + for(conduit::index_t k = 0; k < nz; k++) + { + conduit::index_t offset1 = k * ptsPerPlane; + conduit::index_t offset2 = offset1 + ptsPerPlane; + + for(conduit::index_t j = 0; j < ny; j++) + { + for(conduit::index_t i = 0; i < nx; i++) + { + Tile ¤t = tiles[(j*nx + i)]; + addVolumeElements(current.getPointIds(), offset1, offset2, conn, sizes); + } + } + } + } + + // Make the Blueprint mesh. + res["coordsets/coords/type"] = "explicit"; + res["coordsets/coords/values/x"].set(x); + res["coordsets/coords/values/y"].set(y); + if(!z.empty()) + res["coordsets/coords/values/z"].set(z); + + conduit::Node &topo = res["topologies/" + meshName]; + topo["type"] = "unstructured"; + topo["coordset"] = "coords"; + topo["elements/shape"] = z.empty() ? shape2 : shape3; + conduit::Node tmp; + tmp.set_external(conn.data(), conn.size()); + tmp.to_data_type(indexDT.id(), topo["elements/connectivity"]); + tmp.set_external(sizes.data(), sizes.size()); + tmp.to_data_type(indexDT.id(), topo["elements/sizes"]); + +#ifdef CONDUIT_TILER_DEBUG_FIELDS + // Add fields to test the reordering. + std::vector nodeids, elemids; + auto npts = static_cast(x.size()); + nodeids.reserve(npts); + for(conduit::index_t i = 0; i < npts; i++) + nodeids.push_back(i); + res["fields/nodeids/topology"] = meshName; + res["fields/nodeids/association"] = "vertex"; + res["fields/nodeids/values"].set(nodeids); + + auto nelem = static_cast(sizes.size()); + elemids.reserve(nelem); + for(conduit::index_t i = 0; i < nelem; i++) + elemids.push_back(i); + res["fields/elemids/topology"] = meshName; + res["fields/elemids/association"] = "element"; + res["fields/elemids/values"].set(elemids); + + std::vector dist; + dist.reserve(npts); + if(nz < 1) + { + for(conduit::index_t i = 0; i < npts; i++) + dist.push_back(sqrt(x[i]*x[i] + y[i]*y[i])); + } + else + { + for(conduit::index_t i = 0; i < npts; i++) + dist.push_back(sqrt(x[i]*x[i] + y[i]*y[i] + z[i]*z[i])); + } + res["fields/dist/topology"] = meshName; + res["fields/dist/association"] = "vertex"; + res["fields/dist/values"].set(dist); +#endif + + // Reorder the elements unless it was turned off. + std::vector old2NewPoint; + if(reorder) + { + // 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); + +#ifdef CONDUIT_USE_PARTITIONER_FOR_REORDER + // NOTE: This was an idea I had after I made reorder. Reordering is like + // making an explicit selection for the partitioner. Can we just use + // the partitioner? Kind of, it turns out. + // + // 1. Elements are reordered like we want + // 2. Nodes are not reordered in their order of use by elements. + // 3. Passing the same node as input/output does bad things. + + // Make an explicit selection for partition to do the reordering. + conduit::Node options; + conduit::Node &sel = options["selections"].append(); + sel["type"] = "explicit"; + sel["topology"] = meshName; + sel["elements"].set_external(const_cast(elemOrder.data()), elemOrder.size()); + conduit::Node output; + conduit::blueprint::mesh::partition(res, options, output); + + // Extract the vertex mapping. + auto ids = output.fetch_existing("fields/original_vertex_ids/values/ids").as_index_t_accessor(); + for(conduit::index_t i = 0; i < ids.number_of_elements(); i++) + { + old2NewPoint.push_back(ids[i]); + } + res.reset(); + res.move(output); +#else + conduit::blueprint::mesh::utils::topology::unstructured::reorder( + topo, res["coordsets/coords"], res["fields"], + elemOrder, + topo, res["coordsets/coords"], res["fields"], + old2NewPoint); +#endif + } + + // Boundaries + std::string bshape; + bool flags[6]; + boundaryFlags(options, flags); + + if(nz < 1) + { + // 2D + bshape = "line"; + if(reorder) + { + iterateBoundary2D(tiles, nx, ny, flags, + [&](const conduit::index_t *ids, conduit::index_t npts, int bnd) + { + for(conduit::index_t i = 0; i < npts; i++) + bconn.push_back(old2NewPoint[ids[i]]); // Renumber + bsizes.push_back(npts); + btype.push_back(bnd + 1); // Make 1-origin + }); + } + else + { + iterateBoundary2D(tiles, nx, ny, flags, + [&](const conduit::index_t *ids, conduit::index_t npts, int bnd) + { + for(conduit::index_t i = 0; i < npts; i++) + bconn.push_back(ids[i]); + bsizes.push_back(npts); + btype.push_back(bnd + 1); // Make 1-origin + }); + } + } + else + { + // 3D + bshape = "quad"; + bool anyNonQuads = false; + if(reorder) + { + iterateBoundary3D(tiles, nx, ny, nz, ptsPerPlane, flags, + [&](const conduit::index_t *ids, conduit::index_t npts, int bnd) + { + for(conduit::index_t i = 0; i < npts; i++) + bconn.push_back(old2NewPoint[ids[i]]); // Renumber + bsizes.push_back(npts); + btype.push_back(bnd + 1); // Make 1-origin + anyNonQuads |= (npts != 4); + }); + } + else + { + iterateBoundary3D(tiles, nx, ny, nz, ptsPerPlane, flags, + [&](const conduit::index_t *ids, conduit::index_t npts, int bnd) + { + for(conduit::index_t i = 0; i < npts; i++) + bconn.push_back(ids[i]); + bsizes.push_back(npts); + btype.push_back(bnd + 1); // Make 1-origin + anyNonQuads |= (npts != 4); + }); + } + if(anyNonQuads) + bshape = "polygonal"; + } + if(!bconn.empty()) + { + conduit::Node &btopo = res["topologies/" + boundaryMeshName]; + btopo["type"] = "unstructured"; + btopo["coordset"] = "coords"; + btopo["elements/shape"] = bshape; + + tmp.set_external(bconn.data(), bconn.size()); + tmp.to_data_type(indexDT.id(), btopo["elements/connectivity"]); + + tmp.set_external(bsizes.data(), bsizes.size()); + tmp.to_data_type(indexDT.id(), btopo["elements/sizes"]); + + if(bshape == "polygonal") + conduit::blueprint::mesh::utils::topology::unstructured::generate_offsets(btopo, btopo["elements/offsets"]); + + res["fields/boundary_attribute/topology"] = boundaryMeshName; + res["fields/boundary_attribute/association"] = "element"; + res["fields/boundary_attribute/values"].set(btype); + } + + // Build an adjacency set. + addAdjset(tiles, nx, ny, nz, ptsPerPlane, reorder, old2NewPoint, options, res); + +#if 0 + // Print for debugging. + conduit::Node opts; + opts["num_children_threshold"] = 100000; + opts["num_elements_threshold"] = 500; + std::cout << res.to_summary_string(opts) << std::endl; +#endif +} + +//--------------------------------------------------------------------------- +void +Tiler::boundaryFlags(const conduit::Node &options, bool flags[6]) const +{ + bool handled = false; + if(options.has_path("domain") && options.has_path("domains")) + { + auto domain = options.fetch_existing("domain").as_int_accessor(); + auto domains = options.fetch_existing("domains").as_int_accessor(); + if(domain.number_of_elements() == 3 && + domain.number_of_elements() == domains.number_of_elements()) + { + int ndoms = domains[0] * domains[1] * domains[2]; + if(ndoms > 1) + { + flags[BoundaryLeft] = (domain[0] == 0); + flags[BoundaryRight] = (domain[0] == domains[0]-1); + flags[BoundaryBottom] = (domain[1] == 0); + flags[BoundaryTop] = (domain[1] == domains[1]-1); + flags[BoundaryBack] = (domain[2] == 0); + flags[BoundaryFront] = (domain[2] == domains[2]-1); + + handled = true; + } + } + } + if(!handled) + { + for(int i = 0; i < 6; i++) + flags[i] = true; + } +} + +//--------------------------------------------------------------------------- +template +void +Tiler::iterateBoundary2D(const std::vector &tiles, + conduit::index_t nx, + conduit::index_t ny, + const bool flags[6], + Body &&body) const +{ + conduit::index_t idlist[2]; + + if(flags[BoundaryLeft]) + { + for(conduit::index_t i = 0, j = 0; j < ny; j++) + { + const Tile ¤t = tiles[(j*nx + i)]; + const auto ids = current.getPointIds(left()); + for(size_t bi = 0; bi < ids.size() - 1; bi++) + { + idlist[0] = ids[bi]; + idlist[1] = ids[bi + 1]; + body(idlist, 2, BoundaryLeft); + } + } + } + if(flags[BoundaryRight]) + { + for(conduit::index_t i = nx - 1, j = 0; j < ny; j++) + { + const Tile ¤t = tiles[(j*nx + i)]; + const auto ids = current.getPointIds(right()); + for(size_t bi = 0; bi < ids.size() - 1; bi++) + { + idlist[0] = ids[bi]; + idlist[1] = ids[bi + 1]; + body(idlist, 2, BoundaryRight); + } + } + } + if(flags[BoundaryBottom]) + { + for(conduit::index_t i = 0, j = 0; i < nx; i++) + { + const Tile ¤t = tiles[(j*nx + i)]; + const auto ids = current.getPointIds(bottom()); + for(size_t bi = 0; bi < ids.size() - 1; bi++) + { + idlist[0] = ids[bi]; + idlist[1] = ids[bi + 1]; + body(idlist, 2, BoundaryBottom); + } + } + } + if(flags[BoundaryTop]) + { + for(conduit::index_t i = 0, j = ny - 1; i < nx; i++) + { + const Tile ¤t = tiles[(j*nx + i)]; + const auto ids = current.getPointIds(top()); + for(size_t bi = 0; bi < ids.size() - 1; bi++) + { + idlist[0] = ids[bi]; + idlist[1] = ids[bi + 1]; + body(idlist, 2, BoundaryTop); + } + } + } +} + +//--------------------------------------------------------------------------- +template +void +Tiler::iterateBoundary3D(const std::vector &tiles, + conduit::index_t nx, + conduit::index_t ny, + conduit::index_t nz, + conduit::index_t nPtsPerPlane, + const bool flags[6], + Body &&body) const +{ + conduit::index_t idlist[4]; + + if(flags[BoundaryLeft]) + { + for(conduit::index_t k = 0; k < nz; k++) + { + conduit::index_t offset1 = k * nPtsPerPlane; + conduit::index_t offset2 = (k + 1) * nPtsPerPlane; + for(conduit::index_t i = 0, j = 0; j < ny; j++) + { + const Tile ¤t = tiles[(j*nx + i)]; + const auto ids = current.getPointIds(left()); + for(size_t bi = 0; bi < ids.size() - 1; bi++) + { + idlist[0] = offset1 + ids[bi]; + idlist[1] = offset2 + ids[bi]; + idlist[2] = offset2 + ids[bi + 1]; + idlist[3] = offset1 + ids[bi + 1]; + body(idlist, 4, BoundaryLeft); + } + } + } + } + if(flags[BoundaryRight]) + { + for(conduit::index_t k = 0; k < nz; k++) + { + conduit::index_t offset1 = k * nPtsPerPlane; + conduit::index_t offset2 = (k + 1) * nPtsPerPlane; + for(conduit::index_t i = nx - 1, j = 0; j < ny; j++) + { + const Tile ¤t = tiles[(j*nx + i)]; + const auto ids = current.getPointIds(right()); + for(size_t bi = 0; bi < ids.size() - 1; bi++) + { + idlist[0] = offset1 + ids[bi]; + idlist[1] = offset1 + ids[bi + 1]; + idlist[2] = offset2 + ids[bi + 1]; + idlist[3] = offset2 + ids[bi]; + body(idlist, 4, BoundaryRight); + } + } + } + } + if(flags[BoundaryBottom]) + { + for(conduit::index_t k = 0; k < nz; k++) + { + conduit::index_t offset1 = k * nPtsPerPlane; + conduit::index_t offset2 = (k + 1) * nPtsPerPlane; + for(conduit::index_t i = 0, j = 0; i < nx; i++) + { + const Tile ¤t = tiles[(j*nx + i)]; + const auto ids = current.getPointIds(bottom()); + for(size_t bi = 0; bi < ids.size() - 1; bi++) + { + idlist[0] = offset1 + ids[bi]; + idlist[1] = offset1 + ids[bi + 1]; + idlist[2] = offset2 + ids[bi + 1]; + idlist[3] = offset2 + ids[bi]; + body(idlist, 4, BoundaryBottom); + } + } + } + } + if(flags[BoundaryTop]) + { + for(conduit::index_t k = 0; k < nz; k++) + { + conduit::index_t offset1 = k * nPtsPerPlane; + conduit::index_t offset2 = (k + 1) * nPtsPerPlane; + for(conduit::index_t i = 0, j = ny - 1; i < nx; i++) + { + const Tile ¤t = tiles[(j*nx + i)]; + const auto ids = current.getPointIds(top()); + for(size_t bi = 0; bi < ids.size() - 1; bi++) + { + idlist[0] = offset1 + ids[bi]; + idlist[1] = offset2 + ids[bi]; + idlist[2] = offset2 + ids[bi + 1]; + idlist[3] = offset1 + ids[bi + 1]; + body(idlist, 4, BoundaryTop); + } + } + } + } + if(flags[BoundaryBack]) + { + for(conduit::index_t j = 0; j < ny; j++) + for(conduit::index_t i = 0; i < nx; i++) + { + const Tile ¤t = tiles[(j*nx + i)]; + iterateFaces(current.getPointIds(), 0, true, BoundaryBack, body); + } + } + if(flags[BoundaryFront]) + { + for(conduit::index_t j = 0; j < ny; j++) + for(conduit::index_t i = 0; i < nx; i++) + { + const Tile ¤t = tiles[(j*nx + i)]; + iterateFaces(current.getPointIds(), nz * nPtsPerPlane, false, BoundaryFront, body); + } + } +} + +//--------------------------------------------------------------------------- +conduit::index_t +Tiler::domainIndex(conduit::index_t d0, conduit::index_t d1, conduit::index_t d2, + const std::vector &domain, + const std::vector &domains) const +{ + const conduit::index_t dom[3] = {domain[0] + d0, domain[1] + d1, domain[2] + d2}; + // If the domain exists, make its domain id. + conduit::index_t domainId = InvalidDomain; + if((dom[0] >= 0 && dom[0] < domains[0]) && + (dom[1] >= 0 && dom[1] < domains[1]) && + (dom[2] >= 0 && dom[2] < domains[2])) + { + domainId = dom[2] * (domains[0] * domains[1]) + dom[1] * domains[0] + dom[0]; + } + return domainId; +} + +//--------------------------------------------------------------------------- +conduit::index_t +Tiler::domainIndex(const conduit::index_t delta[3], + const std::vector &domain, + const std::vector &domains) const +{ + return domainIndex(delta[0], delta[1], delta[2], domain, domains); +} + +//--------------------------------------------------------------------------- +void +Tiler::addAdjsetValues(const std::vector &ptids, + conduit::index_t planeOffset, + bool reorder, + const std::vector &old2NewPoint, + conduit::Node &values) const +{ + if(reorder) + { + values.set(conduit::DataType::index_t(ptids.size())); + conduit::index_t *dest = values.as_index_t_ptr(); + if(planeOffset > 0) + { + for(size_t i = 0; i < ptids.size(); i++) + *dest++ = old2NewPoint[ptids[i] + planeOffset]; + } + else + { + for(size_t i = 0; i < ptids.size(); i++) + *dest++ = old2NewPoint[ptids[i]]; + } + } + else if(planeOffset > 0) + { + values.set(conduit::DataType::index_t(ptids.size())); + conduit::index_t *dest = values.as_index_t_ptr(); + for(size_t i = 0; i < ptids.size(); i++) + *dest++ = ptids[i] + planeOffset; + } + else + { + values.set(ptids); + } +} + +//--------------------------------------------------------------------------- +void +Tiler::addAdjset(const std::vector &tiles, + conduit::index_t nx, + conduit::index_t ny, + conduit::index_t nz, + conduit::index_t ptsPerPlane, + bool reorder, + const std::vector &old2NewPoint, + const conduit::Node &options, + conduit::Node &out) const +{ + // Make the adjset name for 2 domains. + auto adjset_name = [](conduit::index_t d0, conduit::index_t d1) { + if(d0 > d1) + std::swap(d0, d1); + std::stringstream ss; + ss << "domain_" << d0 << "_" << d1; + return ss.str(); + }; + + // Build up unique points, given a list of points. Use unique and ptvec by capture. + std::set unique; + std::vector ptvec; + auto addPoints = [&](const conduit::index_t *ids, conduit::index_t npts, int /*bnd*/) + { + for(conduit::index_t i = 0; i < npts; i++) + { + const auto id = ids[i]; + if(unique.find(id) == unique.end()) + { + unique.insert(id); + ptvec.push_back(id); + } + } + }; + + // We need to know where this domain is in the domains to make the adjset. + if(options.has_child("domain") && options.has_child("domains")) + { + auto domain = toIndexVector(options.fetch_existing("domain")); + auto domains = toIndexVector(options.fetch_existing("domains")); + if(domain.size() == 3 && domain.size() == domains.size()) + { + if(domains[0] * domains[1] * domains[2] > 1) + { + auto thisDom = domainIndex(0, 0, 0, domain, domains); + + // Make a state node. + out["state/domain_id"] = thisDom; + + // Make the top level adjset nodes. + conduit::Node &adjset = out["adjsets/" + meshName + "_adjset"]; + adjset["association"] = "vertex"; + adjset["topology"] = meshName; + conduit::Node &groups = adjset["groups"]; + + // Trace around the tile to get the points for the edges of the tile + // in original node order. + //----------------------------------------------------------- + std::vector left, right, bottom, top; + const bool left_flags[] = {true, false, false, false}; + iterateBoundary2D(tiles, nx, ny, left_flags, addPoints); + left.swap(ptvec); // Steal the points. + unique.clear(); + const bool right_flags[] = {false, true, false, false}; + iterateBoundary2D(tiles, nx, ny, right_flags, addPoints); + right.swap(ptvec); // Steal the points. + unique.clear(); + const bool bottom_flags[] = {false, false, true, false}; + iterateBoundary2D(tiles, nx, ny, bottom_flags, addPoints); + bottom.swap(ptvec); // Steal the points. + unique.clear(); + const bool top_flags[] = {false, false, false, true}; + iterateBoundary2D(tiles, nx, ny, top_flags, addPoints); + top.swap(ptvec); // Steal the points. + unique.clear(); + + // Make corner neighbors. + //----------------------------------------------------------- + conduit::index_t frontPlane = ptsPerPlane * nz; + conduit::index_t corner[8]; + corner[0] = left[0]; // back, lower left + corner[1] = right[0]; // back, lower right + corner[2] = left[left.size() - 1]; // back, upper left + corner[3] = right[right.size() - 1]; // back, upper right + corner[4] = left[0] + frontPlane; // front, lower left + corner[5] = right[0] + frontPlane; // front, lower right + corner[6] = left[left.size() - 1] + frontPlane; // front, upper left + corner[7] = right[right.size() - 1] + frontPlane; // front, upper right + int maxCornerNeighbors = (nz < 1) ? 4 : 8; + conduit::index_t z0 = (nz < 1) ? 0 : -1; + conduit::index_t z1 = (nz < 1) ? 0 : 1; + conduit::index_t neighborId = InvalidDomain; + const conduit::index_t cornerNeighbors[][3] = { + {-1, -1, z0}, + {1, -1, z0}, + {-1, 1, z0}, + {1, 1, z0}, + {-1, -1, z1}, + {1, -1, z1}, + {-1, 1, z1}, + {1, 1, z1} + }; + for(int ni = 0; ni < maxCornerNeighbors; ni++) + { + if((neighborId = domainIndex(cornerNeighbors[ni], domain, domains)) != InvalidDomain) + { + auto name = adjset_name(thisDom, neighborId); + conduit::Node &group = groups[name]; + group["neighbors"] = neighborId; + group["values"] = reorder ? old2NewPoint[corner[ni]] : corner[ni]; + } + } + + // Make edge neighbors for 3D + //----------------------------------------------------------- + if(nz > 0) + { + // Back, left edge. + if((neighborId = domainIndex(-1, 0, -1, domain, domains)) != InvalidDomain) + { + auto name = adjset_name(thisDom, neighborId); + conduit::Node &group = groups[name]; + group["neighbors"] = neighborId; + addAdjsetValues(left, 0, reorder, old2NewPoint, group["values"]); + } + // Back, right edge. + if((neighborId = domainIndex(1, 0, -1, domain, domains)) != InvalidDomain) + { + auto name = adjset_name(thisDom, neighborId); + conduit::Node &group = groups[name]; + group["neighbors"] = neighborId; + addAdjsetValues(right, 0, reorder, old2NewPoint, group["values"]); + } + // Back, bottom edge. + if((neighborId = domainIndex(0, -1, -1, domain, domains)) != InvalidDomain) + { + auto name = adjset_name(thisDom, neighborId); + conduit::Node &group = groups[name]; + group["neighbors"] = neighborId; + addAdjsetValues(bottom, 0, reorder, old2NewPoint, group["values"]); + } + // Back, top edge. + if((neighborId = domainIndex(0, 1, -1, domain, domains)) != InvalidDomain) + { + auto name = adjset_name(thisDom, neighborId); + conduit::Node &group = groups[name]; + group["neighbors"] = neighborId; + addAdjsetValues(top, 0, reorder, old2NewPoint, group["values"]); + } + // Lower left edge + if((neighborId = domainIndex(-1, -1, 0, domain, domains)) != InvalidDomain) + { + auto name = adjset_name(thisDom, neighborId); + std::vector values; + values.reserve(nz); + for(conduit::index_t zi = 0; zi <= nz; zi++) + values.push_back(corner[0] + zi * ptsPerPlane); + conduit::Node &group = groups[name]; + group["neighbors"] = neighborId; + addAdjsetValues(values, 0, reorder, old2NewPoint, group["values"]); + } + // Lower right edge + if((neighborId = domainIndex(1, -1, 0, domain, domains)) != InvalidDomain) + { + auto name = adjset_name(thisDom, neighborId); + std::vector values; + values.reserve(nz); + for(conduit::index_t zi = 0; zi <= nz; zi++) + values.push_back(corner[1] + zi * ptsPerPlane); + conduit::Node &group = groups[name]; + group["neighbors"] = neighborId; + addAdjsetValues(values, 0, reorder, old2NewPoint, group["values"]); + } + // Upper left edge + if((neighborId = domainIndex(-1, 1, 0, domain, domains)) != InvalidDomain) + { + auto name = adjset_name(thisDom, neighborId); + std::vector values; + values.reserve(nz); + for(conduit::index_t zi = 0; zi <= nz; zi++) + values.push_back(corner[2] + zi * ptsPerPlane); + conduit::Node &group = groups[name]; + group["neighbors"] = neighborId; + addAdjsetValues(values, 0, reorder, old2NewPoint, group["values"]); + } + // Upper right edge + if((neighborId = domainIndex(1, 1, 0, domain, domains)) != InvalidDomain) + { + auto name = adjset_name(thisDom, neighborId); + std::vector values; + values.reserve(nz); + for(conduit::index_t zi = 0; zi <= nz; zi++) + values.push_back(corner[3] + zi * ptsPerPlane); + conduit::Node &group = groups[name]; + group["neighbors"] = neighborId; + addAdjsetValues(values, 0, reorder, old2NewPoint, group["values"]); + } + // Front, left edge. + if((neighborId = domainIndex(-1, 0, 1, domain, domains)) != InvalidDomain) + { + auto name = adjset_name(thisDom, neighborId); + conduit::Node &group = groups[name]; + group["neighbors"] = neighborId; + addAdjsetValues(left, frontPlane, reorder, old2NewPoint, group["values"]); + } + // Front, right edge. + if((neighborId = domainIndex(1, 0, 1, domain, domains)) != InvalidDomain) + { + auto name = adjset_name(thisDom, neighborId); + conduit::Node &group = groups[name]; + group["neighbors"] = neighborId; + addAdjsetValues(right, frontPlane, reorder, old2NewPoint, group["values"]); + } + // Front, bottom edge. + if((neighborId = domainIndex(0, -1, 1, domain, domains)) != InvalidDomain) + { + auto name = adjset_name(thisDom, neighborId); + conduit::Node &group = groups[name]; + group["neighbors"] = neighborId; + addAdjsetValues(bottom, frontPlane, reorder, old2NewPoint, group["values"]); + } + // Front, top edge. + if((neighborId = domainIndex(0, 1, 1, domain, domains)) != InvalidDomain) + { + auto name = adjset_name(thisDom, neighborId); + conduit::Node &group = groups[name]; + group["neighbors"] = neighborId; + addAdjsetValues(top, frontPlane, reorder, old2NewPoint, group["values"]); + } + } + + // Make "face" neighbors. + //----------------------------------------------------------- + conduit::index_t neighbor[6]; + neighbor[BoundaryLeft] = domainIndex(-1, 0, 0, domain, domains); + neighbor[BoundaryRight] = domainIndex(1, 0, 0, domain, domains); + neighbor[BoundaryBottom] = domainIndex(0, -1, 0, domain, domains); + neighbor[BoundaryTop] = domainIndex(0, 1, 0, domain, domains); + neighbor[BoundaryBack] = domainIndex(0, 0, -1, domain, domains); + neighbor[BoundaryFront] = domainIndex(0, 0, 1, domain, domains); + int maxNeighbors = (nz < 1) ? 4 : 6; + for(int ni = 0; ni < maxNeighbors; ni++) + { + // If this domain has no neighbor in the current direction, skip. + if(neighbor[ni] == InvalidDomain) + continue; + + // Iterate over faces and come up with unique points. + bool flags[6] = {false, false, false, false, false, false}; + flags[ni] = true; + unique.clear(); + ptvec.clear(); + if(nz < 1) + iterateBoundary2D(tiles, nx, ny, flags, addPoints); + else + iterateBoundary3D(tiles, nx, ny, nz, ptsPerPlane, flags, addPoints); + + if(!ptvec.empty()) + { + auto name = adjset_name(thisDom, neighbor[ni]); + conduit::Node &group = groups[name]; + group["neighbors"] = neighbor[ni]; + addAdjsetValues(ptvec, 0, reorder, old2NewPoint, group["values"]); + } + } // end for + } + } + } +} + +} +//----------------------------------------------------------------------------- +// -- end conduit::blueprint::mesh::examples::detail -- +//----------------------------------------------------------------------------- + +void +tiled(conduit::index_t nx, conduit::index_t ny, conduit::index_t nz, + conduit::Node &res, const conduit::Node &options) +{ + detail::Tiler T; + T.generate(nx, ny, nz, res, options); +} + +} +//----------------------------------------------------------------------------- +// -- end conduit::blueprint::mesh::examples -- +//----------------------------------------------------------------------------- + +} +//----------------------------------------------------------------------------- +// -- end conduit::blueprint::mesh -- +//----------------------------------------------------------------------------- + +} +//----------------------------------------------------------------------------- +// -- end conduit::blueprint -- +//----------------------------------------------------------------------------- + +} +//----------------------------------------------------------------------------- +// -- end conduit -- +//----------------------------------------------------------------------------- + + + diff --git a/src/libs/blueprint/conduit_blueprint_mesh_examples_tiled.hpp b/src/libs/blueprint/conduit_blueprint_mesh_examples_tiled.hpp new file mode 100644 index 000000000..1d16d5847 --- /dev/null +++ b/src/libs/blueprint/conduit_blueprint_mesh_examples_tiled.hpp @@ -0,0 +1,78 @@ +// Copyright (c) Lawrence Livermore National Security, LLC and other Conduit +// Project developers. See top-level LICENSE AND COPYRIGHT files for dates and +// other details. No copyright assignment is required to contribute to Conduit. + +//----------------------------------------------------------------------------- +/// +/// file: conduit_blueprint_mesh_examples_tiled.hpp +/// +//----------------------------------------------------------------------------- + +#ifndef CONDUIT_BLUEPRINT_MESH_EXAMPLES_TILED_HPP +#define CONDUIT_BLUEPRINT_MESH_EXAMPLES_TILED_HPP + +//----------------------------------------------------------------------------- +// conduit lib includes +//----------------------------------------------------------------------------- +#include "conduit.hpp" +#include "conduit_blueprint.hpp" +#include "conduit_blueprint_exports.h" + +//----------------------------------------------------------------------------- +// -- begin conduit::-- +//----------------------------------------------------------------------------- +namespace conduit +{ + +//----------------------------------------------------------------------------- +// -- begin conduit::blueprint -- +//----------------------------------------------------------------------------- +namespace blueprint +{ + +//----------------------------------------------------------------------------- +// -- begin conduit::blueprint::mesh -- +//----------------------------------------------------------------------------- +namespace mesh +{ + +//----------------------------------------------------------------------------- +/// Methods that generate example meshes. +//----------------------------------------------------------------------------- +namespace examples +{ + /// Generates a tiled unstructured mesh of quads or hexs. + void CONDUIT_BLUEPRINT_API tiled(conduit::index_t nx, + conduit::index_t ny, + conduit::index_t nz, + conduit::Node &res, + const conduit::Node &options); + +} +//----------------------------------------------------------------------------- +// -- end conduit::blueprint::mesh::examples -- +//----------------------------------------------------------------------------- + + +//----------------------------------------------------------------------------- +} +//----------------------------------------------------------------------------- +// -- end conduit::blueprint::mesh -- +//----------------------------------------------------------------------------- + + +} +//----------------------------------------------------------------------------- +// -- end conduit::blueprint -- +//----------------------------------------------------------------------------- + +} +//----------------------------------------------------------------------------- +// -- end conduit -- +//----------------------------------------------------------------------------- + + +#endif + + + diff --git a/src/libs/blueprint/conduit_blueprint_mesh_kdtree.hpp b/src/libs/blueprint/conduit_blueprint_mesh_kdtree.hpp index 615b0275a..a88136416 100644 --- a/src/libs/blueprint/conduit_blueprint_mesh_kdtree.hpp +++ b/src/libs/blueprint/conduit_blueprint_mesh_kdtree.hpp @@ -61,19 +61,19 @@ namespace utils @tparam Indexable The container that contains elements of type T. This can be raw pointers, a vector, a data_array, data_accessor, etc. - @tparam T The storage type of the coordinates. + @tparam CoordinateType The storage type of the coordinates. @tparam NDIMS The number of dimensions in a coordinate. */ -template +template class kdtree { public: - static const int NoChild; - static const int NotFound; + static const conduit::index_t NoChild; + static const conduit::index_t NotFound; using IndexableType = Indexable; - using BoxType = T[NDIMS][2]; - using PointType = T[NDIMS]; + using BoxType = CoordinateType[NDIMS][2]; + using PointType = CoordinateType[NDIMS]; /// Constructor kdtree(); @@ -82,7 +82,7 @@ class kdtree @brief Set the point tolerance, a distance under which points are the same. @param tolerance The point tolerance distance. */ - void setPointTolerance(T tolerance); + void setPointTolerance(CoordinateType tolerance); /** @brief Initialize the tree so it can look up points. @@ -98,7 +98,7 @@ class kdtree @return The point id of the point or NotFound if the point was not found. */ - int findPoint(const PointType &pt) const; + conduit::index_t findPoint(const PointType &pt) const; /** @brief Return the number of dimensions. @@ -106,6 +106,18 @@ class kdtree */ inline int dims() const { return NDIMS; } + /** + @brief Return the sorted index order. + @return A vector sorted spatially. The elements contain the original indices. + */ + std::vector &getIndices() { return index; } + + /** + @brief Return the sorted index order. + @return A vector sorted spatially. The elements contain the original indices. + */ + const std::vector &getIndices() const { return index; } + /** @brief Print the tree to a stream. @param os A stream to use for printing. @@ -115,8 +127,8 @@ class kdtree /// Represents ranges of the index vector. struct RangeType { - int offset; - int size; + conduit::index_t offset; + conduit::index_t size; }; /// Represents a box in the tree hierarchy. @@ -133,68 +145,68 @@ class kdtree // Internal helpers. void cutBox(const BoxType &input, int dimension, BoxType &A, BoxType &B) const; - void cutRange(const RangeType &input, int dimension, T maxValue, RangeType &A, RangeType &B) const; + void cutRange(const RangeType &input, int dimension, CoordinateType maxValue, RangeType &A, RangeType &B) const; void calculateExtents(); void construct(); - void constructBox(int bp, const RangeType &range, const BoxType &b, int level, int maxlevel); + void constructBox(conduit::index_t bp, const RangeType &range, const BoxType &b, int level, int maxlevel); int longest(const BoxType &b) const; void sortIndexRange(const RangeType &range, int dimension); bool pointInBox(const PointType &pt, const BoxType &b) const; - bool pointEqual(const PointType &pt, int index) const; + bool pointEqual(const PointType &pt, conduit::index_t index) const; void printBox(std::ostream &os, const BoxType &b, const std::string &indent) const; private: - std::vector boxes; //!< Describes boxes - std::vector index; //!< Contains sorted coordinate indices - BoxType box; //!< Overall bounding box - Indexable coords[NDIMS]; //!< Coordinate arrays - size_t coordlen; //!< Length of a coordinate array. - T pointTolerance; //!< Distance to a point before it is same. - T pointTolerance2; //!< pointTolerance^2. + std::vector boxes; //!< Describes boxes + std::vector index; //!< Contains sorted coordinate indices + BoxType box; //!< Overall bounding box + Indexable coords[NDIMS]; //!< Coordinate arrays + size_t coordlen; //!< Length of a coordinate array. + CoordinateType pointTolerance; //!< Distance to a point before it is same. + CoordinateType pointTolerance2; //!< pointTolerance^2. }; //--------------------------------------------------------------------------- -template +template std::ostream & -operator << (std::ostream &os, const kdtree &obj) +operator << (std::ostream &os, const kdtree &obj) { obj.print(os); return os; } //--------------------------------------------------------------------------- -template -const int kdtree::NoChild = -1; +template +const conduit::index_t kdtree::NoChild = -1; -template -const int kdtree::NotFound = -1; +template +const conduit::index_t kdtree::NotFound = -1; //--------------------------------------------------------------------------- -template -kdtree::kdtree() : boxes(), index() +template +kdtree::kdtree() : boxes(), index() { for(int i = 0; i < dims(); i++) { - box[i][0] = T{}; - box[i][1] = T{}; + box[i][0] = CoordinateType{}; + box[i][1] = CoordinateType{}; coords[i] = Indexable{}; } coordlen = 0; - constexpr T DEFAULT_POINT_TOLERANCE = (T)(1.e-9); + constexpr auto DEFAULT_POINT_TOLERANCE = static_cast(1.e-9); setPointTolerance(DEFAULT_POINT_TOLERANCE); } //--------------------------------------------------------------------------- -template -void kdtree::setPointTolerance(T tolerance) +template +void kdtree::setPointTolerance(CoordinateType tolerance) { pointTolerance = tolerance; pointTolerance2 = pointTolerance * pointTolerance; } //--------------------------------------------------------------------------- -template -void kdtree::initialize(Indexable c[NDIMS], size_t len) +template +void kdtree::initialize(Indexable c[NDIMS], size_t len) { boxes.clear(); index.clear(); @@ -209,9 +221,9 @@ void kdtree::initialize(Indexable c[NDIMS], size_t len) } //--------------------------------------------------------------------------- -template -void kdtree::printBox(std::ostream &os, - const kdtree::BoxType &b, const std::string &indent) const +template +void kdtree::printBox(std::ostream &os, + const kdtree::BoxType &b, const std::string &indent) const { os << indent << "box: ["; for(int i = 0; i < dims(); i++) @@ -223,8 +235,8 @@ void kdtree::printBox(std::ostream &os, } //--------------------------------------------------------------------------- -template -void kdtree::print(std::ostream &os) const +template +void kdtree::print(std::ostream &os) const { os << "boxesSize: " << boxes.size() << std::endl; os << "boxes:" << std::endl; @@ -270,33 +282,33 @@ void kdtree::print(std::ostream &os) const } //--------------------------------------------------------------------------- -template -void kdtree::cutBox(const BoxType &input, int dimension, +template +void kdtree::cutBox(const BoxType &input, int dimension, BoxType &A, BoxType &B) const { // NOTE: This just does uniform splitting for now. It could sample the // coordinates to get an idea of their center and store a t value // to lerp the dimension min/max. This would help with long dimensions // where points are clustered unevenly. - constexpr T two = 2.; + constexpr CoordinateType two = 2.; memcpy(A, input, sizeof(BoxType)); memcpy(B, input, sizeof(BoxType)); - T mid = (input[dimension][0] + input[dimension][1]) / two; + CoordinateType mid = (input[dimension][0] + input[dimension][1]) / two; A[dimension][1] = mid; B[dimension][0] = mid; } //--------------------------------------------------------------------------- -template -void kdtree::cutRange(const RangeType &input, int dimension, - T maxValue, RangeType &A, RangeType &B) const +template +void kdtree::cutRange(const RangeType &input, int dimension, + CoordinateType maxValue, RangeType &A, RangeType &B) const { // We can't just partition the range in half because it might mean leaving // coordinate data with the same value on either side of the split. // That doesn't work spatially. Instead, partition the range such that A // contains all indices less than maxValue and B contains the rest. - int sizeA = 0; - for(int idx = 0; idx < input.size; idx++) + conduit::index_t sizeA = 0; + for(conduit::index_t idx = 0; idx < input.size; idx++) { if(coords[dimension][index[input.offset + idx]] < maxValue) { @@ -307,7 +319,7 @@ void kdtree::cutRange(const RangeType &input, int dimension break; } } - int sizeB = input.size - sizeA; + conduit::index_t sizeB = input.size - sizeA; A.offset = input.offset; A.size = sizeA; @@ -322,14 +334,14 @@ void kdtree::cutRange(const RangeType &input, int dimension << ", maxValue=" << maxValue << std::endl; std::cout << "\tA={"; - for(int idx = 0; idx < A.size; idx++) + for(conduit::index_t idx = 0; idx < A.size; idx++) { if(idx > 0) std::cout << ", "; std::cout << coords[dimension][index[A.offset + idx]]; } std::cout << "}" << std::endl; std::cout << "\tB={"; - for(int idx = 0; idx < B.size; idx++) + for(conduit::index_t idx = 0; idx < B.size; idx++) { if(idx > 0) std::cout << ", "; std::cout << coords[dimension][index[B.offset + idx]]; @@ -339,11 +351,11 @@ void kdtree::cutRange(const RangeType &input, int dimension } //--------------------------------------------------------------------------- -template -int -kdtree::findPoint(const kdtree::PointType &pt) const +template +conduit::index_t +kdtree::findPoint(const kdtree::PointType &pt) const { - int foundIndex = NotFound; + conduit::index_t foundIndex = NotFound; #ifdef CONDUIT_DEBUG_KDTREE std::cout << "findPoint("; for(int i = 0; i < dims(); i++) @@ -358,7 +370,7 @@ kdtree::findPoint(const kdtree::PointT // Iterate until we find a box that contains the point. BoxType currentBox; memcpy(currentBox, box, sizeof(BoxType)); - int bp = 0, prevbp = 0; + conduit::index_t bp = 0, prevbp = 0; while(boxes[bp].childOffset != NoChild) { // Cut the box to make 2 child boxes. @@ -390,12 +402,12 @@ kdtree::findPoint(const kdtree::PointT std::cout << " Search box: bp=" << prevbp << std::endl; #endif // Check prevbp box for the points of interest. - int offset = boxes[prevbp].range.offset; - int size = boxes[prevbp].range.size; - for(int i = 0; i < size && (foundIndex == NotFound); i++) + conduit::index_t offset = boxes[prevbp].range.offset; + conduit::index_t size = boxes[prevbp].range.size; + for(conduit::index_t i = 0; i < size && (foundIndex == NotFound); i++) { // This the index of the point we want to look at. - int ptIdx = index[offset + i]; + conduit::index_t ptIdx = index[offset + i]; // See if the point is "equal". if(pointEqual(pt, ptIdx)) @@ -418,27 +430,27 @@ kdtree::findPoint(const kdtree::PointT } //--------------------------------------------------------------------------- -template +template bool -kdtree::pointEqual(const kdtree::PointType &pt, int ptIdx) const +kdtree::pointEqual(const kdtree::PointType &pt, conduit::index_t ptIdx) const { - T dist2 = 0.; + CoordinateType dist2 = 0.; for(int i = 0; i < dims(); i++) { - T delta = pt[i] - coords[i][ptIdx]; + CoordinateType delta = pt[i] - coords[i][ptIdx]; dist2 += delta * delta; } return dist2 < pointTolerance2; } //--------------------------------------------------------------------------- -template -void kdtree::calculateExtents() +template +void kdtree::calculateExtents() { for(int i = 0; i < dims(); i++) { - box[i][0] = std::numeric_limits::max(); - box[i][1] = std::numeric_limits::min(); + box[i][0] = std::numeric_limits::max(); + box[i][1] = std::numeric_limits::min(); for(size_t j = 0; j < coordlen; j++) { box[i][0] = std::min(box[i][0], coords[i][j]); @@ -449,15 +461,15 @@ void kdtree::calculateExtents() // Expand the box a little for(int i = 0; i < dims(); i++) { - T d = (box[i][1] - box[i][0]) / static_cast(200.); + CoordinateType d = (box[i][1] - box[i][0]) / static_cast(200.); box[i][0] -= d; box[i][1] += d; } } //--------------------------------------------------------------------------- -template -void kdtree::construct() +template +void kdtree::construct() { // Fill the index with 0,1,2,3,... index.resize(coordlen); @@ -481,15 +493,15 @@ void kdtree::construct() // Build the box for level 0. RangeType range; range.offset = 0; - range.size = static_cast(coordlen); + range.size = static_cast(coordlen); constructBox(0, range, box, 0, maxLevels); } //--------------------------------------------------------------------------- -template -void kdtree::constructBox(int bp, - const kdtree::RangeType &range, - const kdtree::BoxType &b, +template +void kdtree::constructBox(conduit::index_t bp, + const kdtree::RangeType &range, + const kdtree::BoxType &b, int level, int maxLevels) { @@ -504,14 +516,14 @@ void kdtree::constructBox(int bp, BoxType A, B; cutBox(b, boxes[bp].splitDimension, A, B); - T maxValue = A[boxes[bp].splitDimension][1]; + CoordinateType maxValue = A[boxes[bp].splitDimension][1]; RangeType rangeA, rangeB; cutRange(range, boxes[bp].splitDimension, maxValue, rangeA, rangeB); if(level < maxLevels) { - int child0 = boxes[bp].childOffset; - int child1 = boxes[bp].childOffset + 1; + conduit::index_t child0 = boxes[bp].childOffset; + conduit::index_t child1 = boxes[bp].childOffset + 1; constructBox(child0, rangeA, A, level + 1, maxLevels); constructBox(child1, rangeB, B, level + 1, maxLevels); } @@ -527,14 +539,14 @@ void kdtree::constructBox(int bp, } //--------------------------------------------------------------------------- -template -int kdtree::longest(const kdtree::BoxType &b) const +template +int kdtree::longest(const kdtree::BoxType &b) const { int rv = 0; - T len = b[0][1] - b[0][0]; + CoordinateType len = b[0][1] - b[0][0]; for(int i = 1; i < dims(); i++) { - T newlen = b[i][1] - b[i][0]; + CoordinateType newlen = b[i][1] - b[i][0]; if(newlen > len) rv = i; } @@ -542,9 +554,9 @@ int kdtree::longest(const kdtree::BoxT } //--------------------------------------------------------------------------- -template -void kdtree::sortIndexRange( - const kdtree::RangeType &range, int dimension) +template +void kdtree::sortIndexRange( + const kdtree::RangeType &range, int dimension) { if(range.size > 1) { @@ -552,7 +564,7 @@ void kdtree::sortIndexRange( const Indexable &data = coords[dimension]; std::sort(index.begin() + range.offset, index.begin() + range.offset + range.size, - [&data](int idx1, int idx2) + [&data](conduit::index_t idx1, conduit::index_t idx2) { return data[idx1] < data[idx2]; }); @@ -560,9 +572,9 @@ void kdtree::sortIndexRange( } //--------------------------------------------------------------------------- -template -bool kdtree::pointInBox( - const kdtree::PointType &pt, const BoxType &b) const +template +bool kdtree::pointInBox( + const kdtree::PointType &pt, const BoxType &b) const { bool retval = true; for(int i = 0; i < dims(); i++) diff --git a/src/libs/blueprint/conduit_blueprint_mesh_partition.cpp b/src/libs/blueprint/conduit_blueprint_mesh_partition.cpp index 509aaf5a4..175ed5b66 100644 --- a/src/libs/blueprint/conduit_blueprint_mesh_partition.cpp +++ b/src/libs/blueprint/conduit_blueprint_mesh_partition.cpp @@ -812,6 +812,17 @@ SelectionExplicit::determine_is_whole(const conduit::Node &n_mesh) const for(index_t i = 0; i < n; i++) unique.insert(indices[i]); is_whole = static_cast(unique.size()) == num_elem_in_mesh; + + // If the mesh is whole, check that the indices are all in ascending + // order. If not, we're reordering and we should not consider the + // chunk available for passing through whole. + if(n > 1) + { + for(index_t i = 1; i < n && is_whole; i++) + { + is_whole &= (indices[i - 1] < indices[i]); + } + } } } catch(conduit::Error &) @@ -2768,11 +2779,11 @@ Partitioner::copy_field(const conduit::Node &n_field, for(index_t i = 0; i < n_values.number_of_children(); i++) { const conduit::Node &n_vals = n_values[i]; - slice_array(n_vals, ids, new_values[n_vals.name()]); + conduit::blueprint::mesh::utils::slice_array(n_vals, ids, new_values[n_vals.name()]); } } else - slice_array(n_values, ids, new_values); + conduit::blueprint::mesh::utils::slice_array(n_values, ids, new_values); } else { @@ -2785,130 +2796,11 @@ Partitioner::copy_field(const conduit::Node &n_field, for(index_t i = 0; i < n.number_of_children(); i++) { const conduit::Node &n_vals = n[i]; - slice_array(n_vals, ids, new_values[n_vals.name()]); + conduit::blueprint::mesh::utils::slice_array(n_vals, ids, new_values[n_vals.name()]); } } else - slice_array(n, ids, new_values); - } -} - -//--------------------------------------------------------------------------- -// @brief Slice the n_src array using the indices stored in ids. We use the -// array classes for their [] operators that deal with interleaved -// and non-interleaved arrays. -template -inline void -typed_slice_array(const T &src, const std::vector &ids, T &dest) -{ - size_t n = ids.size(); - for(size_t i = 0; i < n; i++) - dest[i] = src[ids[i]]; -} - -//--------------------------------------------------------------------------- -// @note Should this be part of conduit::Node or DataArray somehow. The number -// of times I've had to slice an array... -void -Partitioner::slice_array(const conduit::Node &n_src_values, - const std::vector &ids, Node &n_dest_values) const -{ - // Copy the DataType of the input conduit::Node but override the number of elements - // before copying it in so assigning to n_dest_values triggers a memory - // allocation. - auto dt = n_src_values.dtype(); - n_dest_values = DataType(n_src_values.dtype().id(), ids.size()); - - // Do the slice. - if(dt.is_int8()) - { - auto dest(n_dest_values.as_int8_array()); - typed_slice_array(n_src_values.as_int8_array(), ids, dest); - } - else if(dt.is_int16()) - { - auto dest(n_dest_values.as_int16_array()); - typed_slice_array(n_src_values.as_int16_array(), ids, dest); - } - else if(dt.is_int32()) - { - auto dest(n_dest_values.as_int32_array()); - typed_slice_array(n_src_values.as_int32_array(), ids, dest); - } - else if(dt.is_int64()) - { - auto dest(n_dest_values.as_int64_array()); - typed_slice_array(n_src_values.as_int64_array(), ids, dest); - } - else if(dt.is_uint8()) - { - auto dest(n_dest_values.as_uint8_array()); - typed_slice_array(n_src_values.as_uint8_array(), ids, dest); - } - else if(dt.is_uint16()) - { - auto dest(n_dest_values.as_uint16_array()); - typed_slice_array(n_src_values.as_uint16_array(), ids, dest); - } - else if(dt.is_uint32()) - { - auto dest(n_dest_values.as_uint32_array()); - typed_slice_array(n_src_values.as_uint32_array(), ids, dest); - } - else if(dt.is_uint64()) - { - auto dest(n_dest_values.as_uint64_array()); - typed_slice_array(n_src_values.as_uint64_array(), ids, dest); - } - else if(dt.is_char()) - { - auto dest(n_dest_values.as_char_array()); - typed_slice_array(n_src_values.as_char_array(), ids, dest); - } - else if(dt.is_short()) - { - auto dest(n_dest_values.as_short_array()); - typed_slice_array(n_src_values.as_short_array(), ids, dest); - } - else if(dt.is_int()) - { - auto dest(n_dest_values.as_int_array()); - typed_slice_array(n_src_values.as_int_array(), ids, dest); - } - else if(dt.is_long()) - { - auto dest(n_dest_values.as_long_array()); - typed_slice_array(n_src_values.as_long_array(), ids, dest); - } - else if(dt.is_unsigned_char()) - { - auto dest(n_dest_values.as_unsigned_char_array()); - typed_slice_array(n_src_values.as_unsigned_char_array(), ids, dest); - } - else if(dt.is_unsigned_short()) - { - auto dest(n_dest_values.as_unsigned_short_array()); - typed_slice_array(n_src_values.as_unsigned_short_array(), ids, dest); - } - else if(dt.is_unsigned_int()) - { - auto dest(n_dest_values.as_unsigned_int_array()); - typed_slice_array(n_src_values.as_unsigned_int_array(), ids, dest); - } - else if(dt.is_unsigned_long()) - { - auto dest(n_dest_values.as_unsigned_long_array()); - typed_slice_array(n_src_values.as_unsigned_long_array(), ids, dest); - } - else if(dt.is_float()) - { - auto dest(n_dest_values.as_float_array()); - typed_slice_array(n_src_values.as_float_array(), ids, dest); - } - else if(dt.is_double()) - { - auto dest(n_dest_values.as_double_array()); - typed_slice_array(n_src_values.as_double_array(), ids, dest); + conduit::blueprint::mesh::utils::slice_array(n, ids, new_values); } } @@ -3270,7 +3162,7 @@ Partitioner::create_new_rectilinear_coordset(const conduit::Node &n_coordset, indices.push_back(i); const conduit::Node &src = n_values[d]; - slice_array(src, indices, n_new_values[src.name()]); + conduit::blueprint::mesh::utils::slice_array(src, indices, n_new_values[src.name()]); } } @@ -3292,7 +3184,7 @@ Partitioner::create_new_explicit_coordset(const conduit::Node &n_coordset, { const conduit::Node &n_axis_values = n_values[axes[i]]; conduit::Node &n_new_axis_values = n_new_values[axes[i]]; - slice_array(n_axis_values, vertex_ids, n_new_axis_values); + conduit::blueprint::mesh::utils::slice_array(n_axis_values, vertex_ids, n_new_axis_values); } } else if(n_coordset["type"].as_string() == "rectilinear") @@ -3306,7 +3198,7 @@ Partitioner::create_new_explicit_coordset(const conduit::Node &n_coordset, { const conduit::Node &n_axis_values = n_values[axes[i]]; conduit::Node &n_new_axis_values = n_new_values[axes[i]]; - slice_array(n_axis_values, vertex_ids, n_new_axis_values); + conduit::blueprint::mesh::utils::slice_array(n_axis_values, vertex_ids, n_new_axis_values); } } else if(n_coordset["type"].as_string() == "explicit") @@ -3318,7 +3210,7 @@ Partitioner::create_new_explicit_coordset(const conduit::Node &n_coordset, { const conduit::Node &n_axis_values = n_values[axes[i]]; conduit::Node &n_new_axis_values = n_new_values[axes[i]]; - slice_array(n_axis_values, vertex_ids, n_new_axis_values); + conduit::blueprint::mesh::utils::slice_array(n_axis_values, vertex_ids, n_new_axis_values); } } } diff --git a/src/libs/blueprint/conduit_blueprint_mesh_partition.hpp b/src/libs/blueprint/conduit_blueprint_mesh_partition.hpp index fe9284ce8..059470a11 100644 --- a/src/libs/blueprint/conduit_blueprint_mesh_partition.hpp +++ b/src/libs/blueprint/conduit_blueprint_mesh_partition.hpp @@ -431,10 +431,6 @@ class CONDUIT_BLUEPRINT_API Partitioner const std::vector &ids, Node &n_output_fields) const; - void slice_array(const conduit::Node &n_src_values, - const std::vector &ids, - Node &n_dest_values) const; - void get_vertex_ids_for_element_ids(const conduit::Node &n_topo, const std::vector &element_ids, std::set &vertex_ids) const; diff --git a/src/libs/blueprint/conduit_blueprint_mesh_utils.cpp b/src/libs/blueprint/conduit_blueprint_mesh_utils.cpp index 0c9a17044..9ca1ca867 100644 --- a/src/libs/blueprint/conduit_blueprint_mesh_utils.cpp +++ b/src/libs/blueprint/conduit_blueprint_mesh_utils.cpp @@ -345,6 +345,227 @@ find_domain_id(const Node &node) return domain_id; } +//--------------------------------------------------------------------------- +// @brief Slice the n_src array using the indices stored in ids. We use the +// array classes for their [] operators that deal with interleaved +// and non-interleaved arrays. +template +inline void +typed_slice_array(const T &src, const std::vector &ids, T &dest) +{ + size_t n = ids.size(); + for(size_t i = 0; i < n; i++) + dest[i] = src[ids[i]]; +} + +//--------------------------------------------------------------------------- +/** + @brief Slice a node as an array using its native data type. + */ +template +void +slice_array_internal(const conduit::Node &n_src_values, + const std::vector &ids, + Node &n_dest_values) +{ + // Copy the DataType of the input conduit::Node but override the number of elements + // before copying it in so assigning to n_dest_values triggers a memory + // allocation. + auto dt = n_src_values.dtype(); + n_dest_values = DataType(n_src_values.dtype().id(), ids.size()); + + // Do the slice. + if(dt.is_int8()) + { + auto dest(n_dest_values.as_int8_array()); + typed_slice_array(n_src_values.as_int8_array(), ids, dest); + } + else if(dt.is_int16()) + { + auto dest(n_dest_values.as_int16_array()); + typed_slice_array(n_src_values.as_int16_array(), ids, dest); + } + else if(dt.is_int32()) + { + auto dest(n_dest_values.as_int32_array()); + typed_slice_array(n_src_values.as_int32_array(), ids, dest); + } + else if(dt.is_int64()) + { + auto dest(n_dest_values.as_int64_array()); + typed_slice_array(n_src_values.as_int64_array(), ids, dest); + } + else if(dt.is_uint8()) + { + auto dest(n_dest_values.as_uint8_array()); + typed_slice_array(n_src_values.as_uint8_array(), ids, dest); + } + else if(dt.is_uint16()) + { + auto dest(n_dest_values.as_uint16_array()); + typed_slice_array(n_src_values.as_uint16_array(), ids, dest); + } + else if(dt.is_uint32()) + { + auto dest(n_dest_values.as_uint32_array()); + typed_slice_array(n_src_values.as_uint32_array(), ids, dest); + } + else if(dt.is_uint64()) + { + auto dest(n_dest_values.as_uint64_array()); + typed_slice_array(n_src_values.as_uint64_array(), ids, dest); + } + else if(dt.is_float32()) + { + auto dest(n_dest_values.as_float32_array()); + typed_slice_array(n_src_values.as_float32_array(), ids, dest); + } + else if(dt.is_float64()) + { + auto dest(n_dest_values.as_float64_array()); + typed_slice_array(n_src_values.as_float64_array(), ids, dest); + } + else if(dt.is_char()) + { + auto dest(n_dest_values.as_char_array()); + typed_slice_array(n_src_values.as_char_array(), ids, dest); + } + else if(dt.is_short()) + { + auto dest(n_dest_values.as_short_array()); + typed_slice_array(n_src_values.as_short_array(), ids, dest); + } + else if(dt.is_int()) + { + auto dest(n_dest_values.as_int_array()); + typed_slice_array(n_src_values.as_int_array(), ids, dest); + } + else if(dt.is_long()) + { + auto dest(n_dest_values.as_long_array()); + typed_slice_array(n_src_values.as_long_array(), ids, dest); + } + else if(dt.is_unsigned_char()) + { + auto dest(n_dest_values.as_unsigned_char_array()); + typed_slice_array(n_src_values.as_unsigned_char_array(), ids, dest); + } + else if(dt.is_unsigned_short()) + { + auto dest(n_dest_values.as_unsigned_short_array()); + typed_slice_array(n_src_values.as_unsigned_short_array(), ids, dest); + } + else if(dt.is_unsigned_int()) + { + auto dest(n_dest_values.as_unsigned_int_array()); + typed_slice_array(n_src_values.as_unsigned_int_array(), ids, dest); + } + else if(dt.is_unsigned_long()) + { + auto dest(n_dest_values.as_unsigned_long_array()); + typed_slice_array(n_src_values.as_unsigned_long_array(), ids, dest); + } + else if(dt.is_float()) + { + auto dest(n_dest_values.as_float_array()); + typed_slice_array(n_src_values.as_float_array(), ids, dest); + } + else if(dt.is_double()) + { + auto dest(n_dest_values.as_double_array()); + typed_slice_array(n_src_values.as_double_array(), ids, dest); + } +} + +//--------------------------------------------------------------------------- +bool same_nodes(const conduit::Node &n1, const conduit::Node &n2) +{ + return (&n1 == &n2) || + (n1.contiguous_data_ptr() != nullptr && + n1.contiguous_data_ptr() == n2.contiguous_data_ptr()); +} + +//--------------------------------------------------------------------------- +void +slice_array(const conduit::Node &n_src_values, + const std::vector &ids, + conduit::Node &n_dest_values) +{ + // Check whether the src and dest nodes are the same. If so, we slice into + // a tmp node and move its contents. + if(same_nodes(n_src_values, n_dest_values)) + { + conduit::Node tmp; + slice_array_internal(n_src_values, ids, tmp); + n_dest_values.move(tmp); + } + else + { + slice_array_internal(n_src_values, ids, n_dest_values); + } +} + +//--------------------------------------------------------------------------- +void +slice_array(const conduit::Node &n_src_values, + const std::vector &ids, + conduit::Node &n_dest_values) +{ + // Check whether the src and dest nodes are the same. If so, we slice into + // a tmp node and move its contents. + if(same_nodes(n_src_values, n_dest_values)) + { + conduit::Node tmp; + slice_array_internal(n_src_values, ids, tmp); + n_dest_values.move(tmp); + } + else + { + slice_array_internal(n_src_values, ids, n_dest_values); + } +} + +//--------------------------------------------------------------------------- +template +void +slice_field_internal(const conduit::Node &n_src_values, + const std::vector &ids, + conduit::Node &n_dest_values) +{ + if(n_src_values.number_of_children() > 0) + { + // Reorder an mcarray + for(conduit::index_t ci = 0; ci < n_src_values.number_of_children(); ci++) + { + const conduit::Node &comp = n_src_values[ci]; + slice_array(comp, ids, n_dest_values[comp.name()]); + } + } + else + { + slice_array(n_src_values, ids, n_dest_values); + } +} + +//--------------------------------------------------------------------------- +void +slice_field(const conduit::Node &n_src_values, + const std::vector &ids, + conduit::Node &n_dest_values) +{ + slice_field_internal(n_src_values, ids, n_dest_values); +} + +//--------------------------------------------------------------------------- +void +slice_field(const conduit::Node &n_src_values, + const std::vector &ids, + conduit::Node &n_dest_values) +{ + slice_field_internal(n_src_values, ids, n_dest_values); +} + + //----------------------------------------------------------------------------- // -- begin conduit::blueprint::mesh::utils::connectivity -- //----------------------------------------------------------------------------- @@ -1158,6 +1379,32 @@ coordset::extents(const Node &n) return cset_extents; } +//----------------------------------------------------------------------------- +std::tuple +coordset::supports_pointer_access(const conduit::Node &coordset) +{ + bool suitable = false; + conduit::DataType dt; + if(coordset.has_child("values")) + { + suitable = true; + const conduit::Node &values = coordset.fetch_existing("values"); + for(conduit::index_t i = 0; i < values.number_of_children(); i++) + { + if(i == 0) + { + suitable &= values[i].dtype().is_compact(); + dt = values[i].dtype(); + } + else + { + suitable &= (values[i].dtype().is_compact() && dt.id() == values[i].dtype().id()); + } + } + } + return std::make_tuple(suitable, dt); +} + //----------------------------------------------------------------------------- // -- begin conduit::blueprint::mesh::utils::coordset::uniform -- //----------------------------------------------------------------------------- @@ -1388,345 +1635,94 @@ topology::reindex_coords(const Node& topo, } //----------------------------------------------------------------------------- -void -topology::unstructured::generate_offsets_inline(Node &topo) +template +static std::vector +spatial_ordering_impl(std::vector &coords, ElementType /*notused*/, conduit::index_t npts) { - // check for polyhedral case - if(topo.has_child("subelements")) + // Sort the coordinates spatially + std::vector reorder; + if(coords.size() == 2) { - // if ele or subelee offsets are missing or empty we want to generate - if( (!topo["elements"].has_child("offsets") || - topo["elements/offsets"].dtype().is_empty()) || - (!topo["subelements"].has_child("offsets") || - topo["subelements/offsets"].dtype().is_empty()) - ) - { - blueprint::mesh::utils::topology::unstructured::generate_offsets(topo, - topo["elements/offsets"], - topo["subelements/offsets"]); - } - + conduit::blueprint::mesh::utils::kdtree spatial_sort; + spatial_sort.initialize(&coords[0], npts); + reorder = std::move(spatial_sort.getIndices()); } - else + else if(coords.size() == 3) { - // if ele offsets is missing or empty we want to generate - if( !topo["elements"].has_child("offsets") || - topo["elements/offsets"].dtype().is_empty()) - { - blueprint::mesh::utils::topology::unstructured::generate_offsets(topo, - topo["elements/offsets"]); - } + conduit::blueprint::mesh::utils::kdtree spatial_sort; + spatial_sort.initialize(&coords[0], npts); + reorder = std::move(spatial_sort.getIndices()); } -} - - -//----------------------------------------------------------------------------- -void -topology::unstructured::generate_offsets(const Node &topo, - Node &ele_offsets) -{ - Node subele_offsets; - generate_offsets(topo,ele_offsets,subele_offsets); + return reorder; } //----------------------------------------------------------------------------- -void -topology::unstructured::generate_offsets(const Node &topo, - Node &dest_ele_offsets, - Node &dest_subele_offsets) +std::vector +topology::spatial_ordering(const conduit::Node &topo) { - const ShapeType topo_shape(topo); - const DataType int_dtype = find_widest_dtype(topo, DEFAULT_INT_DTYPES); - std::string key("elements/connectivity"), stream_key("elements/stream"); - - if(!topo.has_path(key)) - key = stream_key; - const Node &topo_conn = topo[key]; - - const DataType topo_dtype(topo_conn.dtype().id(), 1, 0, 0, - topo_conn.dtype().element_bytes(), topo_conn.dtype().endianness()); - - bool elem_offsets_exist = topo["elements"].has_child("offsets") && - !topo["elements/offsets"].dtype().is_empty(); - bool subelem_offsets_exist = false; + // Make a new centroid topo and coordset. The coordset will contain the + // element centers. This ought to be an explicit coordset. + Node topo_dest, coords_dest, s2dmap, d2smap; + mesh::topology::unstructured::generate_centroids(topo, + topo_dest, + coords_dest, + s2dmap, + d2smap); - // if these have already been generate, use set external to copy out results - if(topo_shape.type == "polyhedral") + std::vector reorder; + conduit::Node &values = coords_dest.fetch_existing("values"); + const auto pa = coordset::supports_pointer_access(coords_dest); + const auto suitable = std::get<0>(pa); + conduit::index_t npts{}; + if(suitable) { - subelem_offsets_exist = topo["subelements"].has_child("offsets") && - !topo["subelements/offsets"].dtype().is_empty(); - if(elem_offsets_exist && subelem_offsets_exist) + // The coordinates will be accessed using pointers. + const auto &dt = std::get<1>(pa); + if(dt.is_double()) { - // they are both already here, set external and return - if(&dest_ele_offsets != &topo["elements/offsets"]) + std::vector coords; + double elem{}; + for(conduit::index_t i = 0; i < values.number_of_children(); i++) { - dest_ele_offsets.set_external(topo["elements/offsets"]); + npts = (i == 0) ? (values[i].as_double_array().number_of_elements()) : npts; + coords.push_back(values[i].as_double_ptr()); } - if(&dest_subele_offsets != &topo["subelements/offsets"]) + reorder = spatial_ordering_impl(coords, elem, npts); + } + else if(dt.is_float()) + { + std::vector coords; + float elem{}; + for(conduit::index_t i = 0; i < values.number_of_children(); i++) { - dest_subele_offsets.set_external(topo["subelements/offsets"]); + npts = (i == 0) ? (values[i].as_float_array().number_of_elements()) : npts; + coords.push_back(values[i].as_float_ptr()); } - // we are done - return; + reorder = spatial_ordering_impl(coords, elem, npts); } } - else // non polyhedral + if(reorder.empty()) { - if(elem_offsets_exist) + // Use a double accessor to access coordinates. + std::vector coords; + double elem{}; + for(conduit::index_t i = 0; i < values.number_of_children(); i++) { - // they are already here, set external and return - if(&dest_ele_offsets != &topo["elements/offsets"]) - { - dest_ele_offsets.set_external(topo["elements/offsets"]); - } - return; + npts = (i == 0) ? (values[i].as_double_accessor().number_of_elements()) : npts; + coords.push_back(values[i].as_double_accessor()); } + reorder = spatial_ordering_impl(coords, elem, npts); } - // Selectively reset now that early returns have happened. We do the - // checks for the polyhedral case since some offsets might exist. - if(!elem_offsets_exist) - dest_ele_offsets.reset(); - if(!subelem_offsets_exist) - dest_subele_offsets.reset(); + return reorder; +} - /// - /// Generate Cases - /// - if(topo.has_path(stream_key)) - { - /// - /// TODO STREAM TOPOS ARE DEPRECATED - /// - // Mixed element types - std::map stream_id_npts; - const conduit::Node &n_element_types = topo["elements/element_types"]; - for(index_t i = 0; i < n_element_types.number_of_children(); i++) - { - const Node &n_et = n_element_types[i]; - auto stream_id = n_et["stream_id"].to_int(); - std::string shape(n_et["shape"].as_string()); - for(size_t j = 0; j < TOPO_SHAPES.size(); j++) - { - if(shape == TOPO_SHAPES[j]) - { - stream_id_npts[stream_id] = TOPO_SHAPE_EMBED_COUNTS[j]; - break; - } - } - } +//--------------------------------------------------------------------------- +topology::TopologyBuilder::TopologyBuilder(const conduit::Node &_topo) : topo(_topo), + old_to_new(), topo_conn(), topo_sizes() +{ - const Node &n_stream_ids = topo["elements/element_index/stream_ids"]; - std::vector offsets; - if(topo.has_path("elements/element_index/element_counts")) - { - const Node &n_element_counts = topo["elements/element_index/element_counts"]; - - index_t offset = 0, elemid = 0; - for(index_t j = 0; j < n_stream_ids.dtype().number_of_elements(); j++) - { - // Get the j'th elements from n_stream_ids, n_element_counts - const Node n_elem_ct_j(int_dtype, - const_cast(n_element_counts.element_ptr(j)), true); - const Node n_stream_ids_j(int_dtype, - const_cast(n_stream_ids.element_ptr(j)), true); - auto ec = static_cast(n_elem_ct_j.to_int64()); - auto sid = static_cast(n_stream_ids_j.to_int64()); - auto npts = stream_id_npts[sid]; - for(index_t i = 0; i < ec; i++) - { - offsets.push_back(offset); - offset += npts; - elemid++; - } - } - } - else if(topo.has_path("elements/element_index/offsets")) - { - const Node &n_stream = topo["elements/stream"]; - const Node &n_element_offsets = topo["elements/element_index/offsets"]; - index_t offset = 0, elemid = 0; - for(index_t j = 0; j < n_stream_ids.dtype().number_of_elements(); j++) - { - // Get the j'th elements from n_stream_ids, n_element_offsets - const Node n_stream_ids_j(int_dtype, - const_cast(n_stream_ids.element_ptr(j)), true); - const Node n_element_offsets_j(int_dtype, - const_cast(n_element_offsets.element_ptr(j)), true); - offset = n_element_offsets.to_index_t(); - index_t next_offset = offset; - if(j == n_stream_ids.dtype().number_of_elements() - 1) - { - next_offset = n_stream.dtype().number_of_elements(); - } - else - { - const Node n_element_offsets_j1(int_dtype, - const_cast(n_element_offsets.element_ptr(j)), true); - next_offset = n_element_offsets_j1.to_index_t(); - } - const auto sid = static_cast(n_stream_ids_j.to_int64()); - const auto npts = stream_id_npts[sid]; - while(offset < next_offset) { - offsets.push_back(offset); - offset += npts; - elemid++; - } - } - } - else - { - CONDUIT_ERROR("Stream based mixed topology has no element_counts or offsets.") - } - - Node off_node; - off_node.set_external(offsets); - off_node.to_data_type(int_dtype.id(), dest_ele_offsets); - } - else if(!topo_shape.is_poly()) - { - // Single element type - const index_t num_topo_shapes = - topo_conn.dtype().number_of_elements() / topo_shape.indices; - - Node shape_node(DataType::int64(num_topo_shapes)); - int64_array shape_array = shape_node.as_int64_array(); - for(index_t s = 0; s < num_topo_shapes; s++) - { - shape_array[s] = s * topo_shape.indices; - } - shape_node.to_data_type(int_dtype.id(), dest_ele_offsets); - } - else if(topo_shape.type == "polygonal") - { - const Node &topo_size = topo["elements/sizes"]; - int64_accessor topo_sizes = topo_size.as_int64_accessor(); - std::vector shape_array; - index_t i = 0; - index_t s = 0; - while(i < topo_size.dtype().number_of_elements()) - { - shape_array.push_back(s); - s += topo_sizes[i]; - i++; - } - - Node shape_node; - shape_node.set_external(shape_array); - shape_node.to_data_type(int_dtype.id(), dest_ele_offsets); - } - else if(topo_shape.type == "polyhedral") - { - // Construct any offsets that do not exist. - if(!elem_offsets_exist) - { - index_t_accessor topo_elem_size = topo["elements/sizes"].value(); - index_t es_count = topo_elem_size.number_of_elements(); - - dest_ele_offsets.set(DataType::index_t(es_count)); - index_t_array shape_array = dest_ele_offsets.value(); - - index_t es = 0; - for (index_t ei = 0; ei < es_count; ++ei) - { - shape_array[ei] = es; - es += topo_elem_size[ei]; - } - } - if(!subelem_offsets_exist) - { - index_t_accessor topo_subelem_size = topo["subelements/sizes"].value(); - index_t ses_count = topo_subelem_size.number_of_elements(); - - dest_subele_offsets.set(DataType::index_t(ses_count)); - index_t_array subshape_array = dest_subele_offsets.value(); - - index_t ses = 0; - for (index_t ei = 0; ei < ses_count; ++ei) - { - subshape_array[ei] = ses; - ses += topo_subelem_size[ei]; - } - } - } -} - - -//----------------------------------------------------------------------------- -std::vector -topology::unstructured::points(const Node &n, - const index_t ei) -{ - // NOTE(JRC): This is a workaround to ensure offsets are generated up-front - // if they don't exist and aren't regenerated for each subcall that needs them. - Node ntemp; - ntemp.set_external(n); - generate_offsets_inline(ntemp); - - const ShapeType topo_shape(ntemp); - - std::set pidxs; - if(!topo_shape.is_poly()) - { - index_t_accessor poff_vals = ntemp["elements/offsets"].value(); - const index_t eoff = poff_vals[ei]; - - index_t_accessor pidxs_vals = ntemp["elements/connectivity"].value(); - for(index_t pi = 0; pi < topo_shape.indices; pi++) - { - pidxs.insert(pidxs_vals[eoff + pi]); - } - } - else // if(topo_shape.is_poly()) - { - Node enode; - std::set eidxs; - if(topo_shape.is_polygonal()) - { - enode.set_external(ntemp["elements"]); - eidxs.insert(ei); - } - else // if(topo_shape.is_polyhedral()) - { - enode.set_external(ntemp["subelements"]); - - index_t_accessor eidxs_vals = ntemp["elements/connectivity"].value(); - o2mrelation::O2MIterator eiter(ntemp["elements"]); - eiter.to(ei, o2mrelation::ONE); - eiter.to_front(o2mrelation::MANY); - while(eiter.has_next(o2mrelation::MANY)) - { - eiter.next(o2mrelation::MANY); - const index_t tmp = eidxs_vals[eiter.index(o2mrelation::DATA)]; - eidxs.insert(tmp); - } - } - - index_t_accessor pidxs_vals = enode["connectivity"].value(); - o2mrelation::O2MIterator piter(enode); - for(const index_t eidx : eidxs) - { - piter.to(eidx, o2mrelation::ONE); - piter.to_front(o2mrelation::MANY); - while(piter.has_next(o2mrelation::MANY)) - { - piter.next(o2mrelation::MANY); - const index_t tmp = pidxs_vals[piter.index(o2mrelation::DATA)]; - pidxs.insert(tmp); - } - } - } - - return std::vector(pidxs.begin(), pidxs.end()); -} - - -//--------------------------------------------------------------------------- -topology::TopologyBuilder::TopologyBuilder(const conduit::Node &_topo) : topo(_topo), - old_to_new(), topo_conn(), topo_sizes() -{ - -} +} //--------------------------------------------------------------------------- topology::TopologyBuilder::TopologyBuilder(const conduit::Node *_topo) : topo(*_topo), @@ -1942,6 +1938,453 @@ topology::search(const conduit::Node &topo1, const conduit::Node &topo2) return exists; } +//----------------------------------------------------------------------------- +void +topology::unstructured::generate_offsets_inline(Node &topo) +{ + // check for polyhedral case + if(topo.has_child("subelements")) + { + // if ele or subelee offsets are missing or empty we want to generate + if( (!topo["elements"].has_child("offsets") || + topo["elements/offsets"].dtype().is_empty()) || + (!topo["subelements"].has_child("offsets") || + topo["subelements/offsets"].dtype().is_empty()) + ) + { + blueprint::mesh::utils::topology::unstructured::generate_offsets(topo, + topo["elements/offsets"], + topo["subelements/offsets"]); + } + + } + else + { + // if ele offsets is missing or empty we want to generate + if( !topo["elements"].has_child("offsets") || + topo["elements/offsets"].dtype().is_empty()) + { + blueprint::mesh::utils::topology::unstructured::generate_offsets(topo, + topo["elements/offsets"]); + } + } +} + + +//----------------------------------------------------------------------------- +void +topology::unstructured::generate_offsets(const Node &topo, + Node &ele_offsets) +{ + Node subele_offsets; + generate_offsets(topo,ele_offsets,subele_offsets); +} + +//----------------------------------------------------------------------------- +void +topology::unstructured::generate_offsets(const Node &topo, + Node &dest_ele_offsets, + Node &dest_subele_offsets) +{ + const ShapeType topo_shape(topo); + const DataType int_dtype = find_widest_dtype(topo, DEFAULT_INT_DTYPES); + std::string key("elements/connectivity"), stream_key("elements/stream"); + + if(!topo.has_path(key)) + key = stream_key; + const Node &topo_conn = topo[key]; + + const DataType topo_dtype(topo_conn.dtype().id(), 1, 0, 0, + topo_conn.dtype().element_bytes(), topo_conn.dtype().endianness()); + + bool elem_offsets_exist = topo["elements"].has_child("offsets") && + !topo["elements/offsets"].dtype().is_empty(); + bool subelem_offsets_exist = false; + + // if these have already been generate, use set external to copy out results + if(topo_shape.type == "polyhedral") + { + subelem_offsets_exist = topo["subelements"].has_child("offsets") && + !topo["subelements/offsets"].dtype().is_empty(); + if(elem_offsets_exist && subelem_offsets_exist) + { + // they are both already here, set external and return + if(&dest_ele_offsets != &topo["elements/offsets"]) + { + dest_ele_offsets.set_external(topo["elements/offsets"]); + } + if(&dest_subele_offsets != &topo["subelements/offsets"]) + { + dest_subele_offsets.set_external(topo["subelements/offsets"]); + } + // we are done + return; + } + } + else // non polyhedral + { + if(elem_offsets_exist) + { + // they are already here, set external and return + if(&dest_ele_offsets != &topo["elements/offsets"]) + { + dest_ele_offsets.set_external(topo["elements/offsets"]); + } + return; + } + } + + // Selectively reset now that early returns have happened. We do the + // checks for the polyhedral case since some offsets might exist. + if(!elem_offsets_exist) + dest_ele_offsets.reset(); + if(!subelem_offsets_exist) + dest_subele_offsets.reset(); + + /// + /// Generate Cases + /// + if(topo.has_path(stream_key)) + { + /// + /// TODO STREAM TOPOS ARE DEPRECATED + /// + // Mixed element types + std::map stream_id_npts; + const conduit::Node &n_element_types = topo["elements/element_types"]; + for(index_t i = 0; i < n_element_types.number_of_children(); i++) + { + const Node &n_et = n_element_types[i]; + auto stream_id = n_et["stream_id"].to_int(); + std::string shape(n_et["shape"].as_string()); + for(size_t j = 0; j < TOPO_SHAPES.size(); j++) + { + if(shape == TOPO_SHAPES[j]) + { + stream_id_npts[stream_id] = TOPO_SHAPE_EMBED_COUNTS[j]; + break; + } + } + } + + const Node &n_stream_ids = topo["elements/element_index/stream_ids"]; + std::vector offsets; + if(topo.has_path("elements/element_index/element_counts")) + { + const Node &n_element_counts = topo["elements/element_index/element_counts"]; + + index_t offset = 0, elemid = 0; + for(index_t j = 0; j < n_stream_ids.dtype().number_of_elements(); j++) + { + // Get the j'th elements from n_stream_ids, n_element_counts + const Node n_elem_ct_j(int_dtype, + const_cast(n_element_counts.element_ptr(j)), true); + const Node n_stream_ids_j(int_dtype, + const_cast(n_stream_ids.element_ptr(j)), true); + auto ec = static_cast(n_elem_ct_j.to_int64()); + auto sid = static_cast(n_stream_ids_j.to_int64()); + auto npts = stream_id_npts[sid]; + for(index_t i = 0; i < ec; i++) + { + offsets.push_back(offset); + offset += npts; + elemid++; + } + } + } + else if(topo.has_path("elements/element_index/offsets")) + { + const Node &n_stream = topo["elements/stream"]; + const Node &n_element_offsets = topo["elements/element_index/offsets"]; + index_t offset = 0, elemid = 0; + for(index_t j = 0; j < n_stream_ids.dtype().number_of_elements(); j++) + { + // Get the j'th elements from n_stream_ids, n_element_offsets + const Node n_stream_ids_j(int_dtype, + const_cast(n_stream_ids.element_ptr(j)), true); + const Node n_element_offsets_j(int_dtype, + const_cast(n_element_offsets.element_ptr(j)), true); + offset = n_element_offsets.to_index_t(); + index_t next_offset = offset; + if(j == n_stream_ids.dtype().number_of_elements() - 1) + { + next_offset = n_stream.dtype().number_of_elements(); + } + else + { + const Node n_element_offsets_j1(int_dtype, + const_cast(n_element_offsets.element_ptr(j)), true); + next_offset = n_element_offsets_j1.to_index_t(); + } + const auto sid = static_cast(n_stream_ids_j.to_int64()); + const auto npts = stream_id_npts[sid]; + while(offset < next_offset) { + offsets.push_back(offset); + offset += npts; + elemid++; + } + } + } + else + { + CONDUIT_ERROR("Stream based mixed topology has no element_counts or offsets.") + } + + Node off_node; + off_node.set_external(offsets); + off_node.to_data_type(int_dtype.id(), dest_ele_offsets); + } + else if(!topo_shape.is_poly()) + { + // Single element type + const index_t num_topo_shapes = + topo_conn.dtype().number_of_elements() / topo_shape.indices; + + Node shape_node(DataType::int64(num_topo_shapes)); + int64_array shape_array = shape_node.as_int64_array(); + for(index_t s = 0; s < num_topo_shapes; s++) + { + shape_array[s] = s * topo_shape.indices; + } + shape_node.to_data_type(int_dtype.id(), dest_ele_offsets); + } + else if(topo_shape.type == "polygonal") + { + const Node &topo_size = topo["elements/sizes"]; + int64_accessor topo_sizes = topo_size.as_int64_accessor(); + std::vector shape_array; + index_t i = 0; + index_t s = 0; + while(i < topo_size.dtype().number_of_elements()) + { + shape_array.push_back(s); + s += topo_sizes[i]; + i++; + } + + Node shape_node; + shape_node.set_external(shape_array); + shape_node.to_data_type(int_dtype.id(), dest_ele_offsets); + } + else if(topo_shape.type == "polyhedral") + { + // Construct any offsets that do not exist. + if(!elem_offsets_exist) + { + index_t_accessor topo_elem_size = topo["elements/sizes"].value(); + index_t es_count = topo_elem_size.number_of_elements(); + + dest_ele_offsets.set(DataType::index_t(es_count)); + index_t_array shape_array = dest_ele_offsets.value(); + + index_t es = 0; + for (index_t ei = 0; ei < es_count; ++ei) + { + shape_array[ei] = es; + es += topo_elem_size[ei]; + } + } + if(!subelem_offsets_exist) + { + index_t_accessor topo_subelem_size = topo["subelements/sizes"].value(); + index_t ses_count = topo_subelem_size.number_of_elements(); + + dest_subele_offsets.set(DataType::index_t(ses_count)); + index_t_array subshape_array = dest_subele_offsets.value(); + + index_t ses = 0; + for (index_t ei = 0; ei < ses_count; ++ei) + { + subshape_array[ei] = ses; + ses += topo_subelem_size[ei]; + } + } + } +} + + +//----------------------------------------------------------------------------- +std::vector +topology::unstructured::points(const Node &n, + const index_t ei) +{ + // NOTE(JRC): This is a workaround to ensure offsets are generated up-front + // if they don't exist and aren't regenerated for each subcall that needs them. + Node ntemp; + ntemp.set_external(n); + generate_offsets_inline(ntemp); + + const ShapeType topo_shape(ntemp); + + std::set pidxs; + if(!topo_shape.is_poly()) + { + index_t_accessor poff_vals = ntemp["elements/offsets"].value(); + const index_t eoff = poff_vals[ei]; + + index_t_accessor pidxs_vals = ntemp["elements/connectivity"].value(); + for(index_t pi = 0; pi < topo_shape.indices; pi++) + { + pidxs.insert(pidxs_vals[eoff + pi]); + } + } + else // if(topo_shape.is_poly()) + { + Node enode; + std::set eidxs; + if(topo_shape.is_polygonal()) + { + enode.set_external(ntemp["elements"]); + eidxs.insert(ei); + } + else // if(topo_shape.is_polyhedral()) + { + enode.set_external(ntemp["subelements"]); + + index_t_accessor eidxs_vals = ntemp["elements/connectivity"].value(); + o2mrelation::O2MIterator eiter(ntemp["elements"]); + eiter.to(ei, o2mrelation::ONE); + eiter.to_front(o2mrelation::MANY); + while(eiter.has_next(o2mrelation::MANY)) + { + eiter.next(o2mrelation::MANY); + const index_t tmp = eidxs_vals[eiter.index(o2mrelation::DATA)]; + eidxs.insert(tmp); + } + } + + index_t_accessor pidxs_vals = enode["connectivity"].value(); + o2mrelation::O2MIterator piter(enode); + for(const index_t eidx : eidxs) + { + piter.to(eidx, o2mrelation::ONE); + piter.to_front(o2mrelation::MANY); + while(piter.has_next(o2mrelation::MANY)) + { + piter.next(o2mrelation::MANY); + const index_t tmp = pidxs_vals[piter.index(o2mrelation::DATA)]; + pidxs.insert(tmp); + } + } + } + + return std::vector(pidxs.begin(), pidxs.end()); +} + +//----------------------------------------------------------------------------- +void +topology::unstructured::reorder(const conduit::Node &topo, + const conduit::Node &coordset, + const conduit::Node &fields, + const std::vector &reorder, + conduit::Node &dest_topo, + conduit::Node &dest_coordset, + conduit::Node &dest_fields, + std::vector &old2NewPoints) +{ + conduit::blueprint::mesh::utils::ShapeType shape(topo); + + // Handle unstructured meshes (but not polyhedral meshes yet) + if(topo.fetch_existing("type").as_string() == "unstructured" && !shape.is_polyhedral()) + { + // Input connectivity information. + const auto &n_conn = topo.fetch_existing("elements/connectivity"); + const auto &n_sizes = topo.fetch_existing("elements/sizes"); + const auto &n_offsets = topo.fetch_existing("elements/offsets"); + const auto conn = n_conn.as_index_t_accessor(); + const auto sizes = n_sizes.as_index_t_accessor(); + const auto offsets = n_offsets.as_index_t_accessor(); + + // Temp vectors to store reordered connectivity. We use temp vectors so + // we can convert to a matching datatype after we've constructed the data. + std::vector newconn, newoffsets, newsizes; + newconn.reserve(conn.number_of_elements()); + newsizes.reserve(sizes.number_of_elements()); + newoffsets.reserve(offsets.number_of_elements()); + + // Mapping information for the points. + constexpr conduit::index_t invalidNode = -1; + auto npts = conduit::blueprint::mesh::coordset::length(coordset); + // Fill in the old2New point mapping. It gets passed out of the function. + old2NewPoints.resize(npts); + for(conduit::index_t i = 0; i < npts; i++) + old2NewPoints[i] = invalidNode; + // ptReorder is used to reorder/slice vertex-associated data. We'll allow + // up to npts values but it might not be that large if we are selecting a + // subset of elements. + std::vector ptReorder; + ptReorder.reserve(npts); + conduit::index_t newPointIndex = 0; + + // We iterate over elements in the specified order. We iterate over the + // points in each element and renumber the points. + conduit::index_t newoffset = 0; + for(const auto cellIndex : reorder) + { + for(conduit::index_t i = 0; i < sizes[cellIndex]; i++) + { + auto id = conn[offsets[cellIndex] + i]; + // if the old point has not been seen, renumber it. + if(old2NewPoints[id] == invalidNode) + { + ptReorder.push_back(id); + old2NewPoints[id] = newPointIndex++; + } + newconn.push_back(old2NewPoints[id]); + } + newsizes.push_back(sizes[cellIndex]); + newoffsets.push_back(newoffset); + newoffset += sizes[cellIndex]; + } + + // Store the new connectivity. + dest_topo["type"] = topo["type"]; + dest_topo["coordset"] = dest_coordset.name(); + dest_topo["elements/shape"] = topo["elements/shape"]; + conduit::Node tmp; + tmp.set_external(newconn.data(), newconn.size()); + tmp.to_data_type(n_conn.dtype().id(), dest_topo["elements/connectivity"]); + tmp.set_external(newsizes.data(), newsizes.size()); + tmp.to_data_type(n_sizes.dtype().id(), dest_topo["elements/sizes"]); + tmp.set_external(newoffsets.data(), newoffsets.size()); + tmp.to_data_type(n_offsets.dtype().id(), dest_topo["elements/offsets"]); + + // Reorder the coordset now, making it explicit if needed. + dest_coordset["type"] = "explicit"; + conduit::Node coordset_explicit; + conduit::blueprint::mesh::coordset::to_explicit(coordset, coordset_explicit); + conduit::blueprint::mesh::utils::slice_field(coordset_explicit["values"], ptReorder, dest_coordset["values"]); + + // Reorder fields that match this topo. + std::vector fieldNames; + for(conduit::index_t fi = 0; fi < fields.number_of_children(); fi++) + { + const conduit::Node &src = fields[fi]; + if(src["topology"].as_string() == topo.name()) + { + fieldNames.push_back(src.name()); + } + } + for(const auto &fieldName : fieldNames) + { + const conduit::Node &src = fields.fetch_existing(fieldName); + conduit::Node &dest = dest_fields[fieldName]; + dest["association"] = src["association"]; + dest["topology"] = dest_topo.name(); + if(dest["association"].as_string() == "element") + { + conduit::blueprint::mesh::utils::slice_field(src["values"], reorder, dest["values"]); + } + else + { + conduit::blueprint::mesh::utils::slice_field(src["values"], ptReorder, dest["values"]); + } + } + + // TODO: renumber adjsets. + } +} + //----------------------------------------------------------------------------- // -- end conduit::blueprint::mesh::utils::topology -- //----------------------------------------------------------------------------- @@ -2491,7 +2934,7 @@ PointQuery::acceleratedSearch(int ndims, float64 searchPt[3] = {static_cast(input_ptr[i * 3 + 0]), static_cast(input_ptr[i * 3 + 1]), static_cast(input_ptr[i * 3 + 2])}; - int found = search.findPoint(searchPt); + auto found = static_cast(search.findPoint(searchPt)); result_ptr[i] = (found != search.NotFound) ? found : NotFound; }); handled = true; @@ -2514,7 +2957,7 @@ PointQuery::acceleratedSearch(int ndims, float32 searchPt[3] = {static_cast(input_ptr[i * 3 + 0]), static_cast(input_ptr[i * 3 + 1]), static_cast(input_ptr[i * 3 + 2])}; - int found = search.findPoint(searchPt); + auto found = static_cast(search.findPoint(searchPt)); result_ptr[i] = (found != search.NotFound) ? found : NotFound; }); handled = true; @@ -2536,7 +2979,7 @@ PointQuery::acceleratedSearch(int ndims, { float64 searchPt[2] = {static_cast(input_ptr[i * 3 + 0]), static_cast(input_ptr[i * 3 + 1])}; - int found = search.findPoint(searchPt); + auto found = static_cast(search.findPoint(searchPt)); result_ptr[i] = (found != search.NotFound) ? found : NotFound; }); handled = true; @@ -2557,7 +3000,7 @@ PointQuery::acceleratedSearch(int ndims, { float32 searchPt[2] = {static_cast(input_ptr[i * 3 + 0]), static_cast(input_ptr[i * 3 + 1])}; - int found = search.findPoint(searchPt); + auto found = static_cast(search.findPoint(searchPt)); result_ptr[i] = (found != search.NotFound) ? found : NotFound; }); handled = true; diff --git a/src/libs/blueprint/conduit_blueprint_mesh_utils.hpp b/src/libs/blueprint/conduit_blueprint_mesh_utils.hpp index 8eac7b37a..081ee845e 100644 --- a/src/libs/blueprint/conduit_blueprint_mesh_utils.hpp +++ b/src/libs/blueprint/conduit_blueprint_mesh_utils.hpp @@ -17,6 +17,7 @@ #include #include #include +#include #include #include @@ -216,6 +217,41 @@ CONDUIT_BLUEPRINT_API const Node * find_reference_node(const Node &node, const s //----------------------------------------------------------------------------- index_t CONDUIT_BLUEPRINT_API find_domain_id(const Node &node); +//----------------------------------------------------------------------------- +/** + @brief Slice a node containing array data and copy data for the supplied ids + to a new node. + + @param n_src_values The node containing the source values. + @param ids The ids that will be extracted from the source values. + @param dest The new node that will contain the sliced data. + */ +void CONDUIT_BLUEPRINT_API slice_array(const conduit::Node &n_src_values, + const std::vector &ids, + Node &n_dest_values); + +/// Same as above. +void CONDUIT_BLUEPRINT_API slice_array(const conduit::Node &n_src_values, + const std::vector &ids, + Node &n_dest_values); + +//----------------------------------------------------------------------------- +/** + @brief Slice a values node for a field where values may be an mcarray. The new + node will contain the data for the supplied indices. + + @param n_src_values The node containing the source values. + @param ids The ids that will be extracted from the source values. + @param dest The new node that will contain the sliced data. + */ +void CONDUIT_BLUEPRINT_API slice_field(const conduit::Node &n_src_values, + const std::vector &ids, + conduit::Node &dest); + +/// Same as above. +void CONDUIT_BLUEPRINT_API slice_field(const conduit::Node &n_src_values, + const std::vector &ids, + conduit::Node &n_dest_values); //----------------------------------------------------------------------------- // -- begin conduit::blueprint::mesh::utils::connectivity -- @@ -312,6 +348,18 @@ namespace coordset */ std::vector CONDUIT_BLUEPRINT_API extents(const Node &n); + //----------------------------------------------------------------------------- + /** + @brief Check whether all component types match and are compact, making it + suitable for pointer access. + + @param cset The coordset we're checking. + + @return A tuple containing 1) whether the components are suitable for pointer + access and 2) the data type. + */ + std::tuple CONDUIT_BLUEPRINT_API supports_pointer_access(const conduit::Node &coordset); + namespace uniform { /** @@ -369,6 +417,17 @@ namespace topology const conduit::Node& old_gvids, const conduit::Node& new_gvids, conduit::Node& out_topo); + //------------------------------------------------------------------------- + /** + * @brief Applies a spatial sorting algorithm (based on a kdtree) to the + * topology's centroids and returns a vector containing the sorted + * order. + * + * @param topo The topology whose elements are being sorted. + * + * @return A vector containing the new element order. + */ + std::vector CONDUIT_BLUEPRINT_API spatial_ordering(const conduit::Node &topo); //------------------------------------------------------------------------- /** @@ -461,6 +520,40 @@ namespace topology //------------------------------------------------------------------------- std::vector CONDUIT_BLUEPRINT_API points(const Node &topo, const index_t i); + + //------------------------------------------------------------------------- + /** + * @brief Reorder the topology's elements and nodes, according to the new + * element \order vector. + * + * @param topo A node containing the topo to be reordered. + * @param coordset A node containing the coordset to be reordered. + * @param fields A node containing the fields to be reordered. Only fields + * that match the topology node name will be modified. + * @param order A vector containing element indices in their new order. + * @param dest_topo A node that will contain reordered topo. It can + * be the same node as topo. + * @param dest_coordset A node that will contain reordered coordset. It can + * be the same node as coordset. + * @param dest_fields A node that will contain the reordered fields. It can + * be the same node as fields. + * @param[out] old2NewPoints A vector that contains new point indices for + * old point indices, which can be used for mapping + * data from the original order to the new order. + * + * @note This reorder function is similar to calling partition with an explicit + * index selection if the indices are provided in a new order. This + * function will also reorder the nodes in their order of use by elements + * in the new order and partition does not currently do that. + */ + void CONDUIT_BLUEPRINT_API reorder(const conduit::Node &topo, + const conduit::Node &coordset, + const conduit::Node &fields, + const std::vector &order, + conduit::Node &dest_topo, + conduit::Node &dest_coordset, + conduit::Node &dest_fields, + std::vector &old2NewPoints); } //------------------------------------------------------------------------- // -- end conduit::blueprint::mesh::utils::topology::unstructured -- diff --git a/src/libs/blueprint/conduit_blueprint_mpi_mesh_utils.cpp b/src/libs/blueprint/conduit_blueprint_mpi_mesh_utils.cpp index 9b853deb3..fc95f0be5 100644 --- a/src/libs/blueprint/conduit_blueprint_mpi_mesh_utils.cpp +++ b/src/libs/blueprint/conduit_blueprint_mpi_mesh_utils.cpp @@ -635,6 +635,114 @@ adjset::validate(const Node &doms, return retval; } +//----------------------------------------------------------------------------- +bool +adjset::compare_pointwise(conduit::Node &mesh, const std::string &adjsetName, MPI_Comm comm) +{ + namespace bputils = conduit::blueprint::mesh::utils; + std::vector domains = conduit::blueprint::mesh::domains(mesh); + + // Determine total number of domains. + conduit::Node nd_local, nd_total; + nd_local = static_cast(domains.size()); + relay::mpi::sum_all_reduce(nd_local, nd_total, comm); + int maxDomains = nd_total.to_int(); + const int par_rank = relay::mpi::rank(comm); + + // Figure out which MPI ranks own each domain. + conduit::Node n_domain2rank; + { + conduit::Node n_local(DataType::int32(maxDomains)); + int *iptr = n_local.as_int_ptr(); + memset(iptr, 0, sizeof(int) * maxDomains); + for(size_t i = 0; i < domains.size(); i++) + { + auto domainId = bputils::find_domain_id(*domains[i]); + iptr[domainId] = par_rank; + } + relay::mpi::sum_all_reduce(n_local, n_domain2rank, comm); + } + const int *domain2rank = n_domain2rank.as_int_ptr(); + + // Iterate over each of the possible adjset relationships. Not all of these + // will have adjset groups. + const int tag = 12211221; + for(int d0 = 0; d0 < maxDomains; d0++) + { + for(int d1 = d0 + 1; d1 < maxDomains; d1++) + { + // make the adjset group name. + std::stringstream ss; + ss << "group_" << d0 << "_" << d1; + std::string groupName(ss.str()); + + // There are up to 2 local meshes and their corresponding remote mesh. + relay::mpi::communicate_using_schema C(comm); + conduit::Node localMesh[2], remoteMesh[2]; + int mi = 0; + for(auto dom_ptr : domains) + { + Node &domain = *dom_ptr; + auto domainId = bputils::find_domain_id(domain); + + // If the domain has the adjset, make a point mesh of its points + // that we can send to the neighbor. + std::string key("adjsets/" + adjsetName + "/groups/" + groupName + "/values"); + if(domain.has_path(key)) + { + // Get the topology that the adjset wants. + std::string tkey("adjsets/" + adjsetName + "/topology"); + std::string topoName = domain.fetch_existing(tkey).as_string(); + const Node &topo = domain.fetch_existing("topologies/" + topoName); + + // Get the group values and add them as points to the topo builder + // so we pull out a point mesh. + std::string key("adjsets/" + adjsetName + "/groups/" + groupName + "/values"); + const Node &n_values = domain.fetch_existing(key); + const auto values = n_values.as_index_t_accessor(); + bputils::topology::TopologyBuilder B(topo); + for(index_t i = 0; i < values.number_of_elements(); i++) + { + index_t ptid = values[i]; + B.add(&ptid, 1); + } + + // Make the local point mesh. + B.execute(localMesh[mi], "vertex"); + + // Get the neighbor for this group + std::string nkey("adjsets/" + adjsetName + "/groups/" + groupName + "/neighbors"); + int neighbor = domain.fetch_existing(nkey).to_int(); + + // Send this local mesh to the neighbor. + C.add_isend(localMesh[mi], domain2rank[neighbor], tag + mi); + // That neighbor will have to send us a mesh too. + C.add_irecv(remoteMesh[mi], domain2rank[neighbor], tag + mi); + + mi++; + } + } + + // Perform the exchange. + C.execute(); + + // Make sure the nodes are not different. + Node different, reducedDiff; + different = 0; + for(int i = 0; i < mi; i++) + { + Node info; + bool d = localMesh[i].diff(remoteMesh[i], info, 1.e-8); + different = different.to_int() + (d ? 1 : 0); + } + relay::mpi::sum_all_reduce(different, reducedDiff, comm); + if(reducedDiff.to_int() > 0) + return false; + } + } + return true; +} + //----------------------------------------------------------------------------- // -- end conduit::blueprint::mpi::mesh::utils::adjset -- //----------------------------------------------------------------------------- diff --git a/src/libs/blueprint/conduit_blueprint_mpi_mesh_utils.hpp b/src/libs/blueprint/conduit_blueprint_mpi_mesh_utils.hpp index 0735cb3dc..d0ca824c0 100644 --- a/src/libs/blueprint/conduit_blueprint_mpi_mesh_utils.hpp +++ b/src/libs/blueprint/conduit_blueprint_mpi_mesh_utils.hpp @@ -166,6 +166,22 @@ namespace adjset Node &info, MPI_Comm comm); + /** + @brief Traverse the adjset groups and make sure that the points are the same + on both sides of the interface. This is more restrictive than just + checking whether they exist on the other domain. Now, they have to be + the same point. + + @param mesh A node that contains one or more mesh domains. + @param adjsetName The name of the adjset to check. This must be a pairwise adjset. + @param comm The MPI communicator to use. + + @return True if the adjset are the same pointwise across each interface; + False otherwise. + */ + bool CONDUIT_BLUEPRINT_API compare_pointwise(conduit::Node &mesh, + const std::string &adjsetName, + MPI_Comm comm); } //----------------------------------------------------------------------------- // -- end conduit::blueprint::mpi::mesh::utils::adjset -- diff --git a/src/tests/blueprint/CMakeLists.txt b/src/tests/blueprint/CMakeLists.txt index dc4288f1e..51a5553c5 100644 --- a/src/tests/blueprint/CMakeLists.txt +++ b/src/tests/blueprint/CMakeLists.txt @@ -85,7 +85,8 @@ set(BLUEPRINT_RELAY_PARMETIS_MPI_TESTS set(BLUEPRINT_RELAY_MPI_TESTS_RANKS_4 t_blueprint_mpi_mesh_flatten t_blueprint_mpi_mesh_partition - t_blueprint_mpi_mesh_distribute_4_ranks) + t_blueprint_mpi_mesh_distribute_4_ranks + t_blueprint_mpi_mesh_tiled) if(MPI_FOUND) message(STATUS "MPI enabled: Adding conduit_blueprint_mpi and conduit_relay_mpi unit tests") diff --git a/src/tests/blueprint/t_blueprint_mpi_mesh_tiled.cpp b/src/tests/blueprint/t_blueprint_mpi_mesh_tiled.cpp new file mode 100644 index 000000000..b1ef81d96 --- /dev/null +++ b/src/tests/blueprint/t_blueprint_mpi_mesh_tiled.cpp @@ -0,0 +1,296 @@ +// Copyright (c) Lawrence Livermore National Security, LLC and other Conduit +// Project developers. See top-level LICENSE AND COPYRIGHT files for dates and +// other details. No copyright assignment is required to contribute to Conduit. + +//----------------------------------------------------------------------------- +/// +/// file: t_blueprint_mpi_mesh_tiled.cpp +/// +//----------------------------------------------------------------------------- + +#include "conduit.hpp" +#include "conduit_blueprint.hpp" +#include "conduit_blueprint_mesh_examples.hpp" +#include "conduit_blueprint_mpi_mesh.hpp" +#include "conduit_blueprint_mpi_mesh_utils.hpp" +#include "conduit_relay.hpp" +#include "conduit_relay_mpi.hpp" +#include "conduit_relay_mpi_io_blueprint.hpp" +#include "conduit_log.hpp" + +#include "blueprint_test_helpers.hpp" + +#include +#include +#include +#include +#include "gtest/gtest.h" + +using namespace conduit; +using namespace conduit::utils; +using namespace generate; + +// Uncomment if we want to write the data files. +//#define CONDUIT_WRITE_TEST_DATA + +//--------------------------------------------------------------------------- +#ifdef CONDUIT_WRITE_TEST_DATA +/** + @brief Save the node to an HDF5 compatible with VisIt or the + conduit_adjset_validate tool. + */ +void save_mesh(const conduit::Node &root, const std::string &filebase) +{ + // NOTE: Enable this to write files for debugging. + const std::string protocol("hdf5"); + conduit::relay::mpi::io::blueprint::save_mesh(root, filebase, protocol, MPI_COMM_WORLD); +} +#endif + +//----------------------------------------------------------------------------- +/** + @brief Make domains for a tiled mesh. + + @param mesh The node that will contain the domains. + @param dims The number of zones in each domain. dims[2] == 0 for 2D meshes. + @param domains The number of domains to make. All values > 0. + @param reorder The reordering method, if any. + @param domainNumbering An optional vector that can reorder the domains. + */ +void +make_tiled(conduit::Node &mesh, const int dims[3], const int domains[3], + const std::string &reorder, const std::vector &domainNumbering) +{ + const int ndoms = domains[0] * domains[1] * domains[2]; + const int par_rank = relay::mpi::rank(MPI_COMM_WORLD); + const int par_size = relay::mpi::size(MPI_COMM_WORLD); + std::vector domains_per_rank(par_size, 0); + for(int di = 0; di < ndoms; di++) + domains_per_rank[di % par_size]++; + int offset = 0; + for(int i = 0; i < par_rank; i++) + offset += domains_per_rank[i]; + + // Make domains. + const double extents[] = {0., 1., 0., 1., 0., 1.}; + int domainIndex = 0; + for(int k = 0; k < domains[2]; k++) + for(int j = 0; j < domains[1]; j++) + for(int i = 0; i < domains[0]; i++, domainIndex++) + { + int domain[] = {i,j,k}; + + int domainId = domainIndex; + if(!domainNumbering.empty()) + domainId = domainNumbering[domainIndex]; + + // See if we need to make the domain on this rank. + if(domainId >= offset && domainId < offset + domains_per_rank[par_rank]) + { + // Determine the size and location of this domain in the whole. + double sideX = (extents[1] - extents[0]) / static_cast(domains[0]); + double sideY = (extents[3] - extents[2]) / static_cast(domains[1]); + double sideZ = (extents[5] - extents[4]) / static_cast(domains[2]); + double domainExt[] = {extents[0] + domain[0] * sideX, + extents[0] + (domain[0]+1) * sideX, + extents[2] + domain[1] * sideY, + extents[2] + (domain[1]+1) * sideY, + extents[4] + domain[2] * sideZ, + extents[4] + (domain[2]+1) * sideZ}; + conduit::Node opts; + opts["extents"].set(domainExt, 6); + opts["domain"].set(domain, 3); + opts["domains"].set(domains, 3); + opts["reorder"] = reorder; + + if(ndoms > 1) + { + char domainName[64]; + if(!domainNumbering.empty()) + sprintf(domainName, "domain_%07d", domainNumbering[domainIndex]); + else + sprintf(domainName, "domain_%07d", domainIndex); + conduit::Node &dom = mesh[domainName]; + conduit::blueprint::mesh::examples::tiled(dims[0], dims[1], dims[2], dom, opts); + } + else + { + conduit::blueprint::mesh::examples::tiled(dims[0], dims[1], dims[2], mesh, opts); + } + } + } +} + +//----------------------------------------------------------------------------- +template +void +foreach_permutation_impl(int level, std::vector &values, Func &&func) +{ + int n = static_cast(values.size()); + if(level < n) + { + int n = static_cast(values.size()); + for(int i = 0; i < n; i++) + { + const int *start = &values[0]; + const int *end = start + level; + if(std::find(start, end, i) == end) + { + values[level] = i; + foreach_permutation_impl(level + 1, values, func); + } + } + } + else + { + func(values); + } +} + +template +void +foreach_permutation(int maxval, Func &&func) +{ + std::vector values(maxval, -1); + foreach_permutation_impl(0, values, func); +} + +//----------------------------------------------------------------------------- +void +test_tiled_adjsets(const int dims[3], const std::string &testName) +{ + // Make 4 domains but alter their locations. + int index = 0; + foreach_permutation(4, [&](const std::vector &domainNumbering) + { + const int domains[] = {2,2,1}; + const int par_rank = relay::mpi::rank(MPI_COMM_WORLD); + const std::vector reorder{"normal", "kdtree"}; + for(const auto &r : reorder) + { + // Make the mesh. + conduit::Node mesh; + make_tiled(mesh, dims, domains, r, domainNumbering); + + // Now, make a corner mesh + conduit::Node s2dmap, d2smap; + conduit::blueprint::mpi::mesh::generate_corners(mesh, + "mesh_adjset", + "corner_adjset", + "corner_mesh", + "corner_coords", + s2dmap, + d2smap, + MPI_COMM_WORLD); + // Convert to pairwise adjset. + std::vector doms = conduit::blueprint::mesh::domains(mesh); + for(auto dom_ptr : doms) + { + conduit::Node &domain = *dom_ptr; + conduit::blueprint::mesh::adjset::to_pairwise(domain["adjsets/corner_adjset"], + domain["adjsets/corner_pairwise_adjset"]); + } + +#ifdef CONDUIT_WRITE_TEST_DATA + // Save the mesh. + std::stringstream ss; + ss << "_r" << r; + for(const auto &value : domainNumbering) + ss << "_" << value; + std::string filebase(testName + ss.str()); + save_mesh(mesh, filebase); +#endif + + // Check that its adjset points are the same along the edges. + bool same = conduit::blueprint::mpi::mesh::utils::adjset::compare_pointwise(mesh, "mesh_adjset", MPI_COMM_WORLD); + if(!same) + { + mesh.print(); + } + EXPECT_TRUE(same); + + // Check that its adjset points are the same along the edges. + same = conduit::blueprint::mpi::mesh::utils::adjset::compare_pointwise(mesh, "corner_pairwise_adjset", MPI_COMM_WORLD); + if(!same) + { + mesh.print(); + } + EXPECT_TRUE(same); + } + }); +} + +//----------------------------------------------------------------------------- +TEST(conduit_blueprint_mpi_mesh_tiled, two_dimensional) +{ + const int dims[]= {3,3,0}; + test_tiled_adjsets(dims, "two_dimensional"); +} + +//----------------------------------------------------------------------------- +TEST(conduit_blueprint_mpi_mesh_tiled, three_dimensional) +{ + const int dims[]= {2,2,2}; + test_tiled_adjsets(dims, "three_dimensional"); +} + +//----------------------------------------------------------------------------- +TEST(conduit_blueprint_mpi_mesh_tiled, three_dimensional_12) +{ + // This 12 domain case was found to cause adjset problems. + const int dims[] = {2,2,2}, domains[] = {3,2,2}; + const int par_rank = relay::mpi::rank(MPI_COMM_WORLD); + const std::vector reorder{"normal", "kdtree"}; + + for(const auto &r : reorder) + { + // Make the mesh. + conduit::Node mesh; + make_tiled(mesh, dims, domains, r, std::vector{}); + + // Now, make a corner mesh + conduit::Node s2dmap, d2smap; + conduit::blueprint::mpi::mesh::generate_corners(mesh, + "mesh_adjset", + "corner_adjset", + "corner_mesh", + "corner_coords", + s2dmap, + d2smap, + MPI_COMM_WORLD); + // Convert to pairwise adjset. + std::vector doms = conduit::blueprint::mesh::domains(mesh); + for(auto dom_ptr : doms) + { + conduit::Node &domain = *dom_ptr; + conduit::blueprint::mesh::adjset::to_pairwise(domain["adjsets/corner_adjset"], + domain["adjsets/corner_pairwise_adjset"]); + } + +#ifdef CONDUIT_WRITE_TEST_DATA + // Save the mesh. + const std::string testName("three_dimensional_12"); + std::stringstream ss; + ss << "_r_" << r; + std::string filebase(testName + ss.str()); + save_mesh(mesh, filebase); +#endif + + // Check that its adjset points are the same along the edges. + bool same = conduit::blueprint::mpi::mesh::utils::adjset::compare_pointwise(mesh, "corner_pairwise_adjset", MPI_COMM_WORLD); + EXPECT_TRUE(same); + } +} + +//----------------------------------------------------------------------------- +int main(int argc, char* argv[]) +{ + int result = 0; + + ::testing::InitGoogleTest(&argc, argv); + MPI_Init(&argc, &argv); + result = RUN_ALL_TESTS(); + MPI_Finalize(); + + return result; +} diff --git a/src/tests/blueprint/t_blueprint_mpi_mesh_transform.cpp b/src/tests/blueprint/t_blueprint_mpi_mesh_transform.cpp index ccf6b676a..98241b806 100644 --- a/src/tests/blueprint/t_blueprint_mpi_mesh_transform.cpp +++ b/src/tests/blueprint/t_blueprint_mpi_mesh_transform.cpp @@ -836,33 +836,14 @@ iterate_adjset(conduit::Node &mesh, const std::string &adjsetName, Func &&func) } //----------------------------------------------------------------------------- -TEST(conduit_blueprint_mesh_examples, generate_corners_wonky) +void +test_adjset_points(conduit::Node &mesh, const std::string &adjsetName) { - // There is a 3x3x3 zone mesh that was giving generate_corners a problem - // due to adjacency sets. Adjacency sets are produced from the original one - // and when making corners. It was giving rise to adjacency sets that - // contained points that do not exist in neighbor domains. We can use the - // PointQuery to test this. - - conduit::Node mesh, s2dmap, d2smap; - generate_wonky_mesh(mesh); - - // Make the corner mesh. - conduit::blueprint::mpi::mesh::generate_corners(mesh, - "main_adjset", - "corner_adjset", - "corner_mesh", - "corner_coords", - s2dmap, - d2smap, - MPI_COMM_WORLD); - - conduit::blueprint::mpi::mesh::utils::query::PointQuery Q(mesh, MPI_COMM_WORLD); // Iterate over the points in the adjset and add them to the - iterate_adjset(mesh, "corner_adjset", - [&](int /*dom*/, int nbr, int val, const conduit::Node *cset, const conduit::Node */*topo*/) + iterate_adjset(mesh, adjsetName, + [&](int /*dom*/, int nbr, int val, const conduit::Node *cset, const conduit::Node * /*topo*/) { // Get the point (it might not be 3D) auto pt = conduit::blueprint::mesh::utils::coordset::_explicit::coords(*cset, val); @@ -884,10 +865,73 @@ TEST(conduit_blueprint_mesh_examples, generate_corners_wonky) const auto &r = Q.results(domainId); auto it = std::find(r.begin(), r.end(), Q.NotFound); bool found = it != r.end(); + + // If NotFound was in the results, print the occurrances for debugging. + if(found) + { + const int par_rank = relay::mpi::rank(MPI_COMM_WORLD); + const auto &p = Q.inputs(domainId); + std::stringstream ss; + ss << "rank" << par_rank << ":\n"; + ss << " domain: " << domainId << "\n"; + ss << " search:\n"; + for(size_t i = 0; i < r.size(); i++) + { + if(r[i] == Q.NotFound) + { + ss << " -\n"; + ss << " point: [" << p[3*i+0] << ", " << p[3*i+1] << ", " << p[3*i+2] << "]\n"; + ss << " found: " << r[i] << "\n"; + } + } + std::cout << ss.str() << std::endl; + } + EXPECT_FALSE(found); } } +//----------------------------------------------------------------------------- +TEST(conduit_blueprint_mesh_examples, generate_corners_wonky) +{ + // There is a 3x3x3 zone mesh that was giving generate_corners a problem + // due to adjacency sets. Adjacency sets are produced from the original one + // and when making corners. It was giving rise to adjacency sets that + // contained points that do not exist in neighbor domains. We can use the + // PointQuery to test this. + + conduit::Node mesh, s2dmap, d2smap; + generate_wonky_mesh(mesh); + + // Make the corner mesh. + conduit::blueprint::mpi::mesh::generate_corners(mesh, + "main_adjset", + "corner_adjset", + "corner_mesh", + "corner_coords", + s2dmap, + d2smap, + MPI_COMM_WORLD); + + std::vector domains = conduit::blueprint::mesh::domains(mesh); + + // Test the adjset points. + test_adjset_points(mesh, "corner_adjset"); + + // Convert the adjset in all domains to pairwise and test the points again. + for(auto dom_ptr : domains) + { + Node &domain = *dom_ptr; + conduit::blueprint::mesh::adjset::to_pairwise(domain["adjsets/corner_adjset"], + domain["adjsets/corner_pairwise_adjset"]); + } + test_adjset_points(mesh, "corner_pairwise_adjset"); + + bool same_pointwise = conduit::blueprint::mpi::mesh::utils::adjset::compare_pointwise( + mesh, "corner_pairwise_adjset", MPI_COMM_WORLD); + EXPECT_TRUE(same_pointwise); +} + //----------------------------------------------------------------------------- TEST(conduit_blueprint_mesh_examples, generate_faces_wonky) {