Skip to content
Snippets Groups Projects
Commit 35d5c18d authored by Stéphane Del Pino's avatar Stéphane Del Pino
Browse files

Add a first stencil construction mechanism

- in parallel, ghost cells have no stencil
- the stencil is composed of cells that share a node with the
  considered cell
- a storage mechanism should be used to build the stencil once for all (soon?)
parent c5392ac9
No related branches found
No related tags found
1 merge request!205High-order polynomial reconstruction
...@@ -43,6 +43,7 @@ add_library( ...@@ -43,6 +43,7 @@ add_library(
MeshTransformer.cpp MeshTransformer.cpp
MeshUtils.cpp MeshUtils.cpp
MeshVariant.cpp MeshVariant.cpp
StencilBuilder.cpp
SynchronizerManager.cpp SynchronizerManager.cpp
) )
......
#ifndef STENCIL_HPP
#define STENCIL_HPP
#include <mesh/ConnectivityMatrix.hpp>
#include <mesh/ItemId.hpp>
class Stencil
{
private:
ConnectivityMatrix m_stencil;
public:
PUGS_INLINE
auto
operator[](CellId cell_id) const
{
return m_stencil[cell_id];
}
Stencil(const ConnectivityMatrix& stencil) : m_stencil(stencil) {}
Stencil(const Stencil&) = default;
Stencil(Stencil&&) = default;
~Stencil() = default;
};
#endif // STENCIL_HPP
#include <mesh/Connectivity.hpp>
#include <mesh/StencilBuilder.hpp>
template <typename ConnectivityType>
Array<const uint32_t>
StencilBuilder::_getRowMap(const ConnectivityType& connectivity) const
{
auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
auto cell_is_owned = connectivity.cellIsOwned();
Array<uint32_t> row_map{connectivity.numberOfCells() + 1};
row_map[0] = 0;
std::vector<CellId> neighbors;
for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
neighbors.resize(0);
// The stencil is not built for ghost cells
if (cell_is_owned[cell_id]) {
auto cell_nodes = cell_to_node_matrix[cell_id];
for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
const NodeId node_id = cell_nodes[i_node];
auto node_cells = node_to_cell_matrix[node_id];
for (size_t i_node_cell = 0; i_node_cell < node_cells.size(); ++i_node_cell) {
const CellId node_cell_id = node_cells[i_node_cell];
if (node_cell_id != cell_id) {
neighbors.push_back(node_cells[i_node_cell]);
}
}
}
std::sort(neighbors.begin(), neighbors.end());
neighbors.erase(std::unique(neighbors.begin(), neighbors.end()), neighbors.end());
}
// The cell itself is not counted
row_map[cell_id + 1] = row_map[cell_id] + neighbors.size();
}
return row_map;
}
template <typename ConnectivityType>
Array<const uint32_t>
StencilBuilder::_getColumnIndices(const ConnectivityType& connectivity, const Array<const uint32_t>& row_map) const
{
auto cell_number = connectivity.cellNumber();
Array<uint32_t> max_index(row_map.size() - 1);
parallel_for(
max_index.size(), PUGS_LAMBDA(size_t i) { max_index[i] = row_map[i]; });
auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
auto cell_is_owned = connectivity.cellIsOwned();
Array<uint32_t> column_indices(row_map[row_map.size() - 1]);
column_indices.fill(std::numeric_limits<uint32_t>::max());
for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
// The stencil is not built for ghost cells
if (cell_is_owned[cell_id]) {
auto cell_nodes = cell_to_node_matrix[cell_id];
for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
const NodeId node_id = cell_nodes[i_node];
auto node_cells = node_to_cell_matrix[node_id];
for (size_t i_node_cell = 0; i_node_cell < node_cells.size(); ++i_node_cell) {
const CellId node_cell_id = node_cells[i_node_cell];
if (node_cell_id != cell_id) {
bool found = false;
for (size_t i_index = row_map[cell_id]; i_index < max_index[cell_id]; ++i_index) {
if (column_indices[i_index] == node_cell_id) {
found = true;
break;
}
}
if (not found) {
int node_cell_number = cell_number[node_cell_id];
size_t i_index = row_map[cell_id];
// search for position for index
while ((i_index < max_index[cell_id])) {
if (node_cell_number > cell_number[CellId(column_indices[i_index])]) {
++i_index;
} else {
break;
}
}
for (size_t i_destination = max_index[cell_id]; i_destination > i_index; --i_destination) {
size_t i_source = i_destination - 1;
column_indices[i_destination] = column_indices[i_source];
}
++max_index[cell_id];
column_indices[i_index] = node_cell_id;
}
}
}
}
}
}
return column_indices;
}
template <typename ConnectivityType>
Stencil
StencilBuilder::build(const ConnectivityType& connectivity) const
{
Array<const uint32_t> row_map = this->_getRowMap(connectivity);
Array<const uint32_t> column_indices = this->_getColumnIndices(connectivity, row_map);
return ConnectivityMatrix{row_map, column_indices};
}
template Stencil StencilBuilder::build(const Connectivity<1>& connectivity) const;
template Stencil StencilBuilder::build(const Connectivity<2>& connectivity) const;
template Stencil StencilBuilder::build(const Connectivity<3>& connectivity) const;
#ifndef STENCIL_BUILDER_HPP
#define STENCIL_BUILDER_HPP
#include <mesh/Stencil.hpp>
class StencilBuilder
{
private:
template <typename ConnectivityType>
Array<const uint32_t> _getRowMap(const ConnectivityType& connectivity) const;
template <typename ConnectivityType>
Array<const uint32_t> _getColumnIndices(const ConnectivityType& connectivity,
const Array<const uint32_t>& row_map) const;
public:
template <typename ConnectivityType>
Stencil build(const ConnectivityType& connectivity) const;
StencilBuilder() = default;
StencilBuilder(const StencilBuilder&) = default;
StencilBuilder(StencilBuilder&&) = default;
~StencilBuilder() = default;
};
#endif // STENCIL_BUILDER_HPP
...@@ -204,6 +204,7 @@ add_executable (mpi_unit_tests ...@@ -204,6 +204,7 @@ add_executable (mpi_unit_tests
test_ParallelChecker_read.cpp test_ParallelChecker_read.cpp
test_Partitioner.cpp test_Partitioner.cpp
test_RandomEngine.cpp test_RandomEngine.cpp
test_StencilBuilder.cpp
test_SubItemArrayPerItemVariant.cpp test_SubItemArrayPerItemVariant.cpp
test_SubItemValuePerItem.cpp test_SubItemValuePerItem.cpp
test_SubItemValuePerItemVariant.cpp test_SubItemValuePerItemVariant.cpp
......
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_all.hpp>
#include <MeshDataBaseForTests.hpp>
#include <mesh/Connectivity.hpp>
#include <mesh/ConnectivityUtils.hpp>
#include <mesh/ItemValue.hpp>
#include <mesh/ItemValueUtils.hpp>
#include <mesh/Mesh.hpp>
#include <mesh/MeshVariant.hpp>
#include <mesh/StencilBuilder.hpp>
#include <utils/Messenger.hpp>
// clazy:excludeall=non-pod-global-static
TEST_CASE("StencilBuilder", "[mesh]")
{
auto is_valid = [](const auto& connectivity, const Stencil& stencil) {
auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
auto cell_is_owned = connectivity.cellIsOwned();
auto cell_number = connectivity.cellNumber();
for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
if (cell_is_owned[cell_id]) {
std::set<CellId, std::function<bool(CellId, CellId)>> cell_set(
[=](CellId cell_0, CellId cell_1) { return cell_number[cell_0] < cell_number[cell_1]; });
auto cell_nodes = cell_to_node_matrix[cell_id];
for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
const NodeId node_id = cell_nodes[i_node];
auto node_cells = node_to_cell_matrix[node_id];
for (size_t i_node_cell = 0; i_node_cell < node_cells.size(); ++i_node_cell) {
const CellId node_cell_id = node_cells[i_node_cell];
if (node_cell_id != cell_id) {
cell_set.insert(node_cell_id);
}
}
}
auto cell_stencil = stencil[cell_id];
auto i_set_cell = cell_set.begin();
for (size_t index = 0; index < cell_stencil.size(); ++index, ++i_set_cell) {
if (*i_set_cell != cell_stencil[index]) {
return false;
}
}
}
}
return true;
};
SECTION("1D")
{
SECTION("cartesian")
{
const auto& mesh = *MeshDataBaseForTests::get().cartesian1DMesh()->get<Mesh<1>>();
const Connectivity<1>& connectivity = mesh.connectivity();
Stencil stencil = StencilBuilder{}.build(connectivity);
REQUIRE(is_valid(connectivity, stencil));
}
SECTION("unordered")
{
const auto& mesh = *MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
const Connectivity<1>& connectivity = mesh.connectivity();
REQUIRE(is_valid(connectivity, StencilBuilder{}.build(connectivity)));
}
}
SECTION("2D")
{
SECTION("cartesian")
{
const auto& mesh = *MeshDataBaseForTests::get().cartesian2DMesh()->get<Mesh<2>>();
const Connectivity<2>& connectivity = mesh.connectivity();
REQUIRE(is_valid(connectivity, StencilBuilder{}.build(connectivity)));
}
SECTION("hybrid")
{
const auto& mesh = *MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>();
const Connectivity<2>& connectivity = mesh.connectivity();
REQUIRE(is_valid(connectivity, StencilBuilder{}.build(connectivity)));
}
}
SECTION("3D")
{
SECTION("carteian")
{
const auto& mesh = *MeshDataBaseForTests::get().cartesian3DMesh()->get<Mesh<3>>();
const Connectivity<3>& connectivity = mesh.connectivity();
REQUIRE(is_valid(connectivity, StencilBuilder{}.build(connectivity)));
}
SECTION("hybrid")
{
const auto& mesh = *MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>();
const Connectivity<3>& connectivity = mesh.connectivity();
REQUIRE(is_valid(connectivity, StencilBuilder{}.build(connectivity)));
}
}
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment