From 70a7e28769561b07a7258af64d48d72cfcb8f6e3 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Fri, 16 Apr 2021 17:06:11 +0200
Subject: [PATCH] Plug VectorDiamondScheme to the language
Still remains to be validated and to be cleaned-up
---
src/language/modules/SchemeModule.cpp | 21 +
src/scheme/CMakeLists.txt | 3 +-
src/scheme/VectorDiamondScheme.cpp | 977 ++++++++++++++++++++++++++
src/scheme/VectorDiamondScheme.hpp | 26 +
4 files changed, 1026 insertions(+), 1 deletion(-)
create mode 100644 src/scheme/VectorDiamondScheme.cpp
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 523ad8729..232f5b18f 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -29,6 +29,7 @@
#include <scheme/NumberedBoundaryDescriptor.hpp>
#include <scheme/ScalarDiamondScheme.hpp>
#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+#include <scheme/VectorDiamondScheme.hpp>
#include <scheme/PerfectGas.hpp>
@@ -459,6 +460,26 @@ SchemeModule::SchemeModule()
));
+ this->_addBuiltinFunction("unsteadyelasticity",
+ std::make_shared<BuiltinFunctionEmbedder<
+ std::shared_ptr<const IDiscreteFunction>(const std::shared_ptr<const IDiscreteFunction>&,
+ const std::shared_ptr<const IDiscreteFunction>&,
+ const std::shared_ptr<const IDiscreteFunction>&,
+ const std::shared_ptr<const IDiscreteFunction>&,
+ const std::vector<std::shared_ptr<
+ const IBoundaryConditionDescriptor>>&)>>(
+
+ [](const std::shared_ptr<const IDiscreteFunction> alpha,
+ const std::shared_ptr<const IDiscreteFunction> lambda,
+ const std::shared_ptr<const IDiscreteFunction> mu,
+ const std::shared_ptr<const IDiscreteFunction> f,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list) -> const std::shared_ptr<const IDiscreteFunction> {
+ return VectorDiamondSchemeHandler{alpha, lambda, mu, f, bc_descriptor_list}.solution();
+ }
+
+ ));
+
this->_addBuiltinFunction("heat2",
std::make_shared<BuiltinFunctionEmbedder<
void(std::shared_ptr<const IMesh>,
diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index 728c640ed..ab13de8d8 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -7,4 +7,5 @@ add_library(
DiscreteFunctionUtils.cpp
DiscreteFunctionVectorInterpoler.cpp
ScalarDiamondScheme.cpp
- PerfectGas.cpp)
+ PerfectGas.cpp
+ VectorDiamondScheme.cpp)
diff --git a/src/scheme/VectorDiamondScheme.cpp b/src/scheme/VectorDiamondScheme.cpp
new file mode 100644
index 000000000..92a262cb7
--- /dev/null
+++ b/src/scheme/VectorDiamondScheme.cpp
@@ -0,0 +1,977 @@
+#include <scheme/VectorDiamondScheme.hpp>
+
+#include <scheme/DiscreteFunctionP0.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
+
+template <size_t Dimension>
+class VectorDiamondSchemeHandler::InterpolationWeightsManager
+{
+ private:
+ std::shared_ptr<const Mesh<Connectivity<Dimension>>> m_mesh;
+ FaceValue<const bool> m_primal_face_is_on_boundary;
+ NodeValue<const bool> m_primal_node_is_on_boundary;
+ FaceValue<const bool> m_primal_face_is_symmetry;
+ CellValuePerNode<double> m_w_rj;
+ FaceValuePerNode<double> m_w_rl;
+
+ public:
+ InterpolationWeightsManager(std::shared_ptr<const Mesh<Connectivity<Dimension>>> mesh,
+ FaceValue<const bool> primal_face_is_on_boundary,
+ NodeValue<const bool> primal_node_is_on_boundary,
+ FaceValue<const bool> primal_face_is_symmetry)
+ : m_mesh(mesh),
+ m_primal_face_is_on_boundary(primal_face_is_on_boundary),
+ m_primal_node_is_on_boundary(primal_node_is_on_boundary),
+ m_primal_face_is_symmetry(primal_face_is_symmetry)
+ {}
+ ~InterpolationWeightsManager() = default;
+ CellValuePerNode<double>&
+ wrj()
+ {
+ return m_w_rj;
+ }
+ FaceValuePerNode<double>&
+ wrl()
+ {
+ return m_w_rl;
+ }
+ void
+ compute()
+ {
+ using MeshDataType = MeshData<Dimension>;
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
+
+ const NodeValue<const TinyVector<Dimension>>& xr = m_mesh->xr();
+
+ const FaceValue<const TinyVector<Dimension>>& xl = mesh_data.xl();
+ const CellValue<const TinyVector<Dimension>>& xj = mesh_data.xj();
+ const auto& node_to_cell_matrix = m_mesh->connectivity().nodeToCellMatrix();
+ const auto& node_to_face_matrix = m_mesh->connectivity().nodeToFaceMatrix();
+ const auto& face_to_cell_matrix = m_mesh->connectivity().faceToCellMatrix();
+
+ CellValuePerNode<double> w_rj{m_mesh->connectivity()};
+ FaceValuePerNode<double> w_rl{m_mesh->connectivity()};
+
+ const NodeValuePerFace<const TinyVector<Dimension>> primal_nlr = mesh_data.nlr();
+ auto project_to_face = [&](const TinyVector<Dimension>& x, const FaceId face_id) -> const TinyVector<Dimension> {
+ TinyVector<Dimension> proj;
+ const TinyVector<Dimension> nil = primal_nlr(face_id, 0);
+ proj = x - ((x - xl[face_id]), nil) * nil;
+ return proj;
+ };
+
+ for (size_t i = 0; i < w_rl.numberOfValues(); ++i) {
+ w_rl[i] = std::numeric_limits<double>::signaling_NaN();
+ }
+
+ for (NodeId i_node = 0; i_node < m_mesh->numberOfNodes(); ++i_node) {
+ Vector<double> b{Dimension + 1};
+ b[0] = 1;
+ for (size_t i = 1; i < Dimension + 1; i++) {
+ b[i] = xr[i_node][i - 1];
+ }
+ const auto& node_to_cell = node_to_cell_matrix[i_node];
+
+ if (not m_primal_node_is_on_boundary[i_node]) {
+ LocalRectangularMatrix<double> A{Dimension + 1, node_to_cell.size()};
+ for (size_t j = 0; j < node_to_cell.size(); j++) {
+ A(0, j) = 1;
+ }
+ for (size_t i = 1; i < Dimension + 1; i++) {
+ for (size_t j = 0; j < node_to_cell.size(); j++) {
+ const CellId J = node_to_cell[j];
+ A(i, j) = xj[J][i - 1];
+ }
+ }
+
+ Vector<double> x{node_to_cell.size()};
+ x = 0;
+
+ LeastSquareSolver ls_solver;
+ ls_solver.solveLocalSystem(A, x, b);
+
+ for (size_t j = 0; j < node_to_cell.size(); j++) {
+ w_rj(i_node, j) = x[j];
+ }
+ } else {
+ int nb_face_used = 0;
+ for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) {
+ FaceId face_id = node_to_face_matrix[i_node][i_face];
+ if (m_primal_face_is_on_boundary[face_id]) {
+ nb_face_used++;
+ }
+ }
+ LocalRectangularMatrix<double> A{Dimension + 1, node_to_cell.size() + nb_face_used};
+ for (size_t j = 0; j < node_to_cell.size() + nb_face_used; j++) {
+ A(0, j) = 1;
+ }
+ for (size_t i = 1; i < Dimension + 1; i++) {
+ for (size_t j = 0; j < node_to_cell.size(); j++) {
+ const CellId J = node_to_cell[j];
+ A(i, j) = xj[J][i - 1];
+ }
+ }
+ for (size_t i = 1; i < Dimension + 1; i++) {
+ int cpt_face = 0;
+ for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) {
+ FaceId face_id = node_to_face_matrix[i_node][i_face];
+ if (m_primal_face_is_on_boundary[face_id]) {
+ if (m_primal_face_is_symmetry[face_id]) {
+ for (size_t j = 0; j < face_to_cell_matrix[face_id].size(); ++j) {
+ const CellId cell_id = face_to_cell_matrix[face_id][j];
+ TinyVector<Dimension> xproj = project_to_face(xj[cell_id], face_id);
+ A(i, node_to_cell.size() + cpt_face) = xproj[i - 1];
+ }
+ } else {
+ A(i, node_to_cell.size() + cpt_face) = xl[face_id][i - 1];
+ }
+ cpt_face++;
+ }
+ }
+ }
+
+ Vector<double> x{node_to_cell.size() + nb_face_used};
+ x = 0;
+
+ LeastSquareSolver ls_solver;
+ ls_solver.solveLocalSystem(A, x, b);
+
+ for (size_t j = 0; j < node_to_cell.size(); j++) {
+ w_rj(i_node, j) = x[j];
+ }
+ int cpt_face = node_to_cell.size();
+ for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) {
+ FaceId face_id = node_to_face_matrix[i_node][i_face];
+ if (m_primal_face_is_on_boundary[face_id]) {
+ w_rl(i_node, i_face) = x[cpt_face++];
+ }
+ }
+ }
+ }
+ m_w_rj = w_rj;
+ m_w_rl = w_rl;
+ }
+};
+
+class VectorDiamondSchemeHandler::IVectorDiamondScheme
+{
+ public:
+ virtual std::shared_ptr<const IDiscreteFunction> getSolution() const = 0;
+
+ IVectorDiamondScheme() = default;
+ virtual ~IVectorDiamondScheme() = default;
+};
+
+template <size_t Dimension>
+class VectorDiamondSchemeHandler::VectorDiamondScheme : public VectorDiamondSchemeHandler::IVectorDiamondScheme
+{
+ private:
+ using ConnectivityType = Connectivity<Dimension>;
+ using MeshType = Mesh<ConnectivityType>;
+ using MeshDataType = MeshData<Dimension>;
+
+ std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>> m_solution;
+
+ class DirichletBoundaryCondition
+ {
+ private:
+ const Array<const TinyVector<Dimension>> m_value_list;
+ const Array<const FaceId> m_face_list;
+
+ public:
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_face_list;
+ }
+
+ const Array<const TinyVector<Dimension>>&
+ valueList() const
+ {
+ return m_value_list;
+ }
+
+ DirichletBoundaryCondition(const Array<const FaceId>& face_list,
+ const Array<const TinyVector<Dimension>>& value_list)
+ : m_value_list{value_list}, m_face_list{face_list}
+ {
+ Assert(m_value_list.size() == m_face_list.size());
+ }
+
+ ~DirichletBoundaryCondition() = default;
+ };
+
+ class NormalStrainBoundaryCondition
+ {
+ private:
+ const Array<const TinyVector<Dimension>> m_value_list;
+ const Array<const FaceId> m_face_list;
+
+ public:
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_face_list;
+ }
+
+ const Array<const TinyVector<Dimension>>&
+ valueList() const
+ {
+ return m_value_list;
+ }
+
+ NormalStrainBoundaryCondition(const Array<const FaceId>& face_list,
+ const Array<const TinyVector<Dimension>>& value_list)
+ : m_value_list{value_list}, m_face_list{face_list}
+ {
+ Assert(m_value_list.size() == m_face_list.size());
+ }
+
+ ~NormalStrainBoundaryCondition() = default;
+ };
+
+ class SymmetryBoundaryCondition
+ {
+ private:
+ const Array<const TinyVector<Dimension>> m_value_list;
+ const Array<const FaceId> m_face_list;
+
+ public:
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_face_list;
+ }
+
+ public:
+ SymmetryBoundaryCondition(const Array<const FaceId>& face_list) : m_face_list{face_list} {}
+
+ ~SymmetryBoundaryCondition() = default;
+ };
+
+ public:
+ std::shared_ptr<const IDiscreteFunction>
+ getSolution() const final
+ {
+ return m_solution;
+ }
+
+ VectorDiamondScheme(const std::shared_ptr<const MeshType>& mesh,
+ const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& alpha,
+ const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& dual_lambda,
+ const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& dual_mu,
+ const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>>& f,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
+ {
+ Assert(mesh == alpha->mesh());
+ Assert(mesh == f->mesh());
+ Assert(dual_lambda->mesh() == dual_mu->mesh());
+ if (DiamondDualMeshManager::instance().getDiamondDualMesh(mesh) != dual_mu->mesh()) {
+ throw NormalError("diffusion coefficient is not defined on the dual mesh!");
+ }
+
+ // using ConnectivityType = Connectivity<Dimension>;
+ // using MeshType = Mesh<ConnectivityType>;
+ using MeshDataType = MeshData<Dimension>;
+
+ constexpr ItemType FaceType = [] {
+ if constexpr (Dimension > 1) {
+ return ItemType::face;
+ } else {
+ return ItemType::node;
+ }
+ }();
+
+ using BoundaryCondition =
+ std::variant<DirichletBoundaryCondition, NormalStrainBoundaryCondition, SymmetryBoundaryCondition>;
+
+ using BoundaryConditionList = std::vector<BoundaryCondition>;
+
+ BoundaryConditionList boundary_condition_list;
+
+ NodeValue<bool> is_dirichlet{mesh->connectivity()};
+ is_dirichlet.fill(false);
+ NodeValue<TinyVector<Dimension>> dirichlet_value{mesh->connectivity()};
+ {
+ TinyVector<Dimension> nan_tiny_vector;
+ for (size_t i = 0; i < Dimension; ++i) {
+ nan_tiny_vector[i] = std::numeric_limits<double>::signaling_NaN();
+ }
+ dirichlet_value.fill(nan_tiny_vector);
+ }
+
+ for (const auto& bc_descriptor : bc_descriptor_list) {
+ bool is_valid_boundary_condition = true;
+
+ switch (bc_descriptor->type()) {
+ case IBoundaryConditionDescriptor::Type::symmetry: {
+ const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
+ dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
+ for (size_t i_ref_face_list = 0;
+ i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
+ const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
+ const RefId& ref = ref_face_list.refId();
+ if (ref == sym_bc_descriptor.boundaryDescriptor()) {
+ if constexpr (Dimension > 1) {
+ Array<const FaceId> face_list = ref_face_list.list();
+ boundary_condition_list.push_back(SymmetryBoundaryCondition{face_list});
+ } else {
+ throw NotImplementedError("Symmetry conditions are not supported in 1d");
+ }
+ }
+ }
+ break;
+ }
+ case IBoundaryConditionDescriptor::Type::dirichlet: {
+ const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
+ dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
+ if (dirichlet_bc_descriptor.name() == "dirichlet") {
+ for (size_t i_ref_face_list = 0;
+ i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
+ const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
+ const RefId& ref = ref_face_list.refId();
+ if (ref == dirichlet_bc_descriptor.boundaryDescriptor()) {
+ const FunctionSymbolId g_id = dirichlet_bc_descriptor.rhsSymbolId();
+
+ if constexpr (Dimension > 1) {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+
+ Array<const FaceId> face_list = ref_face_list.list();
+ Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
+ TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id, mesh_data.xl(), face_list);
+ boundary_condition_list.push_back(DirichletBoundaryCondition{face_list, value_list});
+ } else {
+ throw NotImplementedError("Neumann conditions are not supported in 1d");
+ }
+ }
+ }
+ } else if (dirichlet_bc_descriptor.name() == "normal_strain") {
+ for (size_t i_ref_face_list = 0;
+ i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
+ const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
+ const RefId& ref = ref_face_list.refId();
+ if (ref == dirichlet_bc_descriptor.boundaryDescriptor()) {
+ const FunctionSymbolId g_id = dirichlet_bc_descriptor.rhsSymbolId();
+
+ if constexpr (Dimension > 1) {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+
+ Array<const FaceId> face_list = ref_face_list.list();
+ Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
+ TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id, mesh_data.xl(), face_list);
+ boundary_condition_list.push_back(NormalStrainBoundaryCondition{face_list, value_list});
+ } else {
+ throw NotImplementedError("Neumann conditions are not supported in 1d");
+ }
+ }
+ }
+ } else {
+ is_valid_boundary_condition = false;
+ }
+ break;
+ }
+ default: {
+ is_valid_boundary_condition = false;
+ }
+ }
+ if (not is_valid_boundary_condition) {
+ std::ostringstream error_msg;
+ error_msg << *bc_descriptor << " is an invalid boundary condition for elasticity equation";
+ throw NormalError(error_msg.str());
+ }
+ }
+
+ if constexpr (Dimension > 1) {
+ const CellValue<const size_t> cell_dof_number = [&] {
+ CellValue<size_t> compute_cell_dof_number{mesh->connectivity()};
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { compute_cell_dof_number[cell_id] = cell_id; });
+ return compute_cell_dof_number;
+ }();
+ size_t number_of_dof = mesh->numberOfCells();
+
+ const FaceValue<const size_t> face_dof_number = [&] {
+ FaceValue<size_t> compute_face_dof_number{mesh->connectivity()};
+ compute_face_dof_number.fill(std::numeric_limits<size_t>::max());
+ for (const auto& boundary_condition : boundary_condition_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr ((std::is_same_v<T, NormalStrainBoundaryCondition>) or
+ (std::is_same_v<T, SymmetryBoundaryCondition>) or
+ (std::is_same_v<T, DirichletBoundaryCondition>)) {
+ const auto& face_list = bc.faceList();
+
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+ if (compute_face_dof_number[face_id] != std::numeric_limits<size_t>::max()) {
+ std::ostringstream os;
+ os << "The face " << face_id << " is used at least twice for boundary conditions";
+ throw NormalError(os.str());
+ } else {
+ compute_face_dof_number[face_id] = number_of_dof++;
+ }
+ }
+ }
+ },
+ boundary_condition);
+ }
+
+ return compute_face_dof_number;
+ }();
+
+ const auto& primal_face_to_node_matrix = mesh->connectivity().faceToNodeMatrix();
+ const auto& face_to_cell_matrix = mesh->connectivity().faceToCellMatrix();
+ const FaceValue<const bool> primal_face_is_neumann = [&] {
+ FaceValue<bool> face_is_neumann{mesh->connectivity()};
+ face_is_neumann.fill(false);
+ for (const auto& boundary_condition : boundary_condition_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr ((std::is_same_v<T, NormalStrainBoundaryCondition>)) {
+ const auto& face_list = bc.faceList();
+
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+ face_is_neumann[face_id] = true;
+ }
+ }
+ },
+ boundary_condition);
+ }
+
+ return face_is_neumann;
+ }();
+
+ const FaceValue<const bool> primal_face_is_symmetry = [&] {
+ FaceValue<bool> face_is_symmetry{mesh->connectivity()};
+ face_is_symmetry.fill(false);
+ for (const auto& boundary_condition : boundary_condition_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr ((std::is_same_v<T, SymmetryBoundaryCondition>)) {
+ const auto& face_list = bc.faceList();
+
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+ face_is_symmetry[face_id] = true;
+ }
+ }
+ },
+ boundary_condition);
+ }
+
+ return face_is_symmetry;
+ }();
+
+ NodeValue<bool> primal_node_is_on_boundary(mesh->connectivity());
+ if (parallel::size() > 1) {
+ throw NotImplementedError("Calculation of node_is_on_boundary is incorrect");
+ }
+
+ primal_node_is_on_boundary.fill(false);
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ if (face_to_cell_matrix[face_id].size() == 1) {
+ for (size_t i_node = 0; i_node < primal_face_to_node_matrix[face_id].size(); ++i_node) {
+ NodeId node_id = primal_face_to_node_matrix[face_id][i_node];
+ primal_node_is_on_boundary[node_id] = true;
+ }
+ }
+ }
+
+ FaceValue<bool> primal_face_is_on_boundary(mesh->connectivity());
+ if (parallel::size() > 1) {
+ throw NotImplementedError("Calculation of face_is_on_boundary is incorrect");
+ }
+
+ primal_face_is_on_boundary.fill(false);
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ if (face_to_cell_matrix[face_id].size() == 1) {
+ primal_face_is_on_boundary[face_id] = true;
+ }
+ }
+
+ FaceValue<bool> primal_face_is_dirichlet(mesh->connectivity());
+ if (parallel::size() > 1) {
+ throw NotImplementedError("Calculation of face_is_neumann is incorrect");
+ }
+
+ primal_face_is_dirichlet.fill(false);
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ primal_face_is_dirichlet[face_id] = (primal_face_is_on_boundary[face_id] &&
+ (!primal_face_is_neumann[face_id]) && (!primal_face_is_symmetry[face_id]));
+ }
+
+ InterpolationWeightsManager iwm(mesh, primal_face_is_on_boundary, primal_node_is_on_boundary,
+ primal_face_is_symmetry);
+ iwm.compute();
+ CellValuePerNode<double> w_rj = iwm.wrj();
+ FaceValuePerNode<double> w_rl = iwm.wrl();
+
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+
+ const FaceValue<const TinyVector<Dimension>>& xl = mesh_data.xl();
+ const CellValue<const TinyVector<Dimension>>& xj = mesh_data.xj();
+ // const auto& node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
+ const auto& node_to_face_matrix = mesh->connectivity().nodeToFaceMatrix();
+ const NodeValuePerFace<const TinyVector<Dimension>> primal_nlr = mesh_data.nlr();
+
+ {
+ std::shared_ptr diamond_mesh = DiamondDualMeshManager::instance().getDiamondDualMesh(mesh);
+
+ MeshDataType& diamond_mesh_data = MeshDataManager::instance().getMeshData(*diamond_mesh);
+
+ std::shared_ptr mapper =
+ DiamondDualConnectivityManager::instance().getConnectivityToDiamondDualConnectivityDataMapper(
+ mesh->connectivity());
+
+ CellValue<const double> dual_muj = dual_mu->cellValues();
+ CellValue<const double> dual_lambdaj = dual_lambda->cellValues();
+
+ CellValue<const TinyVector<Dimension>> fj = f->cellValues();
+
+ const CellValue<const double> dual_Vj = diamond_mesh_data.Vj();
+
+ const FaceValue<const double> mes_l = [&] {
+ if constexpr (Dimension == 1) {
+ FaceValue<double> compute_mes_l{mesh->connectivity()};
+ compute_mes_l.fill(1);
+ return compute_mes_l;
+ } else {
+ return mesh_data.ll();
+ }
+ }();
+
+ const CellValue<const double> dual_mes_l_j = [=] {
+ CellValue<double> compute_mes_j{diamond_mesh->connectivity()};
+ mapper->toDualCell(mes_l, compute_mes_j);
+
+ return compute_mes_j;
+ }();
+
+ const CellValue<const double> primal_Vj = mesh_data.Vj();
+ FaceValue<const CellId> face_dual_cell_id = [=]() {
+ FaceValue<CellId> computed_face_dual_cell_id{mesh->connectivity()};
+ CellValue<CellId> dual_cell_id{diamond_mesh->connectivity()};
+ parallel_for(
+ diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { dual_cell_id[cell_id] = cell_id; });
+
+ mapper->fromDualCell(dual_cell_id, computed_face_dual_cell_id);
+
+ return computed_face_dual_cell_id;
+ }();
+
+ NodeValue<const NodeId> dual_node_primal_node_id = [=]() {
+ CellValue<NodeId> cell_ignored_id{mesh->connectivity()};
+ cell_ignored_id.fill(NodeId{std::numeric_limits<unsigned int>::max()});
+
+ NodeValue<NodeId> node_primal_id{mesh->connectivity()};
+
+ parallel_for(
+ mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { node_primal_id[node_id] = node_id; });
+
+ NodeValue<NodeId> computed_dual_node_primal_node_id{diamond_mesh->connectivity()};
+
+ mapper->toDualNode(node_primal_id, cell_ignored_id, computed_dual_node_primal_node_id);
+
+ return computed_dual_node_primal_node_id;
+ }();
+
+ CellValue<NodeId> primal_cell_dual_node_id = [=]() {
+ CellValue<NodeId> cell_id{mesh->connectivity()};
+ NodeValue<NodeId> node_ignored_id{mesh->connectivity()};
+ node_ignored_id.fill(NodeId{std::numeric_limits<unsigned int>::max()});
+
+ NodeValue<NodeId> dual_node_id{diamond_mesh->connectivity()};
+
+ parallel_for(
+ diamond_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { dual_node_id[node_id] = node_id; });
+
+ CellValue<NodeId> computed_primal_cell_dual_node_id{mesh->connectivity()};
+
+ mapper->fromDualNode(dual_node_id, node_ignored_id, cell_id);
+
+ return cell_id;
+ }();
+ const auto& dual_Cjr = diamond_mesh_data.Cjr();
+ FaceValue<TinyVector<Dimension>> dualClj = [&] {
+ FaceValue<TinyVector<Dimension>> computedClj{mesh->connectivity()};
+ const auto& dual_node_to_cell_matrix = diamond_mesh->connectivity().nodeToCellMatrix();
+ const auto& dual_cell_to_node_matrix = diamond_mesh->connectivity().cellToNodeMatrix();
+ parallel_for(
+ mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
+ const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
+ for (size_t i = 0; i < primal_face_to_cell.size(); i++) {
+ CellId cell_id = primal_face_to_cell[i];
+ const NodeId dual_node_id = primal_cell_dual_node_id[cell_id];
+ for (size_t i_dual_cell = 0; i_dual_cell < dual_node_to_cell_matrix[dual_node_id].size();
+ i_dual_cell++) {
+ const CellId dual_cell_id = dual_node_to_cell_matrix[dual_node_id][i_dual_cell];
+ if (face_dual_cell_id[face_id] == dual_cell_id) {
+ for (size_t i_dual_node = 0; i_dual_node < dual_cell_to_node_matrix[dual_cell_id].size();
+ i_dual_node++) {
+ const NodeId final_dual_node_id = dual_cell_to_node_matrix[dual_cell_id][i_dual_node];
+ if (final_dual_node_id == dual_node_id) {
+ computedClj[face_id] = dual_Cjr(dual_cell_id, i_dual_node);
+ }
+ }
+ }
+ }
+ }
+ });
+ return computedClj;
+ }();
+
+ FaceValue<TinyVector<Dimension>> nlj = [&] {
+ FaceValue<TinyVector<Dimension>> computedNlj{mesh->connectivity()};
+ parallel_for(
+ mesh->numberOfFaces(),
+ PUGS_LAMBDA(FaceId face_id) { computedNlj[face_id] = 1. / l2Norm(dualClj[face_id]) * dualClj[face_id]; });
+ return computedNlj;
+ }();
+
+ FaceValue<const double> alpha_lambda_l = [&] {
+ CellValue<double> alpha_j{diamond_mesh->connectivity()};
+
+ parallel_for(
+ diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId diamond_cell_id) {
+ alpha_j[diamond_cell_id] = dual_lambdaj[diamond_cell_id] / dual_Vj[diamond_cell_id];
+ });
+
+ FaceValue<double> computed_alpha_l{mesh->connectivity()};
+ mapper->fromDualCell(alpha_j, computed_alpha_l);
+ return computed_alpha_l;
+ }();
+
+ FaceValue<const double> alpha_mu_l = [&] {
+ CellValue<double> alpha_j{diamond_mesh->connectivity()};
+
+ parallel_for(
+ diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId diamond_cell_id) {
+ alpha_j[diamond_cell_id] = dual_muj[diamond_cell_id] / dual_Vj[diamond_cell_id];
+ });
+
+ FaceValue<double> computed_alpha_l{mesh->connectivity()};
+ mapper->fromDualCell(alpha_j, computed_alpha_l);
+ return computed_alpha_l;
+ }();
+
+ TinyMatrix<Dimension, double> eye = identity;
+
+ SparseMatrixDescriptor S{number_of_dof * Dimension};
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ const double beta_mu_l = l2Norm(dualClj[face_id]) * alpha_mu_l[face_id] * mes_l[face_id];
+ const double beta_lambda_l = l2Norm(dualClj[face_id]) * alpha_lambda_l[face_id] * mes_l[face_id];
+ const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
+ for (size_t i_cell = 0; i_cell < primal_face_to_cell.size(); ++i_cell) {
+ const CellId i_id = primal_face_to_cell[i_cell];
+ const bool is_face_reversed_for_cell_i = ((dualClj[face_id], xl[face_id] - xj[i_id]) < 0);
+
+ const TinyVector<Dimension> nil = [&] {
+ if (is_face_reversed_for_cell_i) {
+ return -nlj[face_id];
+ } else {
+ return nlj[face_id];
+ }
+ }();
+ for (size_t j_cell = 0; j_cell < primal_face_to_cell.size(); ++j_cell) {
+ const CellId j_id = primal_face_to_cell[j_cell];
+ TinyMatrix<Dimension, double> M =
+ beta_mu_l * eye + beta_mu_l * tensorProduct(nil, nil) + beta_lambda_l * tensorProduct(nil, nil);
+ TinyMatrix<Dimension, double> N = tensorProduct(nil, nil);
+
+ if (i_cell == j_cell) {
+ for (size_t i = 0; i < Dimension; ++i) {
+ for (size_t j = 0; j < Dimension; ++j) {
+ S((cell_dof_number[i_id] * Dimension) + i, (cell_dof_number[j_id] * Dimension) + j) += M(i, j);
+ if (primal_face_is_neumann[face_id]) {
+ S(face_dof_number[face_id] * Dimension + i, cell_dof_number[j_id] * Dimension + j) -= M(i, j);
+ }
+ if (primal_face_is_symmetry[face_id]) {
+ S(face_dof_number[face_id] * Dimension + i, cell_dof_number[j_id] * Dimension + j) +=
+ -((i == j) ? 1 : 0) + N(i, j);
+ S(face_dof_number[face_id] * Dimension + i, face_dof_number[face_id] * Dimension + j) +=
+ (i == j) ? 1 : 0;
+ }
+ }
+ }
+ } else {
+ for (size_t i = 0; i < Dimension; ++i) {
+ for (size_t j = 0; j < Dimension; ++j) {
+ S((cell_dof_number[i_id] * Dimension) + i, (cell_dof_number[j_id] * Dimension) + j) -= M(i, j);
+ }
+ }
+ }
+ }
+ }
+ }
+
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ for (size_t i = 0; i < Dimension; ++i) {
+ const size_t j = cell_dof_number[cell_id] * Dimension + i;
+ S(j, j) += (*alpha)[cell_id] * primal_Vj[cell_id];
+ }
+ }
+
+ const auto& dual_cell_to_node_matrix = diamond_mesh->connectivity().cellToNodeMatrix();
+ const auto& primal_node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ const double alpha_mu_face_id = mes_l[face_id] * alpha_mu_l[face_id];
+ const double alpha_lambda_face_id = mes_l[face_id] * alpha_lambda_l[face_id];
+
+ for (size_t i_face_cell = 0; i_face_cell < face_to_cell_matrix[face_id].size(); ++i_face_cell) {
+ CellId i_id = face_to_cell_matrix[face_id][i_face_cell];
+ const bool is_face_reversed_for_cell_i = ((dualClj[face_id], xl[face_id] - xj[i_id]) < 0);
+
+ for (size_t i_node = 0; i_node < primal_face_to_node_matrix[face_id].size(); ++i_node) {
+ NodeId node_id = primal_face_to_node_matrix[face_id][i_node];
+
+ const TinyVector<Dimension> nil = [&] {
+ if (is_face_reversed_for_cell_i) {
+ return -nlj[face_id];
+ } else {
+ return nlj[face_id];
+ }
+ }();
+
+ CellId dual_cell_id = face_dual_cell_id[face_id];
+
+ for (size_t i_dual_node = 0; i_dual_node < dual_cell_to_node_matrix[dual_cell_id].size(); ++i_dual_node) {
+ const NodeId dual_node_id = dual_cell_to_node_matrix[dual_cell_id][i_dual_node];
+ if (dual_node_primal_node_id[dual_node_id] == node_id) {
+ const TinyVector<Dimension> Clr = dual_Cjr(dual_cell_id, i_dual_node);
+
+ TinyMatrix<Dimension, double> M = alpha_mu_face_id * (Clr, nil) * eye +
+ alpha_mu_face_id * tensorProduct(Clr, nil) +
+ alpha_lambda_face_id * tensorProduct(nil, Clr);
+
+ for (size_t j_cell = 0; j_cell < primal_node_to_cell_matrix[node_id].size(); ++j_cell) {
+ CellId j_id = primal_node_to_cell_matrix[node_id][j_cell];
+ for (size_t i = 0; i < Dimension; ++i) {
+ for (size_t j = 0; j < Dimension; ++j) {
+ S((cell_dof_number[i_id] * Dimension) + i, (cell_dof_number[j_id] * Dimension) + j) -=
+ w_rj(node_id, j_cell) * M(i, j);
+ if (primal_face_is_neumann[face_id]) {
+ S(face_dof_number[face_id] * Dimension + i, cell_dof_number[j_id] * Dimension + j) +=
+ w_rj(node_id, j_cell) * M(i, j);
+ }
+ }
+ }
+ }
+ if (primal_node_is_on_boundary[node_id]) {
+ for (size_t l_face = 0; l_face < node_to_face_matrix[node_id].size(); ++l_face) {
+ FaceId l_id = node_to_face_matrix[node_id][l_face];
+ if (primal_face_is_on_boundary[l_id]) {
+ for (size_t i = 0; i < Dimension; ++i) {
+ for (size_t j = 0; j < Dimension; ++j) {
+ S(cell_dof_number[i_id] * Dimension + i, face_dof_number[l_id] * Dimension + j) -=
+ w_rl(node_id, l_face) * M(i, j);
+ }
+ }
+ if (primal_face_is_neumann[face_id]) {
+ for (size_t i = 0; i < Dimension; ++i) {
+ for (size_t j = 0; j < Dimension; ++j) {
+ S(face_dof_number[face_id] * Dimension + i, face_dof_number[l_id] * Dimension + j) +=
+ w_rl(node_id, l_face) * M(i, j);
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ // }
+ }
+ }
+ }
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ if (primal_face_is_dirichlet[face_id]) {
+ for (size_t i = 0; i < Dimension; ++i) {
+ S(face_dof_number[face_id] * Dimension + i, face_dof_number[face_id] * Dimension + i) += 1;
+ }
+ }
+ }
+
+ Vector<double> b{number_of_dof * Dimension};
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ for (size_t i = 0; i < Dimension; ++i) {
+ b[(cell_dof_number[cell_id] * Dimension) + i] = primal_Vj[cell_id] * fj[cell_id][i];
+ }
+ }
+
+ // Dirichlet
+ NodeValue<bool> node_tag{mesh->connectivity()};
+ node_tag.fill(false);
+ for (const auto& boundary_condition : boundary_condition_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr (std::is_same_v<T, DirichletBoundaryCondition>) {
+ const auto& face_list = bc.faceList();
+ const auto& value_list = bc.valueList();
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+
+ for (size_t i = 0; i < Dimension; ++i) {
+ b[(face_dof_number[face_id] * Dimension) + i] += value_list[i_face][i];
+ }
+ }
+ }
+ },
+ boundary_condition);
+ }
+
+ for (const auto& boundary_condition : boundary_condition_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr ((std::is_same_v<T, NormalStrainBoundaryCondition>)) {
+ const auto& face_list = bc.faceList();
+ const auto& value_list = bc.valueList();
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ FaceId face_id = face_list[i_face];
+ for (size_t i = 0; i < Dimension; ++i) {
+ b[face_dof_number[face_id] * Dimension + i] += mes_l[face_id] * value_list[i_face][i]; // sign
+ }
+ }
+ }
+ },
+ boundary_condition);
+ }
+
+ CRSMatrix A{S};
+ Vector<double> U{number_of_dof * Dimension};
+ U = 1;
+ Vector r = A * U - b;
+ std::cout << "initial (real) residu = " << std::sqrt((r, r)) << '\n';
+
+ LinearSolver solver;
+ solver.solveLocalSystem(A, U, b);
+
+ r = A * U - b;
+
+ std::cout << "final (real) residu = " << std::sqrt((r, r)) << '\n';
+
+ m_solution = std::make_shared<DiscreteFunctionP0<Dimension, TinyVector<Dimension>>>(mesh);
+ auto& solution = *m_solution;
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+ for (size_t i = 0; i < Dimension; ++i) {
+ solution[cell_id][i] = U[(cell_dof_number[cell_id] * Dimension) + i];
+ }
+ });
+
+ // NodeValue<TinyVector<3>> ur3d{mesh->connectivity()};
+ // ur3d.fill(zero);
+
+ // parallel_for(
+ // mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ // TinyVector<Dimension> x = zero;
+ // const auto node_cells = node_to_cell_matrix[node_id];
+ // for (size_t i_cell = 0; i_cell < node_cells.size(); ++i_cell) {
+ // CellId cell_id = node_cells[i_cell];
+ // x += w_rj(node_id, i_cell) * Uj[cell_id];
+ // }
+ // const auto node_faces = node_to_face_matrix[node_id];
+ // for (size_t i_face = 0; i_face < node_faces.size(); ++i_face) {
+ // FaceId face_id = node_faces[i_face];
+ // if (primal_face_is_on_boundary[face_id]) {
+ // x += w_rl(node_id, i_face) * Ul[face_id];
+ // }
+ // }
+ // for (size_t i = 0; i < Dimension; ++i) {
+ // ur3d[node_id][i] = x[i];
+ // }
+ // });
+ }
+ } else {
+ throw NotImplementedError("not done in 1d");
+ }
+ }
+};
+
+std::shared_ptr<const IDiscreteFunction>
+VectorDiamondSchemeHandler::solution() const
+{
+ return m_scheme->getSolution();
+}
+
+VectorDiamondSchemeHandler::VectorDiamondSchemeHandler(
+ const std::shared_ptr<const IDiscreteFunction>& alpha,
+ const std::shared_ptr<const IDiscreteFunction>& dual_lambda,
+ const std::shared_ptr<const IDiscreteFunction>& dual_mu,
+ const std::shared_ptr<const IDiscreteFunction>& f,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
+{
+ const std::shared_ptr i_mesh = getCommonMesh({alpha, f});
+ const std::shared_ptr i_dual_mesh = getCommonMesh({dual_lambda, dual_mu});
+ checkDiscretizationType({alpha, dual_lambda, dual_mu, f}, DiscreteFunctionType::P0);
+
+ switch (i_mesh->dimension()) {
+ case 1: {
+ using MeshType = Mesh<Connectivity<1>>;
+ using DiscreteScalarFunctionType = DiscreteFunctionP0<1, double>;
+ using DiscreteVectorFunctionType = DiscreteFunctionP0<1, TinyVector<1>>;
+
+ std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+
+ if (DiamondDualMeshManager::instance().getDiamondDualMesh(mesh) != dual_mu->mesh()) {
+ throw NormalError("mu_dual is not defined on the dual mesh");
+ }
+
+ m_scheme =
+ std::make_unique<VectorDiamondScheme<1>>(mesh, std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
+ std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_lambda),
+ std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mu),
+ std::dynamic_pointer_cast<const DiscreteVectorFunctionType>(f),
+ bc_descriptor_list);
+ break;
+ }
+ case 2: {
+ using MeshType = Mesh<Connectivity<2>>;
+ using DiscreteScalarFunctionType = DiscreteFunctionP0<2, double>;
+ using DiscreteVectorFunctionType = DiscreteFunctionP0<2, TinyVector<2>>;
+
+ std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+
+ if (DiamondDualMeshManager::instance().getDiamondDualMesh(mesh) != dual_mu->mesh()) {
+ throw NormalError("mu_dual is not defined on the dual mesh");
+ }
+
+ m_scheme =
+ std::make_unique<VectorDiamondScheme<2>>(mesh, std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
+ std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_lambda),
+ std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mu),
+ std::dynamic_pointer_cast<const DiscreteVectorFunctionType>(f),
+ bc_descriptor_list);
+ break;
+ }
+ case 3: {
+ using MeshType = Mesh<Connectivity<3>>;
+ using DiscreteScalarFunctionType = DiscreteFunctionP0<3, double>;
+ using DiscreteVectorFunctionType = DiscreteFunctionP0<3, TinyVector<3>>;
+
+ std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+
+ if (DiamondDualMeshManager::instance().getDiamondDualMesh(mesh) != dual_mu->mesh()) {
+ throw NormalError("mu_dual is not defined on the dual mesh");
+ }
+
+ m_scheme =
+ std::make_unique<VectorDiamondScheme<3>>(mesh, std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
+ std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_lambda),
+ std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mu),
+ std::dynamic_pointer_cast<const DiscreteVectorFunctionType>(f),
+ bc_descriptor_list);
+ break;
+ }
+ default: {
+ throw UnexpectedError("invalid mesh dimension");
+ }
+ }
+}
+
+VectorDiamondSchemeHandler::~VectorDiamondSchemeHandler() = default;
diff --git a/src/scheme/VectorDiamondScheme.hpp b/src/scheme/VectorDiamondScheme.hpp
index bb67cd15d..cc2bca603 100644
--- a/src/scheme/VectorDiamondScheme.hpp
+++ b/src/scheme/VectorDiamondScheme.hpp
@@ -21,6 +21,32 @@
#include <scheme/NeumannBoundaryConditionDescriptor.hpp>
#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+class VectorDiamondSchemeHandler
+{
+ private:
+ class IVectorDiamondScheme;
+
+ template <size_t Dimension>
+ class VectorDiamondScheme;
+
+ template <size_t Dimension>
+ class InterpolationWeightsManager;
+
+ public:
+ std::unique_ptr<IVectorDiamondScheme> m_scheme;
+
+ std::shared_ptr<const IDiscreteFunction> solution() const;
+
+ VectorDiamondSchemeHandler(
+ const std::shared_ptr<const IDiscreteFunction>& alpha,
+ const std::shared_ptr<const IDiscreteFunction>& lambda,
+ const std::shared_ptr<const IDiscreteFunction>& mu,
+ const std::shared_ptr<const IDiscreteFunction>& f,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list);
+
+ ~VectorDiamondSchemeHandler();
+};
+
template <size_t Dimension>
class VectorDiamondScheme
{
--
GitLab