From 316d875b3c58cb4fe03e15924645cd6a82461e44 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Fri, 25 Oct 2024 19:20:09 +0200
Subject: [PATCH] [ci-skip] WIP: Rework stencil construction
---
src/mesh/StencilBuilder.cpp | 150 +++++++++++++++++++++++++++++++++++-
src/mesh/StencilBuilder.hpp | 3 +
2 files changed, 149 insertions(+), 4 deletions(-)
diff --git a/src/mesh/StencilBuilder.cpp b/src/mesh/StencilBuilder.cpp
index 1fa54e702..31db762a1 100644
--- a/src/mesh/StencilBuilder.cpp
+++ b/src/mesh/StencilBuilder.cpp
@@ -7,6 +7,104 @@
#include <set>
+template <ItemType item_type>
+class StencilBuilder::Layer
+{
+ using ItemId = ItemIdT<item_type>;
+ std::vector<ItemId> m_item_id_vector;
+ std::vector<int> m_item_number_vector;
+
+ public:
+ size_t
+ size() const
+ {
+ Assert(m_item_id_vector.size() == m_item_number_vector.size());
+ return m_item_id_vector.size();
+ }
+
+ bool
+ hasItemNumber(const int item_number) const
+ {
+ ssize_t begin = 0;
+ ssize_t end = m_item_number_vector.size();
+
+ while (begin < end) {
+ const ssize_t mid = (begin + end) / 2;
+ const auto& mid_number = m_item_number_vector[mid];
+ if (mid_number == item_number) {
+ return true; // We found the value
+ } else if (mid_number < item_number) {
+ if (begin == mid) {
+ break;
+ }
+ begin = mid - 1;
+ } else {
+ if (end == mid) {
+ break;
+ }
+ end = mid + 1;
+ }
+ }
+ return false;
+ }
+
+ const std::vector<ItemId>&
+ itemIdList() const
+ {
+ return m_item_id_vector;
+ }
+
+ void
+ add(const ItemId item_id, const int item_number)
+ {
+ ssize_t begin = 0;
+ ssize_t end = m_item_number_vector.size();
+
+ while (begin < end) {
+ const ssize_t mid = (begin + end) / 2;
+ const auto& mid_number = m_item_number_vector[mid];
+
+ if (mid_number == item_number) {
+ return; // We found the value
+ } else if (mid_number < item_number) {
+ if (begin == mid) {
+ break;
+ }
+ begin = mid;
+ } else {
+ if (end == mid) {
+ break;
+ }
+ end = mid;
+ }
+ }
+
+ m_item_id_vector.push_back(item_id);
+ m_item_number_vector.push_back(item_number);
+
+ const auto& begin_number = m_item_number_vector[begin];
+
+ if (begin_number > item_number) {
+ for (ssize_t i = m_item_number_vector.size() - 2; i >= begin; --i) {
+ std::swap(m_item_number_vector[i], m_item_number_vector[i + 1]);
+ std::swap(m_item_id_vector[i], m_item_id_vector[i + 1]);
+ }
+ } else if (begin_number < item_number) {
+ for (ssize_t i = m_item_number_vector.size() - 2; i > begin; --i) {
+ std::swap(m_item_number_vector[i], m_item_number_vector[i + 1]);
+ std::swap(m_item_id_vector[i], m_item_id_vector[i + 1]);
+ }
+ }
+ }
+
+ Layer() = default;
+
+ Layer(const Layer&) = default;
+ Layer(Layer&&) = default;
+
+ ~Layer() = default;
+};
+
template <ItemType connecting_item_type, typename ConnectivityType>
auto
StencilBuilder::_buildSymmetryConnectingItemList(const ConnectivityType& connectivity,
@@ -203,16 +301,60 @@ StencilBuilder::_build_for_same_source_and_target(const ConnectivityType& connec
auto item_to_connecting_item_matrix = connectivity.template getItemToItemMatrix<item_type, connecting_item_type>();
auto connecting_item_to_item_matrix = connectivity.template getItemToItemMatrix<connecting_item_type, item_type>();
- auto item_is_owned = connectivity.template isOwned<item_type>();
- auto item_number = connectivity.template number<item_type>();
+ auto item_is_owned = connectivity.template isOwned<item_type>();
+ auto item_number = connectivity.template number<item_type>();
+ auto connecting_item_number = connectivity.template number<connecting_item_type>();
using ItemId = ItemIdT<item_type>;
using ConnectingItemId = ItemIdT<connecting_item_type>;
if (symmetry_boundary_descriptor_list.size() == 0) {
if (number_of_layers == 1) {
- Array<const uint32_t> row_map = this->_getRowMap(connectivity);
- Array<const uint32_t> column_indices = this->_getColumnIndices(connectivity, row_map);
+ Array<uint32_t> row_map{connectivity.template numberOf<item_type>() + 1};
+ row_map[0] = 0;
+
+ std::vector<ItemId> column_indices_vector;
+
+ for (ItemId item_id = 0; item_id < connectivity.template numberOf<item_type>(); ++item_id) {
+ // First layer is a special case
+
+ Layer<item_type> item_layer;
+ Layer<connecting_item_type> connecting_layer;
+
+ if (item_is_owned[item_id]) {
+ for (size_t i_connecting_item_1 = 0; i_connecting_item_1 < item_to_connecting_item_matrix[item_id].size();
+ ++i_connecting_item_1) {
+ const ConnectingItemId layer_1_connecting_item_id =
+ item_to_connecting_item_matrix[item_id][i_connecting_item_1];
+ connecting_layer.add(layer_1_connecting_item_id, connecting_item_number[layer_1_connecting_item_id]);
+ }
+
+ for (auto connecting_item_id : connecting_layer.itemIdList()) {
+ for (size_t i_item_1 = 0; i_item_1 < connecting_item_to_item_matrix[connecting_item_id].size();
+ ++i_item_1) {
+ const ItemId layer_1_item_id = connecting_item_to_item_matrix[connecting_item_id][i_item_1];
+ if (layer_1_item_id != item_id) {
+ item_layer.add(layer_1_item_id, item_number[layer_1_item_id]);
+ }
+ }
+ }
+ }
+
+ for (auto layer_item_id : item_layer.itemIdList()) {
+ column_indices_vector.push_back(layer_item_id);
+ }
+ row_map[item_id + 1] = row_map[item_id] + item_layer.itemIdList().size();
+ }
+
+ if (row_map[row_map.size() - 1] != column_indices_vector.size()) {
+ throw UnexpectedError("incorrect stencil size");
+ }
+ Array<uint32_t> column_indices(row_map[row_map.size() - 1]);
+ column_indices.fill(std::numeric_limits<uint32_t>::max());
+
+ for (size_t i = 0; i < column_indices.size(); ++i) {
+ column_indices[i] = column_indices_vector[i];
+ }
return {ConnectivityMatrix{row_map, column_indices}, {}};
} else {
diff --git a/src/mesh/StencilBuilder.hpp b/src/mesh/StencilBuilder.hpp
index 5bfd821ca..4f7ead00b 100644
--- a/src/mesh/StencilBuilder.hpp
+++ b/src/mesh/StencilBuilder.hpp
@@ -16,6 +16,9 @@ class StencilBuilder
using BoundaryDescriptorList = std::vector<std::shared_ptr<const IBoundaryDescriptor>>;
private:
+ template <ItemType item_type>
+ class Layer;
+
template <typename ConnectivityType>
Array<const uint32_t> _getRowMap(const ConnectivityType& connectivity) const;
--
GitLab