Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add neighboring node indices to dual layouts in Draco_Mesh #1078

Merged
merged 4 commits into from
Jun 30, 2021
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
145 changes: 72 additions & 73 deletions src/mesh/Draco_Mesh.cc
Original file line number Diff line number Diff line change
Expand Up @@ -87,10 +87,11 @@ Draco_Mesh::Draco_Mesh(
num_faces_per_cell_, cell_to_node_linkage_, num_nodes_per_face_per_cell_, side_node_count_,
side_to_node_linkage_, ghost_cell_type_, ghost_cell_to_node_linkage_);

// build cell-to-corner-cell layout
compute_node_to_cell_linkage(num_faces_per_cell_, cell_to_node_linkage_,
num_nodes_per_face_per_cell_, ghost_cell_type_,
ghost_cell_to_node_linkage_, global_node_number_);
// build cell-to-corner-cell layout (\todo: extend to 3D)
if (dimension_ == 2) {
compute_node_to_cell_linkage(ghost_cell_type_, ghost_cell_to_node_linkage_,
global_node_number_);
}
}

//------------------------------------------------------------------------------------------------//
Expand Down Expand Up @@ -434,64 +435,42 @@ Draco_Mesh::compute_node_vec_indx_map(const std::vector<unsigned> &indx_type,

//------------------------------------------------------------------------------------------------//
/*!
* \brief Build a cell-to-cell linkage across corners
* \brief Build a dual layout node-cell-(node neighbors) linkage across corners
*
* \param[in] num_faces_per_cell number of faces per cell.
* \param[in] cell_to_node_linkage serial map of cell to face to node indices.
* \param[in] num_nodes_per_face_per_cell number of nodes per face per cell.
* \param[in] ghost_cell_type number of common vertices per ghost cell (sharing a full face).
* \param[in] ghost_cell_to_node_linkage vertices in common per ghost cell (sharing a full face).
* \param[in] global_node_number vector indexed by local node with global node index as values.
*/
void Draco_Mesh::compute_node_to_cell_linkage(
const std::vector<unsigned> &num_faces_per_cell,
const std::vector<unsigned> &cell_to_node_linkage,
const std::vector<unsigned> &num_nodes_per_face_per_cell,
const std::vector<unsigned> &ghost_cell_type,
const std::vector<unsigned> &ghost_cell_to_node_linkage,
const std::vector<unsigned> &global_node_number) {

// (1a) create map of (single) nodes to set of cells

std::map<unsigned, std::set<unsigned>> node_to_cells;

// initialize cell face counter and cell-node iterator
unsigned cf_counter = 0;
unsigned cn_counter = 0;

// convert cell-node linkage to map of cell face to
// convert cell-node linkage to map of node to adjacent cells and corresponding adjacent nodes
for (unsigned cell = 0; cell < num_cells; ++cell) {
for (unsigned face = 0; face < num_faces_per_cell[cell]; ++face) {
for (unsigned node = 0; node < num_nodes_per_face_per_cell[cf_counter]; ++node) {

// add cell to the set for the node at this point in the cell-node linkage vector
node_to_cells[cell_to_node_linkage[cn_counter]].insert({cell});

// each

// increment cell node counter
cn_counter++;
}

// increment cell face counter
cf_counter++;
}
}

Check(cn_counter == safe_convert_from_size_t(cell_to_node_linkage.size()));
// condense layout at this cell to a vector of unique nodes (preserves ordering in 2D)
const std::vector<unsigned> cell_nodes = this->get_cell_nodes(cell);

// (1b) convert set of cell per node to vector ...
// get the size of the unique node vector
const size_t num_cell_nodes = cell_nodes.size();

for (unsigned node = 0; node < num_nodes; ++node) {
// append a cell-(node vector) pair for each node on this cell
for (size_t cnode = 0; cnode < num_cell_nodes; ++cnode) {

// get the set of local cells for the node
const std::set<unsigned> &nc_set = node_to_cells.at(node);
// get index of the current node considered
const unsigned node = cell_nodes[cnode];

// convert set to vector
node_to_cell_linkage[node] = std::vector<unsigned>(nc_set.begin(), nc_set.end());
// get indices of immediate node neighbors
const size_t cnode_nbr1 = (cnode + 1) % num_cell_nodes;
const size_t cnode_nbr2 = (cnode + num_cell_nodes - 1) % num_cell_nodes;
const std::array<unsigned, 2> node_nbrs = {cell_nodes[cnode_nbr1], cell_nodes[cnode_nbr2]};

// each node must be adjacent to at least one local cell
Check(node_to_cell_linkage[node].size() >= 1);
// add the pair of neighbor cells and node indices to the set for this node
node_to_cellnode_linkage[node].push_back(std::make_pair(cell, node_nbrs));
}
}

// avoid populating ghost node map if there are no faces that go off rank
Expand All @@ -518,15 +497,14 @@ void Draco_Mesh::compute_node_to_cell_linkage(
//----------------------------------------------------------------------------------------------//

// create map of global node to vector of adjacent ranks
std::map<unsigned, std::vector<unsigned>> global_node_to_local_cells;

// first append this ranks index onto the map
// initialize cell face counter and cell-node iterator
unsigned gcn_counter = 0;
std::map<unsigned, std::vector<CellNodes_Pair>> global_node_to_local_cellnodes;

// short-cut to number of ghost faces (i.e. cells across a face on a rank boundary)
const unsigned num_ghost_cells = safe_convert_from_size_t(ghost_cell_type.size());

// initialize cell face counter
unsigned gcn_counter = 0;

// create the pre-comm map of global ghost nodes to local cells
for (unsigned ghost = 0; ghost < num_ghost_cells; ++ghost) {
for (unsigned ghost_node = 0; ghost_node < ghost_cell_type[ghost]; ++ghost_node) {
Expand All @@ -536,36 +514,48 @@ void Draco_Mesh::compute_node_to_cell_linkage(
const unsigned global_node = global_node_number[local_node];

// set initial rank and local cell listing at this global node
global_node_to_local_cells[global_node] = node_to_cell_linkage.at(local_node);
global_node_to_local_cellnodes[global_node] = node_to_cellnode_linkage.at(local_node);

// increment ghost-cell-node-linkage counter
gcn_counter++;
}
}

// create a serial array of local cell indices (for cells with at least one off-rank vertex)
std::vector<unsigned> cells_per_serial;
for (const auto &global_node_cell_pair : global_node_to_local_cells) {
cells_per_serial.insert(cells_per_serial.end(), global_node_cell_pair.second.begin(),
global_node_cell_pair.second.end());
//----------------------------------------------------------------------------------------------//
// create a serial arrays of local cell indices and neighbor nodes
std::vector<unsigned> cellnodes_per_serial;
for (const auto &global_node_cellnode_pair : global_node_to_local_cellnodes) {

// get the number of cell neighbors
const size_t num_cell_nbrs = global_node_cellnode_pair.second.size();

// flatten cell+node neighbor data for each cell
for (size_t j = 0; j < num_cell_nbrs; ++j) {
cellnodes_per_serial.emplace_back(global_node_cellnode_pair.second[j].first);
cellnodes_per_serial.emplace_back(global_node_cellnode_pair.second[j].second[0]);
cellnodes_per_serial.emplace_back(global_node_cellnode_pair.second[j].second[1]);
}
}

Check(cellnodes_per_serial.size() % 3 == 0);

// get the serialized vector size
const size_t num_serial = cells_per_serial.size();
const size_t num_serial = cellnodes_per_serial.size() / 3;

//----------------------------------------------------------------------------------------------//
// initialize a serial array of global node indices
std::vector<unsigned> global_node_per_serial(num_serial);

// map global ghost node indices to serial index
size_t serial_count = 0;
for (const auto &global_node_cell_pair : global_node_to_local_cells) {
for (const auto &global_node_cellnode_pair : global_node_to_local_cellnodes) {

// get the number of local cells for this node
const size_t num_node_cells = global_node_cell_pair.second.size();
const size_t num_node_cells = global_node_cellnode_pair.second.size();

// set to the global node at the serial index for this local cell
for (size_t node_cell = 0; node_cell < num_node_cells; ++node_cell)
global_node_per_serial[serial_count + node_cell] = global_node_cell_pair.first;
global_node_per_serial[serial_count + node_cell] = global_node_cellnode_pair.first;

// increment count over serial index
serial_count += num_node_cells;
Expand All @@ -574,44 +564,53 @@ void Draco_Mesh::compute_node_to_cell_linkage(
// check that serial vector has been filled
Check(serial_count == num_serial);

// gather the local cell indices per serial index per rank
std::vector<std::vector<unsigned>> cells_per_serial_per_rank;
rtt_c4::indeterminate_allgatherv(cells_per_serial, cells_per_serial_per_rank);
//----------------------------------------------------------------------------------------------//
// gather the global indices per serial index per rank
std::vector<std::vector<unsigned>> global_node_per_serial_per_rank;
rtt_c4::indeterminate_allgatherv(global_node_per_serial, global_node_per_serial_per_rank);

// get the number of ranks in this communicator group
const unsigned num_ranks = rtt_c4::nodes();

// sanity-check gatherv communicator group
Check(cells_per_serial_per_rank.size() == num_ranks);
Check(global_node_per_serial_per_rank.size() == num_ranks);

// resize gather target for global nodes since the sizes are determined from previous gather
std::vector<std::vector<unsigned>> global_node_per_serial_per_rank(num_ranks);
//----------------------------------------------------------------------------------------------//
// gather the local cell indices and associate neighbor nodes per serial index per rank

// resize gather target for neighbor data, since the sizes are determined from previous gather
std::vector<std::vector<unsigned>> cellnodes_per_serial_per_rank(num_ranks);
for (unsigned rank = 0; rank < num_ranks; ++rank)
global_node_per_serial_per_rank[rank].resize(cells_per_serial_per_rank[rank].size());
cellnodes_per_serial_per_rank[rank].resize(3 * global_node_per_serial_per_rank[rank].size());

// gather the global ghost node indices per serial index per rank
rtt_c4::determinate_allgatherv(global_node_per_serial, global_node_per_serial_per_rank);
rtt_c4::determinate_allgatherv(cellnodes_per_serial, cellnodes_per_serial_per_rank);

//----------------------------------------------------------------------------------------------//
// merge global_node_per_serial_per_rank and cells_per_serial_per_rank into map per rank
std::vector<std::map<unsigned, std::vector<unsigned>>> global_node_cell_map_per_rank(num_ranks);

std::vector<std::map<unsigned, std::vector<CellNodes_Pair>>> ghost_dualmap_per_rank(num_ranks);
for (unsigned rank = 0; rank < num_ranks; ++rank) {

// short-cut to serialized vector size from rank
const size_t num_serial_on_rank = cells_per_serial_per_rank[rank].size();
Check(num_serial_on_rank == global_node_per_serial_per_rank[rank].size());
const size_t num_serial_on_rank = global_node_per_serial_per_rank[rank].size();
Check(3 * num_serial_on_rank == cellnodes_per_serial_per_rank[rank].size());

// generate map for the rank
for (size_t i = 0; i < num_serial_on_rank; ++i) {

// short-cut to key (global node index) and value modifier (local_cell on rank)
const unsigned global_node = global_node_per_serial_per_rank[rank][i];
const unsigned local_cell = cells_per_serial_per_rank[rank][i];
const unsigned local_cell = cellnodes_per_serial_per_rank[rank][3 * i];
const std::array<unsigned, 2> node_nbrs = {cellnodes_per_serial_per_rank[rank][3 * i + 1],
cellnodes_per_serial_per_rank[rank][3 * i + 2]};

// accumulate local cells (for rank) adjacent to this global node
global_node_cell_map_per_rank[rank][global_node].push_back(local_cell);
ghost_dualmap_per_rank[rank][global_node].push_back(std::make_pair(local_cell, node_nbrs));
}
}

//----------------------------------------------------------------------------------------------//
// invert local-to-global node index layout
std::map<unsigned, unsigned> global_to_local_node;
for (unsigned node = 0; node < num_nodes; ++node)
Expand Down Expand Up @@ -653,8 +652,8 @@ void Draco_Mesh::compute_node_to_cell_linkage(
const unsigned node = global_to_local_node.at(gl_node);

// append each local-cell-rank pair to dual ghost layout
for (auto local_cell : global_node_cell_map_per_rank[rank].at(gl_node))
node_to_ghost_cell_linkage[node].push_back(std::make_pair(local_cell, rank));
for (auto local_cellnodes : ghost_dualmap_per_rank[rank].at(gl_node))
node_to_ghost_cell_linkage[node].push_back(std::make_pair(local_cellnodes, rank));
}
}

Expand Down
29 changes: 16 additions & 13 deletions src/mesh/Draco_Mesh.hh
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@

#include "ds++/config.h"
#include "mesh_element/Geometry.hh"
#include <array>
#include <map>
#include <set>
#include <vector>
Expand All @@ -28,12 +29,15 @@ namespace rtt_mesh {
* (cell adjacency) information. This class also provides basic services, including access to cell
* information. This mesh is based on an unstructured mesh implementation by Kent Budge.
*
* Two important features for a fully realized Draco_Mesh are the following:
* 1) Layout, which stores cell connectivity and hence the mesh topology.
* Important features for a fully realized Draco_Mesh are the following:
* 1) Geometry, which implies a metric for distance between points.
* 2) Layout, which stores cell connectivity and hence the mesh topology.
* a) It has an internal layout containing local cell-to-cell linkage,
* b) and a boundary layout with side off-process linkage, and
* b) a boundary layout with true-boundary linkage, and
* c) a ghost layout containing cell-to-ghost-cell linkage.
* 2) Geometry, which implies a metric for distance between points.
* 3) Dual_Layout, which stores node connectivity and effectively inverts Layout
* 4) Dual_Ghost_Layout, which stores node connectivity to off-process adjacent cells and nodes.
* This an has additional field for the MPI rank index the neighboring cell and nodes are on.
*
* Possibly temporary features:
* 1) The num_faces_per_cell_ vector (argument to the constructor) is currently taken to be the
Expand All @@ -53,12 +57,14 @@ public:
using Layout =
std::map<unsigned int, std::vector<std::pair<unsigned int, std::vector<unsigned int>>>>;

