diff --git a/src/mesh/Stencil.hpp b/src/mesh/Stencil.hpp
index d67aa09ed03252b5085c9238d4e968c5947c2678..5dc20dc326fe0ebc716df1b85fb766f50cc790e0 100644
--- a/src/mesh/Stencil.hpp
+++ b/src/mesh/Stencil.hpp
@@ -8,8 +8,17 @@ class Stencil
{
private:
ConnectivityMatrix m_stencil;
+#warning TEMPORARY IMPLEMENTATION
+ std::vector<std::pair<std::string, ConnectivityMatrix>> m_symmetry_name_stencil;
public:
+ PUGS_INLINE
+ const auto&
+ symmetryNameStencil() const
+ {
+ return m_symmetry_name_stencil;
+ }
+
PUGS_INLINE
auto
operator[](CellId cell_id) const
@@ -17,7 +26,10 @@ class Stencil
return m_stencil[cell_id];
}
- Stencil(const ConnectivityMatrix& stencil) : m_stencil(stencil) {}
+ Stencil(const ConnectivityMatrix& stencil,
+ const std::vector<std::pair<std::string, ConnectivityMatrix>>& symmetry_name_stencil)
+ : m_stencil{stencil}, m_symmetry_name_stencil{symmetry_name_stencil}
+ {}
Stencil(const Stencil&) = default;
Stencil(Stencil&&) = default;
~Stencil() = default;
diff --git a/src/mesh/StencilBuilder.cpp b/src/mesh/StencilBuilder.cpp
index 34015395ae6d57c5e0c2850333ad90ef65211ef1..92738a6432fcaf5e86a5e1a5c9cd1d5f467b053d 100644
--- a/src/mesh/StencilBuilder.cpp
+++ b/src/mesh/StencilBuilder.cpp
@@ -1,6 +1,7 @@
#include <mesh/Connectivity.hpp>
#include <mesh/StencilBuilder.hpp>
+#include <mesh/ItemArray.hpp>
#include <utils/Messenger.hpp>
#warning remove after rework
@@ -109,91 +110,238 @@ StencilBuilder::_getColumnIndices(const ConnectivityType& connectivity, const Ar
template <typename ConnectivityType>
Stencil
-StencilBuilder::_build(const ConnectivityType& connectivity, size_t number_of_layers) const
+StencilBuilder::_build(const ConnectivityType& connectivity,
+ size_t number_of_layers,
+ const std::set<std::string>& symmetry_boundary_names) const
{
- 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);
-
- return ConnectivityMatrix{row_map, column_indices};
- } else {
- if (number_of_layers > 2) {
- throw NotImplementedError("number of layers too large");
- }
- if (parallel::size() > 1) {
- throw NotImplementedError("stencils with more than 1 layers are not defined in parallel");
- }
+#warning TEMPORARY
+ if (symmetry_boundary_names.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);
+
+ return {ConnectivityMatrix{row_map, column_indices}, {}};
+ } else {
+ if (number_of_layers > 2) {
+ throw NotImplementedError("number of layers too large");
+ }
+ if (parallel::size() > 1) {
+ throw NotImplementedError("stencils with more than 1 layers are not defined in parallel");
+ }
- auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
- auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
+ 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();
+ auto cell_is_owned = connectivity.cellIsOwned();
+ auto cell_number = connectivity.cellNumber();
- Array<uint32_t> row_map{connectivity.numberOfCells() + 1};
- row_map[0] = 0;
+ Array<uint32_t> row_map{connectivity.numberOfCells() + 1};
+ row_map[0] = 0;
- std::vector<CellId> column_indices_vector;
+ std::vector<CellId> column_indices_vector;
- 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]; });
+ 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]; });
- for (size_t i_node_1 = 0; i_node_1 < cell_to_node_matrix[cell_id].size(); ++i_node_1) {
- const NodeId layer_1_node_id = cell_to_node_matrix[cell_id][i_node_1];
+ for (size_t i_node_1 = 0; i_node_1 < cell_to_node_matrix[cell_id].size(); ++i_node_1) {
+ const NodeId layer_1_node_id = cell_to_node_matrix[cell_id][i_node_1];
- for (size_t i_cell_1 = 0; i_cell_1 < node_to_cell_matrix[layer_1_node_id].size(); ++i_cell_1) {
- CellId cell_1_id = node_to_cell_matrix[layer_1_node_id][i_cell_1];
+ for (size_t i_cell_1 = 0; i_cell_1 < node_to_cell_matrix[layer_1_node_id].size(); ++i_cell_1) {
+ CellId cell_1_id = node_to_cell_matrix[layer_1_node_id][i_cell_1];
- for (size_t i_node_2 = 0; i_node_2 < cell_to_node_matrix[cell_1_id].size(); ++i_node_2) {
- const NodeId layer_2_node_id = cell_to_node_matrix[cell_1_id][i_node_2];
+ for (size_t i_node_2 = 0; i_node_2 < cell_to_node_matrix[cell_1_id].size(); ++i_node_2) {
+ const NodeId layer_2_node_id = cell_to_node_matrix[cell_1_id][i_node_2];
- for (size_t i_cell_2 = 0; i_cell_2 < node_to_cell_matrix[layer_2_node_id].size(); ++i_cell_2) {
- CellId cell_2_id = node_to_cell_matrix[layer_2_node_id][i_cell_2];
+ for (size_t i_cell_2 = 0; i_cell_2 < node_to_cell_matrix[layer_2_node_id].size(); ++i_cell_2) {
+ CellId cell_2_id = node_to_cell_matrix[layer_2_node_id][i_cell_2];
- if (cell_2_id != cell_id) {
- cell_set.insert(cell_2_id);
+ if (cell_2_id != cell_id) {
+ cell_set.insert(cell_2_id);
+ }
}
}
}
}
- }
- for (auto stencil_cell_id : cell_set) {
- column_indices_vector.push_back(stencil_cell_id);
+ for (auto stencil_cell_id : cell_set) {
+ column_indices_vector.push_back(stencil_cell_id);
+ }
+ row_map[cell_id + 1] = row_map[cell_id] + cell_set.size();
}
- row_map[cell_id + 1] = row_map[cell_id] + cell_set.size();
}
- }
- if (row_map[row_map.size() - 1] != column_indices_vector.size()) {
- throw UnexpectedError("incorrect stencil 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());
- 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];
+ }
+ ConnectivityMatrix primal_stencil{row_map, column_indices};
- for (size_t i = 0; i < column_indices.size(); ++i) {
- column_indices[i] = column_indices_vector[i];
+ return {primal_stencil, {}};
}
+ } else {
+ if constexpr (ConnectivityType::Dimension > 1) {
+ std::vector<Array<const FaceId>> boundary_node_list;
+
+ NodeArray<bool> symmetry_node_list(connectivity, symmetry_boundary_names.size());
+ symmetry_node_list.fill(0);
+
+ auto face_to_node_matrix = connectivity.faceToNodeMatrix();
+ auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
+ auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
+
+ {
+ size_t i_boundary_name = 0;
+ for (auto boundary_name : symmetry_boundary_names) {
+ bool found = false;
+ for (size_t i_ref_node_list = 0;
+ i_ref_node_list < connectivity.template numberOfRefItemList<ItemType::face>(); ++i_ref_node_list) {
+ const auto& ref_face_list = connectivity.template refItemList<ItemType::face>(i_ref_node_list);
+ if (ref_face_list.refId().tagName() == boundary_name) {
+ found = true;
+ boundary_node_list.push_back(ref_face_list.list());
+ for (size_t i_face = 0; i_face < ref_face_list.list().size(); ++i_face) {
+ const FaceId face_id = ref_face_list.list()[i_face];
+ auto node_list = face_to_node_matrix[face_id];
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ const NodeId node_id = node_list[i_node];
+
+ symmetry_node_list[node_id][i_boundary_name] = true;
+ }
+ }
+ break;
+ }
+ }
+ ++i_boundary_name;
+ if (not found) {
+#warning IMPROVE ERROR MESSAGE
+ throw NormalError("cannot find boundary " + boundary_name);
+ };
+ }
+ }
+
+ auto cell_is_owned = connectivity.cellIsOwned();
+ auto cell_number = connectivity.cellNumber();
- return ConnectivityMatrix{row_map, column_indices};
+ Array<uint32_t> row_map{connectivity.numberOfCells() + 1};
+ row_map[0] = 0;
+ std::vector<Array<uint32_t>> symmetry_row_map_list(symmetry_boundary_names.size());
+ for (auto&& symmetry_row_map : symmetry_row_map_list) {
+ symmetry_row_map = Array<uint32_t>{connectivity.numberOfCells() + 1};
+ symmetry_row_map[0] = 0;
+ }
+
+ std::vector<uint32_t> column_indices_vector;
+ std::vector<std::vector<uint32_t>> symmetry_column_indices_vector(symmetry_boundary_names.size());
+
+ for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
+ std::set<CellId> cell_set;
+
+ if (cell_is_owned[cell_id]) {
+ auto cell_node_list = cell_to_node_matrix[cell_id];
+ for (size_t i_cell_node = 0; i_cell_node < cell_node_list.size(); ++i_cell_node) {
+ const NodeId cell_node_id = cell_node_list[i_cell_node];
+ auto node_cell_list = node_to_cell_matrix[cell_node_id];
+ for (size_t i_node_cell = 0; i_node_cell < node_cell_list.size(); ++i_node_cell) {
+ const CellId node_cell_id = node_cell_list[i_node_cell];
+ if (cell_id != node_cell_id) {
+ cell_set.insert(node_cell_id);
+ }
+ }
+ }
+ row_map[cell_id + 1] = row_map[cell_id] + cell_set.size();
+
+ {
+ std::vector<CellId> cell_vector;
+ for (auto&& set_cell_id : cell_set) {
+ cell_vector.push_back(set_cell_id);
+ }
+ std::sort(cell_vector.begin(), cell_vector.end(),
+ [&cell_number](const CellId& cell0_id, const CellId& cell1_id) {
+ return cell_number[cell0_id] < cell_number[cell1_id];
+ });
+
+ for (auto&& vector_cell_id : cell_vector) {
+ column_indices_vector.push_back(vector_cell_id);
+ }
+ }
+
+ std::vector<std::set<CellId>> by_boundary_symmetry_cell;
+ for (size_t i = 0; i < symmetry_boundary_names.size(); ++i) {
+ std::set<CellId> symmetry_cell_set;
+ for (size_t i_cell_node = 0; i_cell_node < cell_node_list.size(); ++i_cell_node) {
+ const NodeId cell_node_id = cell_node_list[i_cell_node];
+ if (symmetry_node_list[cell_node_id][i]) {
+ auto node_cell_list = node_to_cell_matrix[cell_node_id];
+ for (size_t i_node_cell = 0; i_node_cell < node_cell_list.size(); ++i_node_cell) {
+ const CellId node_cell_id = node_cell_list[i_node_cell];
+ symmetry_cell_set.insert(node_cell_id);
+ }
+ }
+ }
+ symmetry_row_map_list[i][cell_id + 1] = symmetry_row_map_list[i][cell_id] + symmetry_cell_set.size();
+ by_boundary_symmetry_cell.emplace_back(symmetry_cell_set);
+
+ std::vector<CellId> cell_vector;
+ for (auto&& set_cell_id : symmetry_cell_set) {
+ cell_vector.push_back(set_cell_id);
+ }
+ std::sort(cell_vector.begin(), cell_vector.end(),
+ [&cell_number](const CellId& cell0_id, const CellId& cell1_id) {
+ return cell_number[cell0_id] < cell_number[cell1_id];
+ });
+
+ for (auto&& vector_cell_id : cell_vector) {
+ symmetry_column_indices_vector[i].push_back(vector_cell_id);
+ }
+ }
+ }
+ }
+ ConnectivityMatrix primal_stencil{row_map, convert_to_array(column_indices_vector)};
+
+ std::vector<std::pair<std::string, ConnectivityMatrix>> symmetry_name_stencil;
+ {
+ size_t i = 0;
+ for (auto&& boundary_name : symmetry_boundary_names) {
+ symmetry_name_stencil.emplace_back(
+ std::make_pair(boundary_name, ConnectivityMatrix{symmetry_row_map_list[i],
+ convert_to_array(symmetry_column_indices_vector[i])}));
+ ++i;
+ }
+ }
+
+ return {{primal_stencil}, {symmetry_name_stencil}};
+
+ } else {
+ throw NotImplementedError("Only implemented in 2D");
+ }
}
}
Stencil
-StencilBuilder::build(const IConnectivity& connectivity, size_t number_of_layers) const
+StencilBuilder::build(const IConnectivity& connectivity,
+ size_t number_of_layers,
+ const std::set<std::string>& symmetry_boundary_names) const
{
switch (connectivity.dimension()) {
case 1: {
- return StencilBuilder::_build(dynamic_cast<const Connectivity<1>&>(connectivity), number_of_layers);
+ return StencilBuilder::_build(dynamic_cast<const Connectivity<1>&>(connectivity), number_of_layers,
+ symmetry_boundary_names);
}
case 2: {
- return StencilBuilder::_build(dynamic_cast<const Connectivity<2>&>(connectivity), number_of_layers);
+ return StencilBuilder::_build(dynamic_cast<const Connectivity<2>&>(connectivity), number_of_layers,
+ symmetry_boundary_names);
}
case 3: {
- return StencilBuilder::_build(dynamic_cast<const Connectivity<3>&>(connectivity), number_of_layers);
+ return StencilBuilder::_build(dynamic_cast<const Connectivity<3>&>(connectivity), number_of_layers,
+ symmetry_boundary_names);
}
default: {
throw UnexpectedError("invalid connectivity dimension");
diff --git a/src/mesh/StencilBuilder.hpp b/src/mesh/StencilBuilder.hpp
index 53918e64ac9234c396986f2bfef3f90510f287e0..c52933da26b4969410eba7f8aca090842ba9fdc8 100644
--- a/src/mesh/StencilBuilder.hpp
+++ b/src/mesh/StencilBuilder.hpp
@@ -3,6 +3,9 @@
#include <mesh/Stencil.hpp>
+#include <set>
+#include <string>
+
class IConnectivity;
class StencilBuilder
{
@@ -15,10 +18,14 @@ class StencilBuilder
const Array<const uint32_t>& row_map) const;
template <typename ConnectivityType>
- Stencil _build(const ConnectivityType& connectivity, size_t number_of_layers) const;
+ Stencil _build(const ConnectivityType& connectivity,
+ size_t number_of_layers,
+ const std::set<std::string>& symmetry_boundary_names) const;
friend class StencilManager;
- Stencil build(const IConnectivity& connectivity, size_t number_of_layers) const;
+ Stencil build(const IConnectivity& connectivity,
+ size_t number_of_layers,
+ const std::set<std::string>& symmetry_boundary_names) const;
public:
StencilBuilder() = default;
diff --git a/src/mesh/StencilManager.cpp b/src/mesh/StencilManager.cpp
index 120a95e0dfc54a38f9ea7a4ab94bbe4607014981..618431ee75dab1696c4500a306dd70a9f99e1e7f 100644
--- a/src/mesh/StencilManager.cpp
+++ b/src/mesh/StencilManager.cpp
@@ -32,14 +32,16 @@ StencilManager::destroy()
}
const Stencil&
-StencilManager::getStencil(const IConnectivity& connectivity, size_t degree)
+StencilManager::getStencil(const IConnectivity& connectivity,
+ size_t degree,
+ const std::set<std::string>& symmetry_boundary_names)
{
- if (not m_stored_stencil_map.contains(Key{connectivity.id(), degree})) {
- m_stored_stencil_map[Key{connectivity.id(), degree}] =
- std::make_shared<Stencil>(StencilBuilder{}.build(connectivity, degree));
+ if (not m_stored_stencil_map.contains(Key{connectivity.id(), degree, symmetry_boundary_names})) {
+ m_stored_stencil_map[Key{connectivity.id(), degree, symmetry_boundary_names}] =
+ std::make_shared<Stencil>(StencilBuilder{}.build(connectivity, degree, symmetry_boundary_names));
}
- return *m_stored_stencil_map.at(Key{connectivity.id(), degree});
+ return *m_stored_stencil_map.at(Key{connectivity.id(), degree, symmetry_boundary_names});
}
void
diff --git a/src/mesh/StencilManager.hpp b/src/mesh/StencilManager.hpp
index f2912be1fca197aab6b9e010c600de3d8e507827..8199ce5a324a26879ca2c94c05f8871e9e6d7559 100644
--- a/src/mesh/StencilManager.hpp
+++ b/src/mesh/StencilManager.hpp
@@ -5,6 +5,7 @@
#include <mesh/Stencil.hpp>
#include <memory>
+#include <set>
#include <unordered_map>
class StencilManager
@@ -19,11 +20,26 @@ class StencilManager
{
size_t connectivity_id;
size_t degree;
+ std::set<std::string> symmetry_boundary_names;
PUGS_INLINE bool
operator==(const Key& k) const
{
- return (connectivity_id == k.connectivity_id) and (degree == k.degree);
+ if ((connectivity_id != k.connectivity_id) or (degree != k.degree) or
+ (symmetry_boundary_names.size() != k.symmetry_boundary_names.size())) {
+ return false;
+ }
+
+ auto i_name = symmetry_boundary_names.begin();
+ auto i_k_name = k.symmetry_boundary_names.begin();
+
+ for (; i_name != symmetry_boundary_names.end(); ++i_name, ++i_k_name) {
+ if (*i_name != *i_k_name) {
+ return false;
+ }
+ }
+
+ return true;
}
};
struct HashKey
@@ -31,6 +47,7 @@ class StencilManager
size_t
operator()(const Key& key) const
{
+ // We do not use the symmetry boundary set by now
return (std::hash<decltype(Key::connectivity_id)>()(key.connectivity_id)) ^
(std::hash<decltype(Key::degree)>()(key.degree) << 1);
}
@@ -52,7 +69,9 @@ class StencilManager
void deleteConnectivity(const size_t connectivity_id);
- const Stencil& getStencil(const IConnectivity& i_connectivity, size_t degree = 1);
+ const Stencil& getStencil(const IConnectivity& i_connectivity,
+ size_t degree = 1,
+ const std::set<std::string>& symmetry_boundary_names = {});
StencilManager(const StencilManager&) = delete;
StencilManager(StencilManager&&) = delete;
diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index 7673255bbb3c3e20266a07a0eb6889998c3ae25e..086583c9f5b5fcca55a9321c00a2f51f13934709 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -17,11 +17,29 @@
#include <geometry/TriangleTransformation.hpp>
#include <mesh/MeshData.hpp>
#include <mesh/MeshDataManager.hpp>
+#include <mesh/MeshFlatFaceBoundary.hpp>
+#include <mesh/NamedBoundaryDescriptor.hpp>
#include <mesh/StencilManager.hpp>
#include <scheme/DiscreteFunctionDPkVariant.hpp>
#include <scheme/DiscreteFunctionUtils.hpp>
#include <scheme/DiscreteFunctionVariant.hpp>
+template <size_t Dimension>
+PUGS_INLINE auto
+symmetrize_vector(const TinyVector<Dimension>& normal, const TinyVector<Dimension>& u)
+{
+ return u - 2 * dot(u, normal) * normal;
+}
+
+template <size_t Dimension>
+PUGS_INLINE auto
+symmetrize_coordinates(const TinyVector<Dimension>& origin,
+ const TinyVector<Dimension>& normal,
+ const TinyVector<Dimension>& u)
+{
+ return u - 2 * dot(u - origin, normal) * normal;
+}
+
class PolynomialReconstruction::MutableDiscreteFunctionDPkVariant
{
public:
@@ -304,6 +322,35 @@ _computeEjkMean(const TinyVector<1>& Xj,
}
}
+template <typename NodeListT>
+void
+_computeEjkMeanInSymmetricCell(const TinyVector<1>& origin,
+ const TinyVector<1>& normal,
+ const TinyVector<1>& Xj,
+ const NodeValue<const TinyVector<1>>& xr,
+ const NodeListT& node_list,
+ const CellType& cell_type,
+ const double Vi,
+ const size_t degree,
+ const size_t basis_dimension,
+ SmallArray<double>& inv_Vj_wq_detJ_ek,
+ SmallArray<double>& mean_of_ejk)
+{
+ if (cell_type == CellType::Line) {
+ const auto x0 = symmetrize_coordinates(origin, normal, xr[node_list[1]]);
+ const auto x1 = symmetrize_coordinates(origin, normal, xr[node_list[0]]);
+
+ const LineTransformation<1> T{x0, x1};
+
+ const auto& quadrature = QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor{degree});
+
+ _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
+
+ } else {
+ throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type)});
+ }
+}
+
template <typename ConformTransformationT>
void
_computeEjkBoundaryMean(const QuadratureFormula<1>& quadrature,
@@ -384,6 +431,46 @@ _computeEjkMeanByBoundary(const TinyVector<2>& Xj,
}
}
+void
+_computeEjkMeanByBoundaryInSymmetricCell(const TinyVector<2>& origin,
+ const TinyVector<2>& normal,
+ const TinyVector<2>& Xj,
+ const CellId& cell_id,
+ const NodeValue<const TinyVector<2>>& xr,
+ const auto& cell_to_face_matrix,
+ const auto& face_to_node_matrix,
+ const auto& cell_face_is_reversed,
+ const CellValue<const double>& Vi,
+ const size_t degree,
+ const size_t basis_dimension,
+ SmallArray<double>& inv_Vj_alpha_p_1_wq_X_prime_orth_ek,
+ SmallArray<double>& mean_of_ejk)
+{
+ const auto& quadrature = QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(degree + 1));
+
+ mean_of_ejk.fill(0);
+ auto cell_face_list = cell_to_face_matrix[cell_id];
+ for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
+ bool is_reversed = cell_face_is_reversed[cell_id][i_face];
+
+ const FaceId face_id = cell_face_list[i_face];
+ auto face_node_list = face_to_node_matrix[face_id];
+
+ const auto x0 = symmetrize_coordinates(origin, normal, xr[face_node_list[1]]);
+ const auto x1 = symmetrize_coordinates(origin, normal, xr[face_node_list[0]]);
+
+ if (is_reversed) {
+ const LineTransformation<2> T{x1, x0};
+ _computeEjkBoundaryMean(quadrature, T, Xj, Vi[cell_id], degree, basis_dimension,
+ inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_of_ejk);
+ } else {
+ const LineTransformation<2> T{x0, x1};
+ _computeEjkBoundaryMean(quadrature, T, Xj, Vi[cell_id], degree, basis_dimension,
+ inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_of_ejk);
+ }
+ }
+}
+
template <typename NodeListT>
void
_computeEjkMean(const TinyVector<2>& Xj,
@@ -418,6 +505,51 @@ _computeEjkMean(const TinyVector<2>& Xj,
}
}
+template <typename NodeListT>
+void
+_computeEjkMeanInSymmetricCell(const TinyVector<2>& origin,
+ const TinyVector<2>& normal,
+ const TinyVector<2>& Xj,
+ const NodeValue<const TinyVector<2>>& xr,
+ const NodeListT& node_list,
+ const CellType& cell_type,
+ const double Vi,
+ const size_t degree,
+ const size_t basis_dimension,
+ SmallArray<double>& inv_Vj_wq_detJ_ek,
+ SmallArray<double>& mean_of_ejk)
+{
+ switch (cell_type) {
+ case CellType::Triangle: {
+ const auto x0 = symmetrize_coordinates(origin, normal, xr[node_list[2]]);
+ const auto x1 = symmetrize_coordinates(origin, normal, xr[node_list[1]]);
+ const auto x2 = symmetrize_coordinates(origin, normal, xr[node_list[0]]);
+
+ const TriangleTransformation<2> T{x0, x1, x2};
+ const auto& quadrature = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor{degree});
+
+ _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
+ break;
+ }
+ case CellType::Quadrangle: {
+ const auto x0 = symmetrize_coordinates(origin, normal, xr[node_list[3]]);
+ const auto x1 = symmetrize_coordinates(origin, normal, xr[node_list[2]]);
+ const auto x2 = symmetrize_coordinates(origin, normal, xr[node_list[1]]);
+ const auto x3 = symmetrize_coordinates(origin, normal, xr[node_list[0]]);
+
+ const SquareTransformation<2> T{x0, x1, x2, x3};
+ const auto& quadrature =
+ QuadratureManager::instance().getSquareFormula(GaussLegendreQuadratureDescriptor{degree + 1});
+
+ _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
+ break;
+ }
+ default: {
+ throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type)});
+ }
+ }
+}
+
template <typename NodeListT>
void
_computeEjkMean(const TinyVector<3>& Xj,
@@ -467,6 +599,81 @@ _computeEjkMean(const TinyVector<3>& Xj,
}
}
+template <typename NodeListT>
+void
+_computeEjkMeanInSymmetricCell(const TinyVector<3>& origin,
+ const TinyVector<3>& normal,
+ const TinyVector<3>& Xj,
+ const NodeValue<const TinyVector<3>>& xr,
+ const NodeListT& node_list,
+ const CellType& cell_type,
+ const double Vi,
+ const size_t degree,
+ const size_t basis_dimension,
+ SmallArray<double>& inv_Vj_wq_detJ_ek,
+ SmallArray<double>& mean_of_ejk)
+{
+ if (cell_type == CellType::Tetrahedron) {
+ const auto x0 = symmetrize_coordinates(origin, normal, xr[node_list[1]]);
+ const auto x1 = symmetrize_coordinates(origin, normal, xr[node_list[0]]);
+ const auto x2 = symmetrize_coordinates(origin, normal, xr[node_list[2]]);
+ const auto x3 = symmetrize_coordinates(origin, normal, xr[node_list[3]]);
+
+ const TetrahedronTransformation T{x0, x1, x2, x3};
+
+ const auto& quadrature = QuadratureManager::instance().getTetrahedronFormula(GaussQuadratureDescriptor{degree});
+
+ _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
+
+ } else if (cell_type == CellType::Prism) {
+ const auto x0 = symmetrize_coordinates(origin, normal, xr[node_list[1]]);
+ const auto x1 = symmetrize_coordinates(origin, normal, xr[node_list[0]]);
+ const auto x2 = symmetrize_coordinates(origin, normal, xr[node_list[2]]);
+ const auto x3 = symmetrize_coordinates(origin, normal, xr[node_list[4]]);
+ const auto x4 = symmetrize_coordinates(origin, normal, xr[node_list[3]]);
+ const auto x5 = symmetrize_coordinates(origin, normal, xr[node_list[5]]);
+
+ const PrismTransformation T{x0, x1, x2, //
+ x3, x4, x5};
+
+ const auto& quadrature = QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor{degree + 1});
+
+ _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
+
+ } else if (cell_type == CellType::Pyramid) {
+ const auto x0 = symmetrize_coordinates(origin, normal, xr[node_list[3]]);
+ const auto x1 = symmetrize_coordinates(origin, normal, xr[node_list[2]]);
+ const auto x2 = symmetrize_coordinates(origin, normal, xr[node_list[1]]);
+ const auto x3 = symmetrize_coordinates(origin, normal, xr[node_list[0]]);
+ const auto x4 = symmetrize_coordinates(origin, normal, xr[node_list[4]]);
+ const PyramidTransformation T{x0, x1, x2, x3, x4};
+
+ const auto& quadrature = QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor{degree + 1});
+
+ _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
+
+ } else if (cell_type == CellType::Hexahedron) {
+ const auto x0 = symmetrize_coordinates(origin, normal, xr[node_list[3]]);
+ const auto x1 = symmetrize_coordinates(origin, normal, xr[node_list[2]]);
+ const auto x2 = symmetrize_coordinates(origin, normal, xr[node_list[1]]);
+ const auto x3 = symmetrize_coordinates(origin, normal, xr[node_list[0]]);
+ const auto x4 = symmetrize_coordinates(origin, normal, xr[node_list[7]]);
+ const auto x5 = symmetrize_coordinates(origin, normal, xr[node_list[6]]);
+ const auto x6 = symmetrize_coordinates(origin, normal, xr[node_list[5]]);
+ const auto x7 = symmetrize_coordinates(origin, normal, xr[node_list[4]]);
+
+ const CubeTransformation T{x0, x1, x2, x3, x4, x5, x6, x7};
+
+ const auto& quadrature =
+ QuadratureManager::instance().getCubeFormula(GaussLegendreQuadratureDescriptor{degree + 1});
+
+ _computeEjkMean(quadrature, T, Xj, Vi, degree, basis_dimension, inv_Vj_wq_detJ_ek, mean_of_ejk);
+
+ } else {
+ throw NotImplementedError("unexpected cell type: " + std::string{name(cell_type)});
+ }
+}
+
size_t
PolynomialReconstruction::_getNumberOfColumns(
const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
@@ -550,7 +757,8 @@ PolynomialReconstruction::_build(
const size_t basis_dimension =
DiscreteFunctionDPk<MeshType::Dimension, double>::BasisViewType::dimensionFromDegree(m_descriptor.degree());
- const Stencil& stencil = StencilManager::instance().getStencil(mesh.connectivity(), m_descriptor.degree());
+ const Stencil& stencil = StencilManager::instance().getStencil(mesh.connectivity(), m_descriptor.degree(),
+ m_descriptor.symmetryBoundaryNames());
auto xr = mesh.xr();
auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
@@ -561,17 +769,58 @@ PolynomialReconstruction::_build(
auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
auto cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
+ auto full_stencil_size = [&](const CellId cell_id) {
+ auto stencil_cell_list = stencil[cell_id];
+ size_t stencil_size = stencil_cell_list.size();
+ for (size_t i = 0; i < m_descriptor.symmetryBoundaryNames().size(); ++i) {
+ auto& ghost_stencil = stencil.symmetryNameStencil()[i].second;
+ stencil_size += ghost_stencil[cell_id].size();
+ }
+
+ return stencil_size;
+ };
+
const size_t max_stencil_size = [&]() {
size_t max_size = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- auto stencil_cell_list = stencil[cell_id];
- if (cell_is_owned[cell_id] and stencil_cell_list.size() > max_size) {
- max_size = stencil_cell_list.size();
+ auto stencil_cell_list = stencil[cell_id];
+ const size_t stencil_size = full_stencil_size(cell_id);
+ if (cell_is_owned[cell_id] and stencil_size > max_size) {
+ max_size = stencil_size;
}
}
return max_size;
}();
+ SmallArray<const Rd> symmetry_normal_list = [&] {
+ SmallArray<Rd> normal_list(m_descriptor.symmetryBoundaryNames().size());
+ size_t i_symmetry_boundary = 0;
+ for (auto name : m_descriptor.symmetryBoundaryNames()) {
+ auto symmetry_boundary = getMeshFlatFaceBoundary(mesh, NamedBoundaryDescriptor(name));
+ normal_list[i_symmetry_boundary++] = symmetry_boundary.outgoingNormal();
+ }
+ return normal_list;
+ }();
+ SmallArray<const Rd> symmetry_origin_list = [&] {
+ SmallArray<Rd> origin_list(m_descriptor.symmetryBoundaryNames().size());
+ size_t i_symmetry_boundary = 0;
+ for (auto name : m_descriptor.symmetryBoundaryNames()) {
+ auto symmetry_boundary = getMeshFlatFaceBoundary(mesh, NamedBoundaryDescriptor(name));
+ origin_list[i_symmetry_boundary++] = symmetry_boundary.origin();
+ }
+ return origin_list;
+ }();
+ {
+ size_t i_symmetry_boundary = 0;
+ for (auto name : m_descriptor.symmetryBoundaryNames()) {
+ auto symmetry_boundary = getMeshFlatFaceBoundary(mesh, NamedBoundaryDescriptor(name));
+ std::cout << name << " n=" << symmetry_boundary.outgoingNormal() << " o=" << symmetry_boundary.origin() << '\n';
+ ++i_symmetry_boundary;
+ }
+ }
+ std::cout << "symmetry_origin_list = " << symmetry_origin_list << '\n';
+ std::cout << "symmetry_normal_list = " << symmetry_normal_list << '\n';
+
Kokkos::Experimental::UniqueToken<Kokkos::DefaultExecutionSpace::execution_space,
Kokkos::Experimental::UniqueTokenScope::Global>
tokens;
@@ -611,7 +860,7 @@ PolynomialReconstruction::_build(
auto stencil_cell_list = stencil[cell_j_id];
- ShrinkMatrixView B(B_pool[t], stencil_cell_list.size());
+ ShrinkMatrixView B(B_pool[t], full_stencil_size(cell_j_id));
size_t column_begin = 0;
for (size_t i_discrete_function_variant = 0;
@@ -624,19 +873,58 @@ PolynomialReconstruction::_build(
if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
using DataType = std::decay_t<typename DiscreteFunctionT::data_type>;
const DataType& qj = discrete_function[cell_j_id];
- for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
+ size_t index = 0;
+ for (size_t i = 0; i < stencil_cell_list.size(); ++i, ++index) {
const CellId cell_i_id = stencil_cell_list[i];
const DataType& qi_qj = discrete_function[cell_i_id] - qj;
if constexpr (std::is_arithmetic_v<DataType>) {
- B(i, column_begin) = qi_qj;
+ B(index, column_begin) = qi_qj;
} else if constexpr (is_tiny_vector_v<DataType>) {
for (size_t kB = column_begin, k = 0; k < DataType::Dimension; ++k, ++kB) {
- B(i, kB) = qi_qj[k];
+ B(index, kB) = qi_qj[k];
}
} else if constexpr (is_tiny_matrix_v<DataType>) {
for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
- B(i, column_begin + k * DataType::NumberOfColumns + l) = qi_qj(k, l);
+ B(index, column_begin + k * DataType::NumberOfColumns + l) = qi_qj(k, l);
+ }
+ }
+ }
+ }
+
+ for (size_t i_symmetry = 0; i_symmetry < stencil.symmetryNameStencil().size(); ++i_symmetry) {
+ auto& ghost_stencil = stencil.symmetryNameStencil()[i_symmetry].second;
+ auto ghost_cell_list = ghost_stencil[cell_j_id];
+ for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
+ const CellId cell_i_id = ghost_cell_list[i];
+ if constexpr (std::is_arithmetic_v<DataType>) {
+ const DataType& qi_qj = discrete_function[cell_i_id] - qj;
+ B(index, column_begin) = qi_qj;
+ } else if constexpr (is_tiny_vector_v<DataType>) {
+ if constexpr (DataType::Dimension == MeshType::Dimension) {
+ const Rd& normal = symmetry_normal_list[i_symmetry];
+
+ const DataType& qi = discrete_function[cell_i_id];
+ const DataType& qi_qj = symmetrize_vector(normal, qi) - qj;
+ for (size_t kB = column_begin, k = 0; k < DataType::Dimension; ++k, ++kB) {
+ B(index, kB) = qi_qj[k];
+ }
+ } else {
+ const DataType& qi_qj = discrete_function[cell_i_id] - qj;
+ for (size_t kB = column_begin, k = 0; k < DataType::Dimension; ++k, ++kB) {
+ B(index, kB) = qi_qj[k];
+ }
+ }
+ } else if constexpr (is_tiny_matrix_v<DataType>) {
+ if constexpr ((DataType::NumberOfColumns == DataType::NumberOfRows) and
+ (DataType::NumberOfColumns == MeshType::Dimension)) {
+ throw NotImplementedError("TinyMatrix symmetries for reconstruction");
+ }
+ const DataType& qi_qj = discrete_function[cell_i_id] - qj;
+ for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
+ for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
+ B(index, column_begin + k * DataType::NumberOfColumns + l) = qi_qj(k, l);
+ }
}
}
}
@@ -650,12 +938,26 @@ PolynomialReconstruction::_build(
} else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
using DataType = std::decay_t<typename DiscreteFunctionT::data_type>;
const auto qj_vector = discrete_function[cell_j_id];
- for (size_t l = 0; l < qj_vector.size(); ++l) {
- const DataType& qj = qj_vector[l];
- for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
- const CellId cell_i_id = stencil_cell_list[i];
- const DataType& qi_qj = discrete_function[cell_i_id][l] - qj;
- B(i, column_begin + l) = qi_qj;
+ size_t index = 0;
+ for (size_t i = 0; i < stencil_cell_list.size(); ++i, ++index) {
+ const CellId cell_i_id = stencil_cell_list[i];
+ for (size_t l = 0; l < qj_vector.size(); ++l) {
+ const DataType& qj = qj_vector[l];
+ const DataType& qi_qj = discrete_function[cell_i_id][l] - qj;
+ B(index, column_begin + l) = qi_qj;
+ }
+ }
+
+ for (size_t i_symmetry = 0; i_symmetry < stencil.symmetryNameStencil().size(); ++i_symmetry) {
+ auto& ghost_stencil = stencil.symmetryNameStencil()[i_symmetry].second;
+ auto ghost_cell_list = ghost_stencil[cell_j_id];
+ for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
+ const CellId cell_i_id = ghost_cell_list[i];
+ for (size_t l = 0; l < qj_vector.size(); ++l) {
+ const DataType& qj = qj_vector[l];
+ const DataType& qi_qj = discrete_function[cell_i_id][l] - qj;
+ B(index, column_begin + l) = qi_qj;
+ }
}
}
column_begin += qj_vector.size();
@@ -668,16 +970,32 @@ PolynomialReconstruction::_build(
discrete_function_variant->discreteFunction());
}
- ShrinkMatrixView A(A_pool[t], stencil_cell_list.size());
+ ShrinkMatrixView A(A_pool[t], full_stencil_size(cell_j_id));
if ((m_descriptor.degree() == 1) and
(m_descriptor.integrationMethodType() == IntegrationMethodType::cell_center)) {
const Rd& Xj = xj[cell_j_id];
- for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
+ size_t index = 0;
+ for (size_t i = 0; i < stencil_cell_list.size(); ++i, ++index) {
const CellId cell_i_id = stencil_cell_list[i];
const Rd Xi_Xj = xj[cell_i_id] - Xj;
for (size_t l = 0; l < basis_dimension - 1; ++l) {
- A(i, l) = Xi_Xj[l];
+ A(index, l) = Xi_Xj[l];
+ }
+ }
+ for (size_t i_symmetry = 0; i_symmetry < stencil.symmetryNameStencil().size(); ++i_symmetry) {
+ auto& ghost_stencil = stencil.symmetryNameStencil()[i_symmetry].second;
+ auto ghost_cell_list = ghost_stencil[cell_j_id];
+
+ const Rd& origin = symmetry_origin_list[i_symmetry];
+ const Rd& normal = symmetry_normal_list[i_symmetry];
+
+ for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
+ const CellId cell_i_id = ghost_cell_list[i];
+ const Rd Xi_Xj = symmetrize_coordinates(origin, normal, xj[cell_i_id]) - Xj;
+ for (size_t l = 0; l < basis_dimension - 1; ++l) {
+ A(index, l) = Xi_Xj[l];
+ }
}
}
@@ -698,7 +1016,8 @@ PolynomialReconstruction::_build(
cell_face_is_reversed, Vj, m_descriptor.degree(), basis_dimension,
inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_j_of_ejk);
- for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
+ size_t index = 0;
+ for (size_t i = 0; i < stencil_cell_list.size(); ++i, ++index) {
const CellId cell_i_id = stencil_cell_list[i];
_computeEjkMeanByBoundary(Xj, cell_i_id, xr, cell_to_face_matrix, face_to_node_matrix,
@@ -706,7 +1025,28 @@ PolynomialReconstruction::_build(
inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_i_of_ejk);
for (size_t l = 0; l < basis_dimension - 1; ++l) {
- A(i, l) = mean_i_of_ejk[l] - mean_j_of_ejk[l];
+ A(index, l) = mean_i_of_ejk[l] - mean_j_of_ejk[l];
+ }
+ }
+ for (size_t i_symmetry = 0; i_symmetry < stencil.symmetryNameStencil().size(); ++i_symmetry) {
+ auto& ghost_stencil = stencil.symmetryNameStencil()[i_symmetry].second;
+ auto ghost_cell_list = ghost_stencil[cell_j_id];
+
+ const Rd& origin = symmetry_origin_list[i_symmetry];
+ const Rd& normal = symmetry_normal_list[i_symmetry];
+
+ for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
+ const CellId cell_i_id = ghost_cell_list[i];
+
+ _computeEjkMeanByBoundaryInSymmetricCell(origin, normal, //
+ Xj, cell_i_id, xr, cell_to_face_matrix, face_to_node_matrix,
+ cell_face_is_reversed, Vj, m_descriptor.degree(),
+ basis_dimension, inv_Vj_alpha_p_1_wq_X_prime_orth_ek,
+ mean_i_of_ejk);
+
+ for (size_t l = 0; l < basis_dimension - 1; ++l) {
+ A(index, l) = mean_i_of_ejk[l] - mean_j_of_ejk[l];
+ }
}
}
} else {
@@ -722,14 +1062,35 @@ PolynomialReconstruction::_build(
_computeEjkMean(Xj, xr, cell_to_node_matrix[cell_j_id], cell_type[cell_j_id], Vj[cell_j_id],
m_descriptor.degree(), basis_dimension, inv_Vj_wq_detJ_ek, mean_j_of_ejk);
- for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
+ size_t index = 0;
+ for (size_t i = 0; i < stencil_cell_list.size(); ++i, ++index) {
const CellId cell_i_id = stencil_cell_list[i];
_computeEjkMean(Xj, xr, cell_to_node_matrix[cell_i_id], cell_type[cell_i_id], Vj[cell_i_id],
m_descriptor.degree(), basis_dimension, inv_Vj_wq_detJ_ek, mean_i_of_ejk);
for (size_t l = 0; l < basis_dimension - 1; ++l) {
- A(i, l) = mean_i_of_ejk[l] - mean_j_of_ejk[l];
+ A(index, l) = mean_i_of_ejk[l] - mean_j_of_ejk[l];
+ }
+ }
+ for (size_t i_symmetry = 0; i_symmetry < stencil.symmetryNameStencil().size(); ++i_symmetry) {
+ auto& ghost_stencil = stencil.symmetryNameStencil()[i_symmetry].second;
+ auto ghost_cell_list = ghost_stencil[cell_j_id];
+
+ const Rd& origin = symmetry_origin_list[i_symmetry];
+ const Rd& normal = symmetry_normal_list[i_symmetry];
+
+ for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
+ const CellId cell_i_id = ghost_cell_list[i];
+
+ _computeEjkMeanInSymmetricCell(origin, normal, //
+ Xj, xr, cell_to_node_matrix[cell_i_id], cell_type[cell_i_id],
+ Vj[cell_i_id], m_descriptor.degree(), basis_dimension, inv_Vj_wq_detJ_ek,
+ mean_i_of_ejk);
+
+ for (size_t l = 0; l < basis_dimension - 1; ++l) {
+ A(index, l) = mean_i_of_ejk[l] - mean_j_of_ejk[l];
+ }
}
}
}
diff --git a/src/scheme/PolynomialReconstructionDescriptor.hpp b/src/scheme/PolynomialReconstructionDescriptor.hpp
index 1160724c4840d1457d52925b7edd55689cb3523a..9274a7282ac6347fba0166ea40644db1a75090b0 100644
--- a/src/scheme/PolynomialReconstructionDescriptor.hpp
+++ b/src/scheme/PolynomialReconstructionDescriptor.hpp
@@ -5,6 +5,7 @@
#include <utils/PugsMacros.hpp>
#include <cstddef>
+#include <set>
class PolynomialReconstructionDescriptor
{
@@ -12,6 +13,9 @@ class PolynomialReconstructionDescriptor
private:
IntegrationMethodType m_integration_method;
size_t m_degree;
+
+ std::set<std::string> m_symmetry_boundary_names;
+
bool m_preconditioning = true;
bool m_row_weighting = true;
@@ -29,6 +33,13 @@ class PolynomialReconstructionDescriptor
return m_degree;
}
+ PUGS_INLINE
+ const std::set<std::string>&
+ symmetryBoundaryNames() const
+ {
+ return m_symmetry_boundary_names;
+ }
+
PUGS_INLINE
bool
preconditioning() const
@@ -64,8 +75,11 @@ class PolynomialReconstructionDescriptor
m_row_weighting = row_weighting;
}
- PolynomialReconstructionDescriptor(const IntegrationMethodType integration_method, const size_t degree)
- : m_integration_method{integration_method}, m_degree{degree}
+#warning TEMPORARY
+ PolynomialReconstructionDescriptor(const IntegrationMethodType integration_method,
+ const size_t degree,
+ const std::set<std::string>& symmetry_boundary_names = {})
+ : m_integration_method{integration_method}, m_degree{degree}, m_symmetry_boundary_names(symmetry_boundary_names)
{}
PolynomialReconstructionDescriptor() = delete;
diff --git a/src/scheme/test_reconstruction.cpp b/src/scheme/test_reconstruction.cpp
index 01007d156bf477df78fd6b814f276ce1c5724012..a17965ac51af42672fca6b583e50acdc0a995d07 100644
--- a/src/scheme/test_reconstruction.cpp
+++ b/src/scheme/test_reconstruction.cpp
@@ -8,16 +8,11 @@ void
test_reconstruction(const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list,
uint64_t degree)
{
- std::vector descriptor_list = {PolynomialReconstructionDescriptor{IntegrationMethodType::cell_center, degree},
- PolynomialReconstructionDescriptor{IntegrationMethodType::element, degree},
- PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree}};
-
- [[maybe_unused]] auto x =
- PolynomialReconstruction{PolynomialReconstructionDescriptor{IntegrationMethodType::cell_center, degree}}.build(
- discrete_function_variant_list);
+ std::vector descriptor_list = {
+ PolynomialReconstructionDescriptor{IntegrationMethodType::boundary, degree, {"XMIN", "XMAX", "YMIN", "YMAX"}}};
for (auto&& descriptor : descriptor_list) {
- const size_t nb_loops = 50;
+ const size_t nb_loops = 1;
std::cout << "** variable list at once (" << nb_loops << " times) using " << rang::fgB::yellow
<< name(descriptor.integrationMethodType()) << rang::fg::reset << "\n";
diff --git a/tests/test_StencilBuilder.cpp b/tests/test_StencilBuilder.cpp
index 0cb54f2623541fc97c1ce5058313c4c2f3b2ab5d..54dd9755678bfe5c0128807d6a1f16e65d927e8b 100644
--- a/tests/test_StencilBuilder.cpp
+++ b/tests/test_StencilBuilder.cpp
@@ -15,101 +15,129 @@
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);
+ SECTION("inner stencil")
+ {
+ 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 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;
+ 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;
- };
+ return true;
+ };
- SECTION("1D")
- {
- SECTION("cartesian")
+ SECTION("1D")
{
- const auto& mesh = *MeshDataBaseForTests::get().cartesian1DMesh()->get<Mesh<1>>();
+ SECTION("cartesian")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().cartesian1DMesh()->get<Mesh<1>>();
- const Connectivity<1>& connectivity = mesh.connectivity();
+ const Connectivity<1>& connectivity = mesh.connectivity();
- Stencil stencil = StencilManager::instance().getStencil(connectivity);
+ Stencil stencil = StencilManager::instance().getStencil(connectivity);
- REQUIRE(is_valid(connectivity, stencil));
- }
+ REQUIRE(is_valid(connectivity, stencil));
+ }
- SECTION("unordered")
- {
- const auto& mesh = *MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
+ SECTION("unordered")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
- const Connectivity<1>& connectivity = mesh.connectivity();
- Stencil stencil = StencilManager::instance().getStencil(connectivity);
+ const Connectivity<1>& connectivity = mesh.connectivity();
+ Stencil stencil = StencilManager::instance().getStencil(connectivity);
- REQUIRE(is_valid(connectivity, stencil));
+ REQUIRE(is_valid(connectivity, stencil));
+ }
}
- }
- SECTION("2D")
- {
- SECTION("cartesian")
+ SECTION("2D")
{
- const auto& mesh = *MeshDataBaseForTests::get().cartesian2DMesh()->get<Mesh<2>>();
+ SECTION("cartesian")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().cartesian2DMesh()->get<Mesh<2>>();
- const Connectivity<2>& connectivity = mesh.connectivity();
- REQUIRE(is_valid(connectivity, StencilManager::instance().getStencil(connectivity)));
+ const Connectivity<2>& connectivity = mesh.connectivity();
+ REQUIRE(is_valid(connectivity, StencilManager::instance().getStencil(connectivity)));
+ }
+
+ SECTION("hybrid")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>();
+
+ const Connectivity<2>& connectivity = mesh.connectivity();
+ REQUIRE(is_valid(connectivity, StencilManager::instance().getStencil(connectivity)));
+ }
}
- SECTION("hybrid")
+ SECTION("3D")
{
- const auto& mesh = *MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>();
+ SECTION("carteian")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().cartesian3DMesh()->get<Mesh<3>>();
+
+ const Connectivity<3>& connectivity = mesh.connectivity();
+ REQUIRE(is_valid(connectivity, StencilManager::instance().getStencil(connectivity)));
+ }
+
+ SECTION("hybrid")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>();
- const Connectivity<2>& connectivity = mesh.connectivity();
- REQUIRE(is_valid(connectivity, StencilManager::instance().getStencil(connectivity)));
+ const Connectivity<3>& connectivity = mesh.connectivity();
+ REQUIRE(is_valid(connectivity, StencilManager::instance().getStencil(connectivity)));
+ }
}
}
- SECTION("3D")
+ SECTION("Stencil using symmetries")
{
- SECTION("carteian")
+ SECTION("2D")
{
- const auto& mesh = *MeshDataBaseForTests::get().cartesian3DMesh()->get<Mesh<3>>();
+ SECTION("cartesian")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().cartesian2DMesh()->get<Mesh<2>>();
- const Connectivity<3>& connectivity = mesh.connectivity();
- REQUIRE(is_valid(connectivity, StencilManager::instance().getStencil(connectivity)));
+ const Connectivity<2>& connectivity = mesh.connectivity();
+ StencilManager::instance().getStencil(connectivity, 1, {"XMIN", "XMAX", "YMIN", "YMAX"});
+ }
}
- SECTION("hybrid")
+ SECTION("3D")
{
- const auto& mesh = *MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>();
+ SECTION("cartesian")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().cartesian3DMesh()->get<Mesh<3>>();
- const Connectivity<3>& connectivity = mesh.connectivity();
- REQUIRE(is_valid(connectivity, StencilManager::instance().getStencil(connectivity)));
+ const Connectivity<3>& connectivity = mesh.connectivity();
+ StencilManager::instance().getStencil(connectivity, 1, {"XMIN", "XMAX", "YMIN", "YMAX", "ZMIN", "ZMAX"});
+ }
}
}
}