diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index ca49d0b16565877fc526dc9c5f7a8740975582d4..d1fa4bf7ce816e054d36a216cf66fc1084a4c6c7 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -42,6 +42,7 @@
#include <scheme/IDiscreteFunctionDescriptor.hpp>
#include <scheme/NeumannBoundaryConditionDescriptor.hpp>
#include <scheme/ScalarDiamondScheme.hpp>
+#include <scheme/ScalarHybridScheme.hpp>
#include <scheme/ScalarNodalScheme.hpp>
#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
#include <scheme/VectorDiamondScheme.hpp>
@@ -637,6 +638,42 @@ SchemeModule::SchemeModule()
.solution();
}));
+ this->_addBuiltinFunction("hybridheat",
+ std::make_shared<BuiltinFunctionEmbedder<std::tuple<
+ std::shared_ptr<const IDiscreteFunction>,
+ 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::vector<std::shared_ptr<
+ const IBoundaryConditionDescriptor>>&,
+ const std::shared_ptr<const IDiscreteFunction>&,
+ const double&, const double&, const double&)>>(
+
+ [](const std::shared_ptr<const IDiscreteFunction>& kappa,
+ const std::shared_ptr<const IDiscreteFunction>& f,
+ const std::shared_ptr<const IDiscreteFunction>& f_prev,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_prev_descriptor_list,
+ const std::shared_ptr<const IDiscreteFunction>& P, const double& dt,
+ const double& cn_coeff,
+ const double& lambda) -> const std::tuple<std::shared_ptr<const IDiscreteFunction>,
+ std::shared_ptr<const IDiscreteFunction>> {
+ return ScalarHybridSchemeHandler{kappa,
+ f,
+ f_prev,
+ bc_descriptor_list,
+ bc_prev_descriptor_list,
+ P,
+ dt,
+ cn_coeff,
+ lambda}
+ .solution();
+ }));
+
this->_addBuiltinFunction("unsteadyelasticity",
std::make_shared<BuiltinFunctionEmbedder<
void(std::shared_ptr<const IMesh>,
diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index b7f170fb75eca85662f60696a7c58b5c798246dd..810c615f3e6cf39047bc868cb614c9bd6fe6e30a 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -10,4 +10,5 @@ add_library(
DiscreteFunctionVectorInterpoler.cpp
ScalarDiamondScheme.cpp
VectorDiamondScheme.cpp
- ScalarNodalScheme.cpp)
+ ScalarNodalScheme.cpp
+ ScalarHybridScheme.cpp)
diff --git a/src/scheme/ScalarHybridScheme.cpp b/src/scheme/ScalarHybridScheme.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..10d41c948a8e3ccacb6b3acca3e28b7e61f96b78
--- /dev/null
+++ b/src/scheme/ScalarHybridScheme.cpp
@@ -0,0 +1,1525 @@
+#include <scheme/ScalarHybridScheme.hpp>
+
+#include <scheme/DiscreteFunctionP0.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
+
+class ScalarHybridSchemeHandler::IScalarHybridScheme
+{
+ public:
+ virtual std::tuple<std::shared_ptr<const IDiscreteFunction>, std::shared_ptr<const IDiscreteFunction>> getSolution()
+ const = 0;
+
+ // virtual std::shared_ptr<const IDiscreteFunction> getSolution() const = 0;
+
+ IScalarHybridScheme() = default;
+ virtual ~IScalarHybridScheme() = default;
+};
+
+template <size_t Dimension>
+class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeHandler::IScalarHybridScheme
+{
+ private:
+ using ConnectivityType = Connectivity<Dimension>;
+ using MeshType = Mesh<ConnectivityType>;
+ using MeshDataType = MeshData<Dimension>;
+
+ std::shared_ptr<const DiscreteFunctionP0<Dimension, double>> m_cell_temperature;
+ std::shared_ptr<const DiscreteFunctionP0<Dimension, double>> m_node_temperature;
+ // std::shared_ptr<const DiscreteFunctionP0<Dimension, double>> m_solution;
+
+ class DirichletBoundaryCondition
+ {
+ private:
+ const Array<const double> m_node_value_list;
+ const Array<const double> m_face_value_list;
+ const Array<const FaceId> m_face_list;
+ const Array<const NodeId> m_node_list;
+
+ public:
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_node_list;
+ }
+
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_face_list;
+ }
+
+ const Array<const double>&
+ nodeValueList() const
+ {
+ return m_node_value_list;
+ }
+
+ const Array<const double>&
+ faceValueList() const
+ {
+ return m_face_value_list;
+ }
+
+ DirichletBoundaryCondition(const Array<const FaceId>& face_list,
+ const Array<const NodeId>& node_list,
+ const Array<const double>& node_value_list,
+ const Array<const double>& face_value_list)
+ : m_node_value_list{node_value_list},
+ m_face_value_list{face_value_list},
+ m_face_list{face_list},
+ m_node_list{node_list}
+ {
+ Assert(m_node_value_list.size() == m_node_list.size());
+ Assert(m_face_value_list.size() == m_face_list.size());
+ }
+
+ ~DirichletBoundaryCondition() = default;
+ };
+
+ class NeumannBoundaryCondition
+ {
+ private:
+ const Array<const double> m_value_list;
+ const Array<const FaceId> m_face_list;
+ const Array<const NodeId> m_node_list;
+
+ public:
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_node_list;
+ }
+
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_face_list;
+ }
+
+ const Array<const double>&
+ valueList() const
+ {
+ return m_value_list;
+ }
+
+ NeumannBoundaryCondition(const Array<const FaceId>& face_list,
+ const Array<const NodeId>& node_list,
+ const Array<const double>& value_list)
+ : m_value_list{value_list}, m_face_list{face_list}, m_node_list{node_list}
+ {
+ Assert(m_value_list.size() == m_face_list.size());
+ }
+
+ ~NeumannBoundaryCondition() = default;
+ };
+
+ class FourierBoundaryCondition
+ {
+ private:
+ const Array<const double> m_coef_list;
+ const Array<const double> m_value_list;
+ const Array<const FaceId> m_face_list;
+
+ public:
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_face_list;
+ }
+
+ const Array<const double>&
+ valueList() const
+ {
+ return m_value_list;
+ }
+
+ const Array<const double>&
+ coefList() const
+ {
+ return m_coef_list;
+ }
+
+ public:
+ FourierBoundaryCondition(const Array<const FaceId>& face_list,
+ const Array<const double>& coef_list,
+ const Array<const double>& value_list)
+ : m_coef_list{coef_list}, m_value_list{value_list}, m_face_list{face_list}
+ {
+ Assert(m_coef_list.size() == m_face_list.size());
+ Assert(m_value_list.size() == m_face_list.size());
+ }
+
+ ~FourierBoundaryCondition() = default;
+ };
+
+ public:
+ std::tuple<std::shared_ptr<const IDiscreteFunction>, std::shared_ptr<const IDiscreteFunction>>
+ getSolution() const
+ {
+ return {m_cell_temperature, m_node_temperature};
+ }
+ // std::shared_ptr<const IDiscreteFunction>
+ // getSolution() const final
+ // {
+ // return m_solution;
+ // }
+
+ ScalarHybridScheme(const std::shared_ptr<const MeshType>& mesh,
+ const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyMatrix<Dimension>>>& cell_k,
+ const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& f,
+ const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& f_prev,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_prev_descriptor_list,
+ const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& P,
+ const double& dt,
+ const double& cn_coeff,
+ const double& lambda)
+ {
+ using BoundaryCondition =
+ std::variant<DirichletBoundaryCondition, FourierBoundaryCondition, NeumannBoundaryCondition>;
+
+ using BoundaryConditionList = std::vector<BoundaryCondition>;
+
+ BoundaryConditionList boundary_condition_list;
+
+ for (const auto& bc_descriptor : bc_descriptor_list) {
+ bool is_valid_boundary_condition = true;
+
+ switch (bc_descriptor->type()) {
+ case IBoundaryConditionDescriptor::Type::dirichlet: {
+ const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
+ dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
+ if constexpr (Dimension > 1) {
+ MeshFaceBoundary<Dimension> mesh_face_boundary =
+ getMeshFaceBoundary(*mesh, dirichlet_bc_descriptor.boundaryDescriptor());
+ MeshNodeBoundary<Dimension> mesh_node_boundary =
+ getMeshNodeBoundary(*mesh, dirichlet_bc_descriptor.boundaryDescriptor());
+
+ const FunctionSymbolId g_id = dirichlet_bc_descriptor.rhsSymbolId();
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+
+ Array<const double> node_value_list =
+ InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::node>(g_id, mesh->xr(),
+ mesh_node_boundary
+ .nodeList());
+ Array<const double> face_value_list =
+ InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id,
+ mesh_data.xl(),
+ mesh_face_boundary
+ .faceList());
+ boundary_condition_list.push_back(DirichletBoundaryCondition{mesh_face_boundary.faceList(),
+ mesh_node_boundary.nodeList(), node_value_list,
+ face_value_list});
+ } else {
+ throw NotImplementedError("Dirichlet BC in 1d");
+ }
+ break;
+ }
+ case IBoundaryConditionDescriptor::Type::fourier: {
+ throw NotImplementedError("NIY");
+ break;
+ }
+ case IBoundaryConditionDescriptor::Type::neumann: {
+ const NeumannBoundaryConditionDescriptor& neumann_bc_descriptor =
+ dynamic_cast<const NeumannBoundaryConditionDescriptor&>(*bc_descriptor);
+
+ if constexpr (Dimension > 1) {
+ MeshFaceBoundary<Dimension> mesh_face_boundary =
+ getMeshFaceBoundary(*mesh, neumann_bc_descriptor.boundaryDescriptor());
+ MeshNodeBoundary<Dimension> mesh_node_boundary =
+ getMeshNodeBoundary(*mesh, neumann_bc_descriptor.boundaryDescriptor());
+
+ const FunctionSymbolId g_id = neumann_bc_descriptor.rhsSymbolId();
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+
+ Array<const double> value_list =
+ InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id,
+ mesh_data.xl(),
+ mesh_face_boundary
+ .faceList());
+
+ boundary_condition_list.push_back(
+ NeumannBoundaryCondition{mesh_face_boundary.faceList(), mesh_node_boundary.nodeList(), value_list});
+
+ } else {
+ throw NotImplementedError("Neumann BC in 1d");
+ }
+ 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 heat equation";
+ throw NormalError(error_msg.str());
+ }
+ }
+
+ BoundaryConditionList boundary_prev_condition_list;
+
+ for (const auto& bc_prev_descriptor : bc_prev_descriptor_list) {
+ bool is_valid_boundary_condition = true;
+
+ switch (bc_prev_descriptor->type()) {
+ case IBoundaryConditionDescriptor::Type::dirichlet: {
+ const DirichletBoundaryConditionDescriptor& dirichlet_bc_prev_descriptor =
+ dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_prev_descriptor);
+ if constexpr (Dimension > 1) {
+ MeshFaceBoundary<Dimension> mesh_face_boundary =
+ getMeshFaceBoundary(*mesh, dirichlet_bc_prev_descriptor.boundaryDescriptor());
+ MeshNodeBoundary<Dimension> mesh_node_boundary =
+ getMeshNodeBoundary(*mesh, dirichlet_bc_prev_descriptor.boundaryDescriptor());
+
+ const FunctionSymbolId g_id = dirichlet_bc_prev_descriptor.rhsSymbolId();
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+
+ Array<const double> node_value_list =
+ InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::node>(g_id, mesh->xr(),
+ mesh_node_boundary
+ .nodeList());
+ Array<const double> face_value_list =
+ InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id,
+ mesh_data.xl(),
+ mesh_face_boundary
+ .faceList());
+ boundary_prev_condition_list.push_back(DirichletBoundaryCondition{mesh_face_boundary.faceList(),
+ mesh_node_boundary.nodeList(),
+ node_value_list, face_value_list});
+
+ } else {
+ throw NotImplementedError("Dirichlet BC in 1d");
+ }
+ break;
+ }
+ case IBoundaryConditionDescriptor::Type::fourier: {
+ throw NotImplementedError("NIY");
+ break;
+ }
+ case IBoundaryConditionDescriptor::Type::neumann: {
+ const NeumannBoundaryConditionDescriptor& neumann_bc_prev_descriptor =
+ dynamic_cast<const NeumannBoundaryConditionDescriptor&>(*bc_prev_descriptor);
+
+ if constexpr (Dimension > 1) {
+ MeshFaceBoundary<Dimension> mesh_face_boundary =
+ getMeshFaceBoundary(*mesh, neumann_bc_prev_descriptor.boundaryDescriptor());
+ MeshNodeBoundary<Dimension> mesh_node_boundary =
+ getMeshNodeBoundary(*mesh, neumann_bc_prev_descriptor.boundaryDescriptor());
+
+ const FunctionSymbolId g_id = neumann_bc_prev_descriptor.rhsSymbolId();
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+
+ Array<const double> value_list =
+ InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id,
+ mesh_data.xl(),
+ mesh_face_boundary
+ .faceList());
+
+ boundary_prev_condition_list.push_back(
+ NeumannBoundaryCondition{mesh_face_boundary.faceList(), mesh_node_boundary.nodeList(), value_list});
+
+ } else {
+ throw NotImplementedError("Neumann BC in 1d");
+ }
+ break;
+ }
+ default: {
+ is_valid_boundary_condition = false;
+ }
+ }
+ if (not is_valid_boundary_condition) {
+ std::ostringstream error_msg;
+ error_msg << *bc_prev_descriptor << " is an invalid boundary condition for heat equation";
+ throw NormalError(error_msg.str());
+ }
+ }
+
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+
+ std::shared_ptr dual_mesh = DualMeshManager::instance().getMedianDualMesh(*mesh);
+
+ std::shared_ptr mapper =
+ DualConnectivityManager::instance().getPrimalToMedianDualConnectivityDataMapper(mesh->connectivity());
+
+ const NodeValue<const TinyVector<Dimension>>& xr = mesh->xr();
+
+ const FaceValue<const TinyVector<Dimension>>& xl = mesh_data.xl();
+ const CellValue<const TinyVector<Dimension>>& xj = mesh_data.xj();
+
+ const NodeValuePerCell<const TinyVector<Dimension>>& Cjr = mesh_data.Cjr();
+ const CellValue<const TinyMatrix<Dimension>> cell_kappaj = cell_k->cellValues();
+
+ const auto is_boundary_node = mesh->connectivity().isBoundaryNode();
+ const auto is_boundary_face = mesh->connectivity().isBoundaryFace();
+
+ const auto& face_to_cell_matrix = mesh->connectivity().faceToCellMatrix();
+ const auto& node_to_face_matrix = mesh->connectivity().nodeToFaceMatrix();
+ const auto& face_to_node_matrix = mesh->connectivity().faceToNodeMatrix();
+ const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
+ const auto& cell_to_face_matrix = mesh->connectivity().cellToFaceMatrix();
+ const auto& node_local_numbers_in_their_cells = mesh->connectivity().nodeLocalNumbersInTheirCells();
+ const auto& cell_local_numbers_in_their_faces = mesh->connectivity().cellLocalNumbersInTheirFaces();
+ const auto& face_local_numbers_in_their_cells = mesh->connectivity().faceLocalNumbersInTheirCells();
+ const CellValue<const double> Vj = mesh_data.Vj();
+ const auto& node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
+ const auto& nl = mesh_data.nl();
+
+ const NodeValue<const bool> node_is_neumann = [&] {
+ NodeValue<bool> compute_node_is_neumann{mesh->connectivity()};
+ compute_node_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, NeumannBoundaryCondition>) {
+ 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];
+ const auto& face_nodes = face_to_node_matrix[face_id];
+ for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
+ const NodeId node_id = face_nodes[i_node];
+ compute_node_is_neumann[node_id] = true;
+ }
+ }
+ }
+ },
+ boundary_condition);
+ }
+ return compute_node_is_neumann;
+ }();
+
+ const FaceValue<const bool> face_is_neumann = [&] {
+ FaceValue<bool> compute_face_is_neumann{mesh->connectivity()};
+ compute_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, NeumannBoundaryCondition>) {
+ 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];
+ compute_face_is_neumann[face_id] = true;
+ }
+ }
+ },
+ boundary_condition);
+ }
+ return compute_face_is_neumann;
+ }();
+
+ const NodeValue<const bool> node_is_dirichlet = [&] {
+ NodeValue<bool> compute_node_is_dirichlet{mesh->connectivity()};
+ compute_node_is_dirichlet.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& node_list = bc.nodeList();
+
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ const NodeId node_id = node_list[i_node];
+
+ compute_node_is_dirichlet[node_id] = true;
+ }
+ }
+ },
+ boundary_condition);
+ }
+ return compute_node_is_dirichlet;
+ }();
+
+ const FaceValue<const bool> face_is_dirichlet = [&] {
+ FaceValue<bool> compute_face_is_dirichlet{mesh->connectivity()};
+ compute_face_is_dirichlet.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();
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+ compute_face_is_dirichlet[face_id] = true;
+ }
+ }
+ },
+ boundary_condition);
+ }
+ return compute_face_is_dirichlet;
+ }();
+
+ const NodeValue<const bool> node_is_corner = [&] {
+ NodeValue<bool> compute_node_is_corner{mesh->connectivity()};
+ compute_node_is_corner.fill(false);
+ parallel_for(
+ mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ if (is_boundary_node[node_id]) {
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+ const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
+ for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+ const unsigned int i_node = node_local_number_in_its_cell[i_cell];
+ const CellId cell_id = node_to_cell[i_cell];
+ const auto& cell_nodes = cell_to_node_matrix[cell_id];
+ NodeId i_node_prev;
+ NodeId i_node_next;
+ if (i_node == 0) {
+ i_node_prev = cell_nodes.size() - 1;
+ } else {
+ i_node_prev = i_node - 1;
+ }
+ if (i_node == cell_nodes.size() - 1) {
+ i_node_next = 0;
+ } else {
+ i_node_next = i_node + 1;
+ }
+ const NodeId prev_node_id = cell_to_node_matrix[cell_id][i_node_prev];
+ const NodeId next_node_id = cell_to_node_matrix[cell_id][i_node_next];
+ if (is_boundary_node[prev_node_id] and is_boundary_node[next_node_id]) {
+ compute_node_is_corner[node_id] = true;
+ }
+ }
+ }
+ });
+ return compute_node_is_corner;
+ }();
+
+ const NodeValue<const bool> node_is_angle = [&] {
+ NodeValue<bool> compute_node_is_angle{mesh->connectivity()};
+ compute_node_is_angle.fill(false);
+ // for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
+ parallel_for(
+ mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ if (is_boundary_node[node_id]) {
+ const auto& node_to_face = node_to_face_matrix[node_id];
+
+ TinyVector<Dimension> n1 = zero;
+ TinyVector<Dimension> n2 = zero;
+
+ for (size_t i_face = 0; i_face < node_to_face.size(); ++i_face) {
+ FaceId face_id = node_to_face[i_face];
+ if (is_boundary_face[face_id]) {
+ if (l2Norm(n1) == 0) {
+ n1 = nl[face_id];
+ } else {
+ n2 = nl[face_id];
+ }
+ }
+ }
+ if (l2Norm(n1 - n2) > (1.E-15) and l2Norm(n1 + n2) > (1.E-15) and not node_is_corner[node_id]) {
+ compute_node_is_angle[node_id] = true;
+ }
+ // std::cout << node_id << " " << n1 << " " << n2 << " " << compute_node_is_angle[node_id] << "\n";
+ }
+ });
+ return compute_node_is_angle;
+ }();
+
+ const NodeValue<const TinyVector<Dimension>> exterior_normal = [&] {
+ NodeValue<TinyVector<Dimension>> compute_exterior_normal{mesh->connectivity()};
+ compute_exterior_normal.fill(zero);
+ for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
+ if (is_boundary_node[node_id]) {
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+ const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
+ for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+ const CellId cell_id = node_to_cell[i_cell];
+ const unsigned int i_node = node_local_number_in_its_cell[i_cell];
+ compute_exterior_normal[node_id] += Cjr(cell_id, i_node);
+ }
+ const double norm_exterior_normal = l2Norm(compute_exterior_normal[node_id]);
+ compute_exterior_normal[node_id] *= 1. / norm_exterior_normal;
+ // std::cout << node_id << " " << compute_exterior_normal[node_id] << "\n";
+ }
+ }
+ return compute_exterior_normal;
+ }();
+
+ // for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
+ // std::cout << node_id << " " << exterior_normal[node_id] << " " << node_is_angle[node_id] << "\n";
+ //};
+
+ 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 NodeValue<const TinyVector<Dimension>> normal_angle_base = [&] {
+ NodeValue<TinyVector<Dimension>> compute_normal_base{mesh->connectivity()};
+ compute_normal_base.fill(zero);
+ for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
+ if (node_is_angle[node_id]) {
+ const auto& node_to_face = node_to_face_matrix[node_id];
+ TinyMatrix<Dimension> A;
+ A(0, 0) = 1000;
+ A(0, 1) = 1000;
+ A(1, 0) = 1000;
+ A(1, 1) = 1000;
+ for (size_t i_face = 0; i_face < node_to_face.size(); ++i_face) {
+ const FaceId face_id = node_to_face[i_face];
+ if (is_boundary_face[face_id]) {
+ if (A(0, 0) == 1000) {
+ if (dot(nl[face_id], exterior_normal[node_id]) >= 0) {
+ A(0, 0) = nl[face_id][0] * mes_l[face_id];
+ A(1, 0) = nl[face_id][1] * mes_l[face_id];
+ } else {
+ A(0, 0) = -nl[face_id][0] * mes_l[face_id];
+ A(1, 0) = -nl[face_id][1] * mes_l[face_id];
+ }
+ } else {
+ if (dot(nl[face_id], exterior_normal[node_id]) >= 0) {
+ A(0, 1) = nl[face_id][0] * mes_l[face_id];
+ A(1, 1) = nl[face_id][1] * mes_l[face_id];
+ } else {
+ A(0, 1) = -nl[face_id][0] * mes_l[face_id];
+ A(1, 1) = -nl[face_id][1] * mes_l[face_id];
+ }
+ }
+ }
+ }
+ compute_normal_base[node_id] = inverse(A) * exterior_normal[node_id];
+ }
+ }
+ return compute_normal_base;
+ }();
+
+ const NodeValue<const TinyMatrix<Dimension>> node_kappar = [&] {
+ NodeValue<TinyMatrix<Dimension>> kappa{mesh->connectivity()};
+ kappa.fill(zero);
+ parallel_for(
+ mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ kappa[node_id] = zero;
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+ double weight = 0;
+ for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+ const CellId cell_id = node_to_cell[i_cell];
+ double local_weight = 1. / l2Norm(xr[node_id] - xj[cell_id]);
+ kappa[node_id] += local_weight * cell_kappaj[cell_id];
+ weight += local_weight;
+ }
+ kappa[node_id] = 1. / weight * kappa[node_id];
+ });
+ return kappa;
+ }();
+
+ const FaceValue<const TinyMatrix<Dimension>> face_kappal = [&] {
+ FaceValue<TinyMatrix<Dimension>> kappa{mesh->connectivity()};
+ kappa.fill(zero);
+ parallel_for(
+ mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
+ const auto& face_to_node = face_to_node_matrix[face_id];
+ kappa[face_id] = 0.5 * (node_kappar[face_to_node[0]] + node_kappar[face_to_node[1]]);
+ });
+ return kappa;
+ }();
+
+ const FaceValuePerCell<const TinyVector<Dimension>> xj1_xj2 = [&] {
+ FaceValuePerCell<TinyVector<Dimension>> compute_xj1_xj2{mesh->connectivity()};
+ compute_xj1_xj2.fill(zero);
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ const auto& cell_to_face = cell_to_face_matrix[cell_id];
+
+ for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
+ const FaceId face_id = cell_to_face[i_face];
+ const auto& face_to_cell = face_to_cell_matrix[face_id];
+ TinyVector<Dimension> xj1;
+ TinyVector<Dimension> xj2;
+
+ if (not is_boundary_face[face_id]) {
+ if (face_to_cell[0] == cell_id) {
+ xj1 = xj[face_to_cell[1]];
+ xj2 = xj[cell_id];
+ } else {
+ xj1 = xj[face_to_cell[0]];
+ xj2 = xj[cell_id];
+ }
+ } else {
+ if (face_to_cell[0] == cell_id) {
+ xj1 = xl[face_id];
+ xj2 = xj[cell_id];
+ } else {
+ xj1 = xl[face_id];
+ xj2 = xj[cell_id];
+ }
+ }
+ compute_xj1_xj2(cell_id, i_face) = xj1 - xj2;
+ }
+ }
+ return compute_xj1_xj2;
+ }();
+
+ const FaceValuePerCell<const TinyVector<Dimension>> xj1_xj2_orth = [&] {
+ FaceValuePerCell<TinyVector<Dimension>> compute_xj1_xj2_orth{mesh->connectivity()};
+ compute_xj1_xj2_orth.fill(zero);
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ const auto& cell_to_face = cell_to_face_matrix[cell_id];
+
+ for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
+ const FaceId face_id = cell_to_face[i_face];
+
+ compute_xj1_xj2_orth(cell_id, i_face)[0] = -xj1_xj2(cell_id, i_face)[1];
+ compute_xj1_xj2_orth(cell_id, i_face)[1] = xj1_xj2(cell_id, i_face)[0];
+ }
+ }
+ return compute_xj1_xj2_orth;
+ }();
+
+ const FaceValuePerCell<const TinyVector<Dimension>> normal_face_base = [&] {
+ FaceValuePerCell<TinyVector<Dimension>> compute_normal_base{mesh->connectivity()};
+ compute_normal_base.fill(zero);
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ const auto& cell_to_face = cell_to_face_matrix[cell_id];
+
+ for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
+ const FaceId face_id = cell_to_face[i_face];
+ TinyMatrix<Dimension> A;
+ TinyVector<Dimension> face_normal;
+
+ A(0, 0) = xj1_xj2_orth(cell_id, i_face)[0];
+ A(0, 1) = xj1_xj2(cell_id, i_face)[0];
+ A(1, 0) = xj1_xj2_orth(cell_id, i_face)[1];
+ A(1, 1) = xj1_xj2(cell_id, i_face)[1];
+
+ if (dot(nl[face_id], xj1_xj2(cell_id, i_face)) > 0) {
+ face_normal = nl[face_id];
+ } else {
+ face_normal = -nl[face_id];
+ }
+
+ compute_normal_base(cell_id, i_face) = inverse(A) * (face_kappal[face_id] * face_normal);
+ }
+ }
+
+ return compute_normal_base;
+ }();
+
+ const FaceValue<const double> face_boundary_values = [&] {
+ FaceValue<double> compute_face_boundary_values{mesh->connectivity()};
+ FaceValue<double> sum_mes_l{mesh->connectivity()};
+ compute_face_boundary_values.fill(0);
+ sum_mes_l.fill(0);
+ 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& face_value_list = bc.faceValueList();
+
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+
+ compute_face_boundary_values[face_id] = face_value_list[i_face];
+ }
+ } else if constexpr (std::is_same_v<T, NeumannBoundaryCondition>) {
+ 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];
+
+ compute_face_boundary_values[face_id] = value_list[i_face];
+ }
+ }
+ },
+ boundary_condition);
+ }
+ return compute_face_boundary_values;
+ }();
+
+ const FaceValue<const double> face_boundary_prev_values = [&] {
+ FaceValue<double> compute_face_boundary_values{mesh->connectivity()};
+ FaceValue<double> sum_mes_l{mesh->connectivity()};
+ compute_face_boundary_values.fill(0);
+ sum_mes_l.fill(0);
+ for (const auto& boundary_condition : boundary_prev_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& face_value_list = bc.faceValueList();
+
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+
+ compute_face_boundary_values[face_id] = face_value_list[i_face];
+ }
+ } else if constexpr (std::is_same_v<T, NeumannBoundaryCondition>) {
+ 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];
+
+ compute_face_boundary_values[face_id] = value_list[i_face];
+ }
+ }
+ },
+ boundary_condition);
+ }
+ return compute_face_boundary_values;
+ }();
+
+ const NodeValue<const double> node_boundary_values = [&] {
+ NodeValue<double> compute_node_boundary_values{mesh->connectivity()};
+ NodeValue<double> sum_mes_l{mesh->connectivity()};
+ compute_node_boundary_values.fill(0);
+ sum_mes_l.fill(0);
+ 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, NeumannBoundaryCondition>) {
+ 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];
+ const auto& face_nodes = face_to_node_matrix[face_id];
+ for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
+ const NodeId node_id = face_nodes[i_node];
+ if (not node_is_dirichlet[node_id]) {
+ if (node_is_angle[node_id]) {
+ const auto& node_to_face = node_to_face_matrix[node_id];
+ for (size_t i_face_node = 0; i_face_node < node_to_face.size(); ++i_face_node) {
+ const FaceId face_node_id = node_to_face[i_face_node];
+ if (face_id == face_node_id) {
+ compute_node_boundary_values[node_id] +=
+ value_list[i_face] * mes_l[face_id] * normal_angle_base[node_id][i_face_node];
+ }
+ }
+ } else {
+ compute_node_boundary_values[node_id] += value_list[i_face] * mes_l[face_id];
+ sum_mes_l[node_id] += mes_l[face_id];
+ }
+ }
+ }
+ }
+
+ } else if constexpr (std::is_same_v<T, DirichletBoundaryCondition>) {
+ const auto& node_list = bc.nodeList();
+ const auto& node_value_list = bc.nodeValueList();
+
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ const NodeId node_id = node_list[i_node];
+
+ compute_node_boundary_values[node_id] = node_value_list[i_node];
+ }
+ }
+ },
+ boundary_condition);
+ }
+ for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
+ if ((not node_is_dirichlet[node_id]) && (node_is_neumann[node_id]) and not node_is_angle[node_id]) {
+ compute_node_boundary_values[node_id] /= sum_mes_l[node_id];
+ }
+ // std::cout << node_id << " " << compute_node_boundary_values[node_id] << "\n";
+ }
+ return compute_node_boundary_values;
+ }();
+
+ const NodeValue<const double> node_boundary_prev_values = [&] {
+ NodeValue<double> compute_node_boundary_values{mesh->connectivity()};
+ NodeValue<double> sum_mes_l{mesh->connectivity()};
+ compute_node_boundary_values.fill(0);
+ sum_mes_l.fill(0);
+ for (const auto& boundary_condition : boundary_prev_condition_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr (std::is_same_v<T, NeumannBoundaryCondition>) {
+ 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];
+ const auto& face_nodes = face_to_node_matrix[face_id];
+ for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
+ const NodeId node_id = face_nodes[i_node];
+ if (not node_is_dirichlet[node_id]) {
+ if (node_is_angle[node_id]) {
+ const auto& node_to_face = node_to_face_matrix[node_id];
+ for (size_t i_face_node = 0; i_face_node < node_to_face.size(); ++i_face_node) {
+ const FaceId face_node_id = node_to_face[i_face_node];
+ if (face_id == face_node_id) {
+ compute_node_boundary_values[node_id] +=
+ value_list[i_face] * mes_l[face_id] * normal_angle_base[node_id][i_face_node];
+ }
+ }
+ } else {
+ compute_node_boundary_values[node_id] += value_list[i_face] * mes_l[face_id];
+ sum_mes_l[node_id] += mes_l[face_id];
+ }
+ }
+ }
+ }
+ } else if constexpr (std::is_same_v<T, DirichletBoundaryCondition>) {
+ const auto& node_list = bc.nodeList();
+ const auto& node_value_list = bc.nodeValueList();
+
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ const NodeId node_id = node_list[i_node];
+
+ compute_node_boundary_values[node_id] = node_value_list[i_node];
+ }
+ }
+ },
+ boundary_condition);
+ }
+ for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
+ if ((not node_is_dirichlet[node_id]) && (node_is_neumann[node_id]) and not node_is_angle[node_id]) {
+ compute_node_boundary_values[node_id] /= sum_mes_l[node_id];
+ }
+ }
+ return compute_node_boundary_values;
+ }();
+
+ const NodeValue<const TinyMatrix<Dimension>> node_betar = [&] {
+ NodeValue<TinyMatrix<Dimension>> beta{mesh->connectivity()};
+ beta.fill(zero);
+ parallel_for(
+ mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+ beta[node_id] = zero;
+ for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+ const CellId cell_id = node_to_cell[i_cell];
+ const unsigned int i_node = node_local_number_in_its_cell[i_cell];
+ beta[node_id] += tensorProduct(Cjr(cell_id, i_node), xr[node_id] - xj[cell_id]);
+ }
+ });
+ return beta;
+ }();
+
+ const NodeValue<const TinyMatrix<Dimension>> kappar_invBetar = [&] {
+ NodeValue<TinyMatrix<Dimension>> kappa_invBeta{mesh->connectivity()};
+ parallel_for(
+ mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ if (not node_is_corner[node_id]) {
+ kappa_invBeta[node_id] = node_kappar[node_id] * inverse(node_betar[node_id]);
+ }
+ });
+ return kappa_invBeta;
+ }();
+
+ const NodeValue<const TinyMatrix<Dimension>> corner_betar = [&] {
+ NodeValue<TinyMatrix<Dimension>> beta{mesh->connectivity()};
+ parallel_for(
+ mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ if (node_is_corner[node_id]) {
+ const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+
+ size_t i_cell = 0;
+ const CellId cell_id = node_to_cell[i_cell];
+ const unsigned int i_node = node_local_number_in_its_cell[i_cell];
+
+ const auto& cell_nodes = cell_to_node_matrix[cell_id];
+ const NodeId node_id_p1 = cell_to_node_matrix[cell_id][(i_node + 1) % cell_nodes.size()];
+ const NodeId node_id_p2 = cell_to_node_matrix[cell_id][(i_node + 2) % cell_nodes.size()];
+ const NodeId node_id_m1 = cell_to_node_matrix[cell_id][(i_node - 1) % cell_nodes.size()];
+
+ const TinyVector<Dimension> xj1 = 1. / 3. * (xr[node_id] + xr[node_id_m1] + xr[node_id_p2]);
+ const TinyVector<Dimension> xj2 = 1. / 3. * (xr[node_id] + xr[node_id_p2] + xr[node_id_p1]);
+
+ const TinyVector<Dimension> xrm1 = xr[node_id_m1];
+ const TinyVector<Dimension> xrp1 = xr[node_id_p1];
+ const TinyVector<Dimension> xrp2 = xr[node_id_p2];
+
+ TinyVector<Dimension> Cjr1;
+ TinyVector<Dimension> Cjr2;
+
+ if (Dimension == 2) {
+ Cjr1[0] = -0.5 * (xrm1[1] - xrp2[1]);
+ Cjr1[1] = 0.5 * (xrm1[0] - xrp2[0]);
+ Cjr2[0] = -0.5 * (xrp2[1] - xrp1[1]);
+ Cjr2[1] = 0.5 * (xrp2[0] - xrp1[0]);
+ } else {
+ throw NotImplementedError("The scheme is not implemented in this dimension.");
+ }
+
+ beta[node_id] = tensorProduct(Cjr1, (xr[node_id] - xj1)) + tensorProduct(Cjr2, (xr[node_id] - xj2));
+ }
+ });
+ return beta;
+ }();
+
+ const NodeValue<const TinyMatrix<Dimension>> corner_kappar_invBetar = [&] {
+ NodeValue<TinyMatrix<Dimension>> kappa_invBeta{mesh->connectivity()};
+ parallel_for(
+ mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ if (node_is_corner[node_id]) {
+ kappa_invBeta[node_id] = node_kappar[node_id] * inverse(corner_betar[node_id]);
+ }
+ });
+ return kappa_invBeta;
+ }();
+
+ const NodeValue<const TinyVector<Dimension>> sum_Cjr = [&] {
+ NodeValue<TinyVector<Dimension>> compute_sum_Cjr{mesh->connectivity()};
+ compute_sum_Cjr.fill(zero);
+ parallel_for(
+ mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+ const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
+ compute_sum_Cjr[node_id] = zero;
+ for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+ const CellId cell_id = node_to_cell[i_cell];
+ const size_t i_node = node_local_number_in_its_cell[i_cell];
+ compute_sum_Cjr[node_id] += Cjr(cell_id, i_node);
+ }
+ });
+ return compute_sum_Cjr;
+ }();
+
+ const NodeValue<const TinyVector<Dimension>> v = [&] {
+ NodeValue<TinyVector<Dimension>> compute_v{mesh->connectivity()};
+ parallel_for(
+ mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ if (is_boundary_node[node_id]) {
+ compute_v[node_id] = 1. / l2Norm(node_kappar[node_id] * exterior_normal[node_id]) * node_kappar[node_id] *
+ exterior_normal[node_id];
+ }
+ });
+ return compute_v;
+ }();
+
+ const NodeValuePerCell<const double> theta = [&] {
+ NodeValuePerCell<double> compute_theta{mesh->connectivity()};
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+ const 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];
+ if (is_boundary_node[node_id] && not node_is_corner[node_id]) {
+ compute_theta(cell_id, i_node) = dot(inverse(node_betar[node_id]) * Cjr(cell_id, i_node), v[node_id]);
+ }
+ }
+ });
+ return compute_theta;
+ }();
+
+ const NodeValue<const double> sum_theta = [&] {
+ NodeValue<double> compute_sum_theta{mesh->connectivity()};
+ compute_sum_theta.fill(0);
+ parallel_for(
+ mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ if ((is_boundary_node[node_id]) && (not node_is_corner[node_id])) {
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+ const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
+ compute_sum_theta[node_id] = 0;
+ for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+ const CellId cell_id = node_to_cell[i_cell];
+ const size_t i_node = node_local_number_in_its_cell[i_cell];
+ compute_sum_theta[node_id] += theta(cell_id, i_node);
+ }
+ }
+ });
+ return compute_sum_theta;
+ }();
+
+ const Array<int> non_zeros{mesh->numberOfCells()};
+ non_zeros.fill(4); // Modif pour optimiser
+ CRSMatrixDescriptor<double> S1(mesh->numberOfCells(), mesh->numberOfCells(), non_zeros);
+ CRSMatrixDescriptor<double> S2(mesh->numberOfCells(), mesh->numberOfCells(), non_zeros);
+
+ for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+ for (size_t i_cell_j = 0; i_cell_j < node_to_cell.size(); ++i_cell_j) {
+ const CellId cell_id_j = node_to_cell[i_cell_j];
+ for (size_t i_cell_k = 0; i_cell_k < node_to_cell.size(); ++i_cell_k) {
+ const CellId cell_id_k = node_to_cell[i_cell_k];
+ S1(cell_id_j, cell_id_k) = 0;
+ S2(cell_id_j, cell_id_k) = 0;
+ }
+ }
+ }
+
+ for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+ const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
+
+ if (not is_boundary_node[node_id]) {
+ for (size_t i_cell_j = 0; i_cell_j < node_to_cell.size(); ++i_cell_j) {
+ const CellId cell_id_j = node_to_cell[i_cell_j];
+ const size_t i_node_j = node_local_number_in_its_cells[i_cell_j];
+
+ for (size_t i_cell_k = 0; i_cell_k < node_to_cell.size(); ++i_cell_k) {
+ const CellId cell_id_k = node_to_cell[i_cell_k];
+ const size_t i_node_k = node_local_number_in_its_cells[i_cell_k];
+
+ S1(cell_id_j, cell_id_k) +=
+ (1 - lambda) * dt * (1. - cn_coeff / 2.) *
+ dot(kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j));
+ S2(cell_id_j, cell_id_k) -=
+ (1 - lambda) * dt * (cn_coeff / 2.) *
+ dot(kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j));
+
+ const auto& cell_to_face = cell_to_face_matrix[cell_id_j];
+
+ for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
+ const FaceId face_id = cell_to_face[i_face];
+ const auto& face_to_node = face_to_node_matrix[face_id];
+
+ if (face_to_node[0] == node_id or face_to_node[1] == node_id) {
+ S1(cell_id_j, cell_id_k) +=
+ lambda * dt * (1. - cn_coeff / 2.) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
+ dot(inverse(node_betar[node_id]) * Cjr(cell_id_k, i_node_k), xj1_xj2_orth(cell_id_j, i_node_j));
+ S2(cell_id_j, cell_id_k) -=
+ lambda * dt * (cn_coeff / 2.) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
+ dot(inverse(node_betar[node_id]) * Cjr(cell_id_k, i_node_k), xj1_xj2_orth(cell_id_j, i_node_j));
+ }
+ }
+ }
+ }
+ } else if ((node_is_neumann[node_id]) && (not node_is_corner[node_id]) && (not node_is_dirichlet[node_id])) {
+ const TinyMatrix<Dimension> I = identity;
+ const TinyVector<Dimension> n = exterior_normal[node_id];
+ const TinyMatrix<Dimension> nxn = tensorProduct(n, n);
+ const TinyMatrix<Dimension> Q = I - nxn;
+
+ for (size_t i_cell_j = 0; i_cell_j < node_to_cell.size(); ++i_cell_j) {
+ const CellId cell_id_j = node_to_cell[i_cell_j];
+ const size_t i_node_j = node_local_number_in_its_cells[i_cell_j];
+
+ for (size_t i_cell_k = 0; i_cell_k < node_to_cell.size(); ++i_cell_k) {
+ const CellId cell_id_k = node_to_cell[i_cell_k];
+ const size_t i_node_k = node_local_number_in_its_cells[i_cell_k];
+
+ S1(cell_id_j, cell_id_k) +=
+ (1 - lambda) * dt * (1. - cn_coeff / 2.) *
+ dot(kappar_invBetar[node_id] *
+ (Cjr(cell_id_k, i_node_k) - theta(cell_id_k, i_node_k) / sum_theta[node_id] * sum_Cjr[node_id]),
+ Q * Cjr(cell_id_j, i_node_j));
+ S2(cell_id_j, cell_id_k) -=
+ (1 - lambda) * dt * (cn_coeff / 2.) *
+ dot(kappar_invBetar[node_id] *
+ (Cjr(cell_id_k, i_node_k) - theta(cell_id_k, i_node_k) / sum_theta[node_id] * sum_Cjr[node_id]),
+ Q * Cjr(cell_id_j, i_node_j));
+ }
+ }
+
+ } else if ((node_is_dirichlet[node_id]) && (not node_is_corner[node_id])) {
+ for (size_t i_cell_j = 0; i_cell_j < node_to_cell.size(); ++i_cell_j) {
+ const CellId cell_id_j = node_to_cell[i_cell_j];
+ const size_t i_node_j = node_local_number_in_its_cells[i_cell_j];
+
+ for (size_t i_cell_k = 0; i_cell_k < node_to_cell.size(); ++i_cell_k) {
+ const CellId cell_id_k = node_to_cell[i_cell_k];
+ const size_t i_node_k = node_local_number_in_its_cells[i_cell_k];
+
+ S1(cell_id_j, cell_id_k) +=
+ (1 - lambda) * dt * (1. - cn_coeff / 2.) *
+ dot(kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j));
+ S2(cell_id_j, cell_id_k) -=
+ (1 - lambda) * dt * (cn_coeff / 2.) *
+ dot(kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j));
+
+ const auto& cell_to_face = cell_to_face_matrix[cell_id_j];
+
+ for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
+ const FaceId face_id = cell_to_face[i_face];
+ const auto& face_to_node = face_to_node_matrix[face_id];
+
+ if (face_to_node[0] == node_id or face_to_node[1] == node_id) {
+ S1(cell_id_j, cell_id_k) +=
+ lambda * dt * (1. - cn_coeff / 2.) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
+ dot(inverse(node_betar[node_id]) * Cjr(cell_id_k, i_node_k), xj1_xj2_orth(cell_id_j, i_node_j));
+ S2(cell_id_j, cell_id_k) -=
+ lambda * dt * (cn_coeff / 2.) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
+ dot(inverse(node_betar[node_id]) * Cjr(cell_id_k, i_node_k), xj1_xj2_orth(cell_id_j, i_node_j));
+ }
+ }
+ }
+ }
+
+ } else if (node_is_dirichlet[node_id] && node_is_corner[node_id]) {
+ for (size_t i_cell_j = 0; i_cell_j < node_to_cell.size(); ++i_cell_j) {
+ const CellId cell_id_j = node_to_cell[i_cell_j];
+ const size_t i_node_j = node_local_number_in_its_cells[i_cell_j];
+
+ for (size_t i_cell_k = 0; i_cell_k < node_to_cell.size(); ++i_cell_k) {
+ const CellId cell_id_k = node_to_cell[i_cell_k];
+ const size_t i_node_k = node_local_number_in_its_cells[i_cell_k];
+
+ S1(cell_id_j, cell_id_k) +=
+ (1 - lambda) * dt * (1. - cn_coeff / 2.) *
+ dot(corner_kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j));
+ S2(cell_id_j, cell_id_k) -=
+ (1 - lambda) * dt * (cn_coeff / 2.) *
+ dot(corner_kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j));
+
+ const auto& cell_to_face = cell_to_face_matrix[cell_id_j];
+
+ for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
+ const FaceId face_id = cell_to_face[i_face];
+ const auto& face_to_node = face_to_node_matrix[face_id];
+
+ if (face_to_node[0] == node_id or face_to_node[1] == node_id) {
+ S1(cell_id_j, cell_id_k) +=
+ lambda * dt * (1. - cn_coeff / 2.) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
+ dot(inverse(corner_betar[node_id]) * Cjr(cell_id_k, i_node_k), xj1_xj2_orth(cell_id_j, i_node_j));
+
+ S2(cell_id_j, cell_id_k) -=
+ lambda * dt * (cn_coeff / 2.) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
+ dot(inverse(corner_betar[node_id]) * Cjr(cell_id_k, i_node_k), xj1_xj2_orth(cell_id_j, i_node_j));
+ }
+ }
+ }
+ }
+ }
+ }
+
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ const auto& face_to_cell = face_to_cell_matrix[face_id];
+ const auto& face_local_number_in_its_cells = face_local_numbers_in_their_cells.itemArray(face_id);
+ if (not is_boundary_face[face_id]) {
+ for (size_t i_cell_j = 0; i_cell_j < face_to_cell.size(); ++i_cell_j) {
+ const CellId cell_id_j = face_to_cell[i_cell_j];
+ const size_t i_face = face_local_number_in_its_cells[i_cell_j];
+
+ for (size_t i_cell_k = 0; i_cell_k < face_to_cell.size(); ++i_cell_k) {
+ const CellId cell_id_k = face_to_cell[i_cell_k];
+
+ if (cell_id_j == cell_id_k) {
+ S1(cell_id_j, cell_id_k) -=
+ lambda * dt * (1. - cn_coeff / 2.) * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[1];
+ S2(cell_id_j, cell_id_k) +=
+ lambda * dt * (cn_coeff / 2.) * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[1];
+ } else {
+ S1(cell_id_j, cell_id_k) +=
+ lambda * dt * (1. - cn_coeff / 2.) * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[1];
+ S2(cell_id_j, cell_id_k) -=
+ lambda * dt * (cn_coeff / 2.) * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[1];
+ }
+ }
+ }
+ } else if (face_is_dirichlet[face_id]) {
+ for (size_t i_cell_j = 0; i_cell_j < face_to_cell.size(); ++i_cell_j) {
+ const CellId cell_id_j = face_to_cell[i_cell_j];
+ const size_t i_face = face_local_number_in_its_cells[i_cell_j];
+ S1(cell_id_j, cell_id_j) +=
+ lambda * dt * (1. - cn_coeff / 2.) * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[1];
+ S2(cell_id_j, cell_id_j) -=
+ lambda * dt * (cn_coeff / 2.) * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[1];
+ }
+ }
+ };
+
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ S1(cell_id, cell_id) += Vj[cell_id];
+ S2(cell_id, cell_id) += Vj[cell_id];
+ }
+
+ CellValue<const double> fj = f->cellValues();
+ CellValue<const double> fj_prev = f_prev->cellValues();
+ CellValue<const double> Pj = P->cellValues();
+
+ Vector<double> Ph{mesh->numberOfCells()};
+ Ph = zero;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ Ph[cell_id] = Pj[cell_id];
+ };
+
+ std::cout << "P = " << Ph << "\n";
+
+ CRSMatrix A1{S1.getCRSMatrix()};
+ CRSMatrix A2{S2.getCRSMatrix()};
+
+ // std::cout << A1 << "\n";
+
+ Vector<double> b{mesh->numberOfCells()};
+ b = zero;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ b[cell_id] = Vj[cell_id] * fj[cell_id];
+ };
+
+ Vector<double> b_prev{mesh->numberOfCells()};
+ b_prev = zero;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ b_prev[cell_id] = Vj[cell_id] * fj_prev[cell_id];
+ };
+
+ for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+ const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
+ const TinyMatrix<Dimension> I = identity;
+ const TinyVector<Dimension> n = exterior_normal[node_id];
+ const TinyMatrix<Dimension> nxn = tensorProduct(n, n);
+ const TinyMatrix<Dimension> Q = I - nxn;
+ if ((node_is_neumann[node_id]) && (not node_is_dirichlet[node_id])) {
+ if ((is_boundary_node[node_id]) and (not node_is_corner[node_id])) {
+ for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+ const CellId cell_id = node_to_cell[i_cell];
+ const size_t i_node = node_local_number_in_its_cells[i_cell];
+
+ b[cell_id] += (1 - lambda) * node_boundary_values[node_id] *
+ (1. / (sum_theta[node_id] * l2Norm(node_kappar[node_id] * n)) *
+ dot(kappar_invBetar[node_id] * sum_Cjr[node_id], Q * Cjr(cell_id, i_node)) +
+ dot(Cjr(cell_id, i_node), n));
+ b_prev[cell_id] += (1 - lambda) * node_boundary_prev_values[node_id] *
+ (1. / (sum_theta[node_id] * l2Norm(node_kappar[node_id] * n)) *
+ dot(kappar_invBetar[node_id] * sum_Cjr[node_id], Q * Cjr(cell_id, i_node)) +
+ dot(Cjr(cell_id, i_node), n));
+ }
+ } else if (node_is_corner[node_id]) {
+ const auto& node_to_face = node_to_face_matrix[node_id];
+ const CellId cell_id = node_to_cell[0];
+ double sum_mes_l = 0;
+ for (size_t i_face = 0; i_face < node_to_face.size(); ++i_face) {
+ FaceId face_id = node_to_face[i_face];
+ sum_mes_l += mes_l[face_id];
+ }
+ b[cell_id] += (1 - lambda) * 0.5 * node_boundary_values[node_id] * sum_mes_l;
+ b_prev[cell_id] += (1 - lambda) * 0.5 * node_boundary_prev_values[node_id] * sum_mes_l;
+ }
+
+ } else if (node_is_dirichlet[node_id]) {
+ if (not node_is_corner[node_id]) {
+ for (size_t i_cell_j = 0; i_cell_j < node_to_cell.size(); ++i_cell_j) {
+ const CellId cell_id_j = node_to_cell[i_cell_j];
+ const size_t i_node_j = node_local_number_in_its_cells[i_cell_j];
+
+ for (size_t i_cell_k = 0; i_cell_k < node_to_cell.size(); ++i_cell_k) {
+ const CellId cell_id_k = node_to_cell[i_cell_k];
+ const size_t i_node_k = node_local_number_in_its_cells[i_cell_k];
+
+ b[cell_id_j] +=
+ (1 - lambda) * dot(node_boundary_values[node_id] * kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k),
+ Cjr(cell_id_j, i_node_j));
+ b_prev[cell_id_j] += (1 - lambda) * dot(node_boundary_prev_values[node_id] * kappar_invBetar[node_id] *
+ Cjr(cell_id_k, i_node_k),
+ Cjr(cell_id_j, i_node_j));
+ }
+ const auto& cell_to_face = cell_to_face_matrix[cell_id_j];
+
+ for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
+ const FaceId face_id = cell_to_face[i_face];
+ const auto& face_to_node = face_to_node_matrix[face_id];
+
+ if (face_to_node[0] == node_id or face_to_node[1] == node_id) {
+ b[cell_id_j] += (1 - lambda) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
+ face_boundary_values[face_id] *
+ dot(inverse(node_betar[node_id]) * sum_Cjr[node_id], xj1_xj2_orth(cell_id_j, i_node_j));
+ b_prev[cell_id_j] +=
+ (1 - lambda) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
+ face_boundary_prev_values[face_id] *
+ dot(inverse(node_betar[node_id]) * sum_Cjr[node_id], xj1_xj2_orth(cell_id_j, i_node_j));
+ }
+ }
+ }
+
+ } else if (node_is_corner[node_id]) {
+ for (size_t i_cell_j = 0; i_cell_j < node_to_cell.size(); ++i_cell_j) {
+ const CellId cell_id_j = node_to_cell[i_cell_j];
+ const size_t i_node_j = node_local_number_in_its_cells[i_cell_j];
+
+ for (size_t i_cell_k = 0; i_cell_k < node_to_cell.size(); ++i_cell_k) {
+ const CellId cell_id_k = node_to_cell[i_cell_k];
+ const size_t i_node_k = node_local_number_in_its_cells[i_cell_k];
+
+ b[cell_id_j] += (1 - lambda) * dot(node_boundary_values[node_id] * corner_kappar_invBetar[node_id] *
+ Cjr(cell_id_k, i_node_k),
+ Cjr(cell_id_j, i_node_j));
+ b_prev[cell_id_j] += (1 - lambda) * dot(node_boundary_prev_values[node_id] *
+ corner_kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k),
+ Cjr(cell_id_j, i_node_j));
+ }
+ const auto& cell_to_face = cell_to_face_matrix[cell_id_j];
+
+ for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
+ const FaceId face_id = cell_to_face[i_face];
+ const auto& face_to_node = face_to_node_matrix[face_id];
+
+ if (face_to_node[0] == node_id or face_to_node[1] == node_id) {
+ b[cell_id_j] +=
+ (1 - lambda) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
+ face_boundary_values[face_id] *
+ dot(inverse(corner_betar[node_id]) * sum_Cjr[node_id], xj1_xj2_orth(cell_id_j, i_node_j));
+ b_prev[cell_id_j] +=
+ (1 - lambda) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
+ face_boundary_prev_values[face_id] *
+ dot(inverse(corner_betar[node_id]) * sum_Cjr[node_id], xj1_xj2_orth(cell_id_j, i_node_j));
+ }
+ }
+ }
+ }
+ }
+ };
+
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ const auto& face_to_cell = face_to_cell_matrix[face_id];
+ const auto& face_local_number_in_its_cells = face_local_numbers_in_their_cells.itemArray(face_id);
+ if (face_is_dirichlet[face_id]) {
+ for (size_t i_cell_j = 0; i_cell_j < face_to_cell.size(); ++i_cell_j) {
+ const CellId cell_id_j = face_to_cell[i_cell_j];
+ const size_t i_face = face_local_number_in_its_cells[i_cell_j];
+
+ b[cell_id_j] +=
+ lambda * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[1] * face_boundary_values[face_id];
+ b_prev[cell_id_j] +=
+ lambda * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[1] * face_boundary_prev_values[face_id];
+ }
+ } else if (face_is_neumann[face_id]) {
+ for (size_t i_cell_j = 0; i_cell_j < face_to_cell.size(); ++i_cell_j) {
+ const CellId cell_id_j = face_to_cell[i_cell_j];
+
+ b[cell_id_j] += lambda * mes_l[face_id] * face_boundary_values[face_id];
+ b_prev[cell_id_j] += lambda * mes_l[face_id] * face_boundary_prev_values[face_id];
+ }
+ }
+ }
+
+ Vector<double> T{mesh->numberOfCells()};
+ T = zero;
+
+ // Array<double> Ph2 = Ph;
+ Vector<double> A2P = A2 * Ph;
+
+ Vector<double> B{mesh->numberOfCells()};
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+ B[cell_id] = dt * ((1. - cn_coeff / 2.) * b[cell_id] + cn_coeff / 2. * b_prev[cell_id]) + A2P[cell_id];
+ });
+
+ std::cout << "B = " << B << "\n"
+ << "A2P = " << A2P << "\n";
+
+ NodeValue<TinyVector<Dimension>> ur{mesh->connectivity()};
+ ur.fill(zero);
+ for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+ const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
+ if (not is_boundary_node[node_id]) {
+ for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+ const CellId cell_id = node_to_cell[i_cell];
+ const size_t i_node = node_local_number_in_its_cell[i_cell];
+ ur[node_id] += Pj[cell_id] * Cjr(cell_id, i_node);
+ }
+ ur[node_id] = -inverse(node_betar[node_id]) * ur[node_id];
+ // std::cout << node_id << "; ur = " << ur[node_id] << "\n";
+ } else if (is_boundary_node[node_id] and node_is_dirichlet[node_id]) {
+ for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+ const CellId cell_id = node_to_cell[i_cell];
+ const size_t i_node = node_local_number_in_its_cell[i_cell];
+ ur[node_id] += (node_boundary_values[node_id] - Pj[cell_id]) * Cjr(cell_id, i_node);
+ }
+ if (not node_is_corner[node_id]) {
+ ur[node_id] = inverse(node_betar[node_id]) * ur[node_id];
+ } else {
+ ur[node_id] = inverse(corner_betar[node_id]) * ur[node_id];
+ }
+ // std::cout << "bord : " << node_id << "; ur = " << ur[node_id] << "\n";
+ }
+ }
+
+ LinearSolver solver;
+ solver.solveLocalSystem(A1, T, B);
+
+ const NodeValue<const double> primal_node_temperature = [&] {
+ NodeValue<double> compute_primal_node_temperature{mesh->connectivity()};
+ compute_primal_node_temperature.fill(0);
+ parallel_for(
+ mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ // for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); node_id++) {
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+ double sum_Vj = 0;
+ for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+ const CellId cell_id = node_to_cell[i_cell];
+ compute_primal_node_temperature[node_id] += Vj[cell_id] * T[cell_id];
+ sum_Vj += Vj[cell_id];
+ }
+ compute_primal_node_temperature[node_id] /= sum_Vj;
+ });
+ return compute_primal_node_temperature;
+ }();
+
+ const CellValue<const double> dual_node_temperature = [=] {
+ CellValue<double> my_dual_temperature{dual_mesh->connectivity()};
+ mapper->toDualCell(primal_node_temperature, my_dual_temperature);
+ return my_dual_temperature;
+ }();
+
+ // m_solution = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh);
+ m_cell_temperature = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh);
+ m_node_temperature = std::make_shared<DiscreteFunctionP0<Dimension, double>>(dual_mesh);
+ // auto& solution = *m_solution;
+ auto& cell_temperature = *m_cell_temperature;
+ auto& node_temperature = *m_node_temperature;
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { cell_temperature[cell_id] = T[cell_id]; });
+ parallel_for(
+ dual_mesh->numberOfCells(),
+ PUGS_LAMBDA(CellId cell_id) { node_temperature[cell_id] = dual_node_temperature[cell_id]; });
+ }
+};
+
+// std::shared_ptr<const IDiscreteFunction>
+// ScalarHybridSchemeHandler::solution() const
+// {
+// return m_scheme->getSolution();
+// }
+std::tuple<std::shared_ptr<const IDiscreteFunction>, std::shared_ptr<const IDiscreteFunction>>
+ScalarHybridSchemeHandler::solution() const
+{
+ return m_scheme->getSolution();
+}
+
+ScalarHybridSchemeHandler::ScalarHybridSchemeHandler(
+ const std::shared_ptr<const IDiscreteFunction>& cell_k,
+ const std::shared_ptr<const IDiscreteFunction>& f,
+ const std::shared_ptr<const IDiscreteFunction>& f_prev,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_prev_descriptor_list,
+ const std::shared_ptr<const IDiscreteFunction>& P,
+ const double& dt,
+ const double& cn_coeff,
+ const double& lambda)
+{
+ const std::shared_ptr i_mesh = getCommonMesh({cell_k, f});
+ if (not i_mesh) {
+ throw NormalError("primal discrete functions are not defined on the same mesh");
+ }
+
+ checkDiscretizationType({cell_k, f}, DiscreteFunctionType::P0);
+
+ switch (i_mesh->dimension()) {
+ case 1: {
+ throw NormalError("The scheme is not implemented in dimension 1");
+ break;
+ }
+ case 2: {
+ using MeshType = Mesh<Connectivity<2>>;
+ using DiscreteTensorFunctionType = DiscreteFunctionP0<2, TinyMatrix<2>>;
+ using DiscreteScalarFunctionType = DiscreteFunctionP0<2, double>;
+
+ std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+
+ m_scheme =
+ std::make_unique<ScalarHybridScheme<2>>(mesh, std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k),
+ std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
+ std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f_prev),
+ bc_descriptor_list, bc_prev_descriptor_list,
+ std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(P), dt,
+ cn_coeff, lambda);
+ break;
+ }
+ case 3: {
+ throw NormalError("The scheme is not implemented in dimension 3");
+ break;
+ }
+ default: {
+ throw UnexpectedError("invalid mesh dimension");
+ }
+ }
+}
+
+ScalarHybridSchemeHandler::~ScalarHybridSchemeHandler() = default;
diff --git a/src/scheme/ScalarHybridScheme.hpp b/src/scheme/ScalarHybridScheme.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..b110887d7bcce04d4456f8185bb6201ed983a8ae
--- /dev/null
+++ b/src/scheme/ScalarHybridScheme.hpp
@@ -0,0 +1,59 @@
+#ifndef SCALAR_HYBRID_SCHEME_HPP
+#define SCALAR_HYBRID_SCHEME_HPP
+
+#include <algebra/CRSMatrix.hpp>
+#include <algebra/CRSMatrixDescriptor.hpp>
+#include <algebra/LeastSquareSolver.hpp>
+#include <algebra/LinearSolver.hpp>
+#include <algebra/TinyVector.hpp>
+#include <algebra/Vector.hpp>
+#include <language/utils/InterpolateItemValue.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/DualConnectivityManager.hpp>
+#include <mesh/DualMeshManager.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <mesh/MeshFaceBoundary.hpp>
+#include <mesh/MeshNodeBoundary.hpp>
+#include <mesh/PrimalToDiamondDualConnectivityDataMapper.hpp>
+#include <mesh/PrimalToMedianDualConnectivityDataMapper.hpp>
+#include <mesh/SubItemValuePerItem.hpp>
+#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
+#include <scheme/FourierBoundaryConditionDescriptor.hpp>
+#include <scheme/IDiscreteFunction.hpp>
+#include <scheme/NeumannBoundaryConditionDescriptor.hpp>
+#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+
+class ScalarHybridSchemeHandler
+{
+ private:
+ class IScalarHybridScheme;
+
+ template <size_t Dimension>
+ class ScalarHybridScheme;
+
+ template <size_t Dimension>
+ class InterpolationWeightsManager;
+
+ public:
+ std::unique_ptr<IScalarHybridScheme> m_scheme;
+
+ // std::shared_ptr<const IDiscreteFunction> solution() const;
+ std::tuple<std::shared_ptr<const IDiscreteFunction>, std::shared_ptr<const IDiscreteFunction>> solution() const;
+
+ ScalarHybridSchemeHandler(
+ const std::shared_ptr<const IDiscreteFunction>& cell_k,
+ const std::shared_ptr<const IDiscreteFunction>& f,
+ const std::shared_ptr<const IDiscreteFunction>& f_prec,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_prev_descriptor_list,
+ const std::shared_ptr<const IDiscreteFunction>& P,
+ const double& dt,
+ const double& cn_coeff,
+ const double& lambda);
+
+ ~ScalarHybridSchemeHandler();
+};
+
+#endif // SCALAR_NODAL_SCHEME_HPP
diff --git a/src/scheme/ScalarNodalScheme.cpp b/src/scheme/ScalarNodalScheme.cpp
index c7fe9279c6fa0237f07cf803ae78dae518141595..0897f6d8323c2cccd368f5cedb23fcdf33b3c7a3 100644
--- a/src/scheme/ScalarNodalScheme.cpp
+++ b/src/scheme/ScalarNodalScheme.cpp
@@ -30,7 +30,8 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
class DirichletBoundaryCondition
{
private:
- const Array<const double> m_value_list;
+ const Array<const double> m_node_value_list;
+ const Array<const double> m_face_value_list;
const Array<const FaceId> m_face_list;
const Array<const NodeId> m_node_list;
@@ -48,17 +49,28 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
}
const Array<const double>&
- valueList() const
+ nodeValueList() const
{
- return m_value_list;
+ return m_node_value_list;
+ }
+
+ const Array<const double>&
+ faceValueList() const
+ {
+ return m_face_value_list;
}
DirichletBoundaryCondition(const Array<const FaceId>& face_list,
const Array<const NodeId>& node_list,
- const Array<const double>& value_list)
- : m_value_list{value_list}, m_face_list{face_list}, m_node_list{node_list}
+ const Array<const double>& node_value_list,
+ const Array<const double>& face_value_list)
+ : m_node_value_list{node_value_list},
+ m_face_value_list{face_value_list},
+ m_face_list{face_list},
+ m_node_list{node_list}
{
- Assert(m_value_list.size() == m_node_list.size());
+ Assert(m_node_value_list.size() == m_node_list.size());
+ Assert(m_face_value_list.size() == m_face_list.size());
}
~DirichletBoundaryCondition() = default;
@@ -183,15 +195,20 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
getMeshNodeBoundary(*mesh, dirichlet_bc_descriptor.boundaryDescriptor());
const FunctionSymbolId g_id = dirichlet_bc_descriptor.rhsSymbolId();
- // MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
- Array<const double> value_list =
+ Array<const double> node_value_list =
InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::node>(g_id, mesh->xr(),
mesh_node_boundary
.nodeList());
- boundary_condition_list.push_back(
- DirichletBoundaryCondition{mesh_face_boundary.faceList(), mesh_node_boundary.nodeList(), value_list});
-
+ Array<const double> face_value_list =
+ InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id,
+ mesh_data.xl(),
+ mesh_face_boundary
+ .faceList());
+ boundary_condition_list.push_back(DirichletBoundaryCondition{mesh_face_boundary.faceList(),
+ mesh_node_boundary.nodeList(), node_value_list,
+ face_value_list});
} else {
throw NotImplementedError("Dirichlet BC in 1d");
}
@@ -255,14 +272,20 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
getMeshNodeBoundary(*mesh, dirichlet_bc_prev_descriptor.boundaryDescriptor());
const FunctionSymbolId g_id = dirichlet_bc_prev_descriptor.rhsSymbolId();
- // MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
- Array<const double> value_list =
+ Array<const double> node_value_list =
InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::node>(g_id, mesh->xr(),
mesh_node_boundary
.nodeList());
- boundary_prev_condition_list.push_back(
- DirichletBoundaryCondition{mesh_face_boundary.faceList(), mesh_node_boundary.nodeList(), value_list});
+ Array<const double> face_value_list =
+ InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id,
+ mesh_data.xl(),
+ mesh_face_boundary
+ .faceList());
+ boundary_prev_condition_list.push_back(DirichletBoundaryCondition{mesh_face_boundary.faceList(),
+ mesh_node_boundary.nodeList(),
+ node_value_list, face_value_list});
} else {
throw NotImplementedError("Dirichlet BC in 1d");
@@ -559,13 +582,13 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
}
} else if constexpr (std::is_same_v<T, DirichletBoundaryCondition>) {
- const auto& node_list = bc.nodeList();
- const auto& value_list = bc.valueList();
+ const auto& node_list = bc.nodeList();
+ const auto& node_value_list = bc.nodeValueList();
for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
const NodeId node_id = node_list[i_node];
- compute_node_boundary_values[node_id] = value_list[i_node];
+ compute_node_boundary_values[node_id] = node_value_list[i_node];
}
}
},
@@ -596,30 +619,40 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
const auto& face_nodes = face_to_node_matrix[face_id];
-
for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
const NodeId node_id = face_nodes[i_node];
if (not node_is_dirichlet[node_id]) {
- compute_node_boundary_values[node_id] += value_list[i_face] * mes_l[face_id];
- sum_mes_l[node_id] += mes_l[face_id];
+ if (node_is_angle[node_id]) {
+ const auto& node_to_face = node_to_face_matrix[node_id];
+ for (size_t i_face_node = 0; i_face_node < node_to_face.size(); ++i_face_node) {
+ const FaceId face_node_id = node_to_face[i_face_node];
+ if (face_id == face_node_id) {
+ compute_node_boundary_values[node_id] +=
+ value_list[i_face] * mes_l[face_id] * normal_base[node_id][i_face_node];
+ }
+ }
+ } else {
+ compute_node_boundary_values[node_id] += value_list[i_face] * mes_l[face_id];
+ sum_mes_l[node_id] += mes_l[face_id];
+ }
}
}
}
} else if constexpr (std::is_same_v<T, DirichletBoundaryCondition>) {
- const auto& node_list = bc.nodeList();
- const auto& value_list = bc.valueList();
+ const auto& node_list = bc.nodeList();
+ const auto& node_value_list = bc.nodeValueList();
for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
const NodeId node_id = node_list[i_node];
- compute_node_boundary_values[node_id] = value_list[i_node];
+ compute_node_boundary_values[node_id] = node_value_list[i_node];
}
}
},
boundary_condition);
}
for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
- if ((not node_is_dirichlet[node_id]) && (node_is_neumann[node_id])) {
+ if ((not node_is_dirichlet[node_id]) && (node_is_neumann[node_id]) and not node_is_angle[node_id]) {
compute_node_boundary_values[node_id] /= sum_mes_l[node_id];
}
}
@@ -792,6 +825,8 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
return compute_sum_theta;
}();
+ // const double lambda = 0;
+
const Array<int> non_zeros{mesh->numberOfCells()};
non_zeros.fill(4); // Modif pour optimiser
CRSMatrixDescriptor<double> S1(mesh->numberOfCells(), mesh->numberOfCells(), non_zeros);