// e.g.: (key: node, value: vector of local or ghost cell indices)
using Dual_Layout = std::map<unsigned int, std::vector<unsigned int>>;
// e.g.: (key: node, value: vector of local or ghost cell indices), only 2D
using CellNodes_Pair = std::pair<unsigned int, std::array<unsigned int, 2>>;
using Dual_Layout = std::map<unsigned int, std::vector<CellNodes_Pair>>;

// e.g.: (key: node, value: vector of pairs of rank and local cell index on the rank)
// vectors of pairs are ordered in increasing rank, by construction (see Draco_Mesh.cc)
using Dual_Ghost_Layout =
std::map<unsigned int, std::vector<std::pair<unsigned int, unsigned int>>>;
std::map<unsigned int, std::vector<std::pair<CellNodes_Pair, unsigned int>>>;

protected:
// >>> DATA
Expand Down Expand Up @@ -107,7 +113,7 @@ protected:
Layout cell_to_ghost_cell_linkage;

// Node map to vector of local cells
Dual_Layout node_to_cell_linkage;
Dual_Layout node_to_cellnode_linkage;

// Node map to vector of ghost cells
Dual_Ghost_Layout node_to_ghost_cell_linkage;
Expand Down Expand Up @@ -150,7 +156,7 @@ public:
Layout get_cc_linkage() const { return cell_to_cell_linkage; }
Layout get_cs_linkage() const { return cell_to_side_linkage; }
Layout get_cg_linkage() const { return cell_to_ghost_cell_linkage; }
Dual_Layout get_nc_linkage() const { return node_to_cell_linkage; }
Dual_Layout get_nc_linkage() const { return node_to_cellnode_linkage; }
Dual_Ghost_Layout get_ngc_linkage() const { return node_to_ghost_cell_linkage; }

// >>> SERVICES
Expand Down Expand Up @@ -191,10 +197,7 @@ private:
const std::vector<unsigned> &indx_to_node_linkage) const;

//! Calculate cell-corner-cell layouts (adjacent cells not sharing a face)
void compute_node_to_cell_linkage(const std::vector<unsigned> &num_faces_per_cell,
const std::vector<unsigned> &cell_to_node_linkage,
const std::vector<unsigned> &num_nodes_per_face_per_cell,
const std::vector<unsigned> &ghost_cell_type,
void compute_node_to_cell_linkage(const std::vector<unsigned> &ghost_cell_type,
const std::vector<unsigned> &ghost_cell_to_node_linkage,
const std::vector<unsigned> &global_node_number);
};
Expand Down
29 changes: 23 additions & 6 deletions src/mesh/test/tstDraco_Mesh.cc
Original file line number Diff line number Diff line change
Expand Up @@ -291,17 +291,34 @@ void cartesian_mesh_2d(rtt_c4::ParallelUnitTest &ut) {
FAIL_IF(nc_layout.at(node).size() == 3);
FAIL_IF_NOT(nc_layout.at(node).size() <= 4);

// these check for an emergent orthogonal cell stride in the dual layout
// these check for emergent orthogonal strides in the dual layout
if (nc_layout.at(node).size() == 4) {
FAIL_IF_NOT(nc_layout.at(node)[1] == nc_layout.at(node)[0] + 1);
FAIL_IF_NOT(nc_layout.at(node)[3] == nc_layout.at(node)[2] + 1);
FAIL_IF_NOT(nc_layout.at(node)[2] == nc_layout.at(node)[0] + num_xdir);

// check cell stride
FAIL_IF_NOT(nc_layout.at(node)[1].first == nc_layout.at(node)[0].first + 1);
FAIL_IF_NOT(nc_layout.at(node)[3].first == nc_layout.at(node)[2].first + 1);
FAIL_IF_NOT(nc_layout.at(node)[2].first == nc_layout.at(node)[0].first + num_xdir);

// check node stride
FAIL_IF_NOT(nc_layout.at(node)[1].second[0] == nc_layout.at(node)[0].second[1]);
FAIL_IF_NOT(nc_layout.at(node)[1].second[1] == nc_layout.at(node)[3].second[0]);
FAIL_IF_NOT(nc_layout.at(node)[2].second[0] == nc_layout.at(node)[3].second[1]);
FAIL_IF_NOT(nc_layout.at(node)[2].second[1] == nc_layout.at(node)[0].second[0]);
}

if (nc_layout.at(node).size() == 2) {
const bool has_xdir_stride = nc_layout.at(node)[1] == nc_layout.at(node)[0] + 1;
const bool has_ydir_stride = nc_layout.at(node)[1] == nc_layout.at(node)[0] + num_xdir;

// check cell stride
const bool has_xdir_stride = nc_layout.at(node)[1].first == nc_layout.at(node)[0].first + 1;
const bool has_ydir_stride =
nc_layout.at(node)[1].first == nc_layout.at(node)[0].first + num_xdir;
FAIL_IF_NOT(has_xdir_stride || has_ydir_stride);

// check node stride
const bool has_node_stride =
nc_layout.at(node)[1].second[1] == nc_layout.at(node)[0].second[0] ||
nc_layout.at(node)[1].second[0] == nc_layout.at(node)[0].second[1];
FAIL_IF_NOT(has_node_stride);
}
}

Expand Down
Loading