diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index dc9502f90106ee888d7d5b5f0ae7ca687b2aab66..ca49d0b16565877fc526dc9c5f7a8740975582d4 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -604,16 +604,17 @@ SchemeModule::SchemeModule()
));
this->_addBuiltinFunction("nodalheat",
- 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::vector<
- std::shared_ptr<const IBoundaryConditionDescriptor>>&,
- const std::vector<
- std::shared_ptr<const IBoundaryConditionDescriptor>>&,
- const std::shared_ptr<const IDiscreteFunction>&, const double&,
- const double&)>>(
+ 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 std::shared_ptr<const IDiscreteFunction>& kappa,
const std::shared_ptr<const IDiscreteFunction>& f,
@@ -623,7 +624,8 @@ SchemeModule::SchemeModule()
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 std::shared_ptr<const IDiscreteFunction> {
+ const double& cn_coeff) -> const std::tuple<std::shared_ptr<const IDiscreteFunction>,
+ std::shared_ptr<const IDiscreteFunction>> {
return ScalarNodalSchemeHandler{kappa,
f,
f_prev,
@@ -633,9 +635,7 @@ SchemeModule::SchemeModule()
dt,
cn_coeff}
.solution();
- }
-
- ));
+ }));
this->_addBuiltinFunction("unsteadyelasticity",
std::make_shared<BuiltinFunctionEmbedder<
diff --git a/src/scheme/ScalarNodalScheme.cpp b/src/scheme/ScalarNodalScheme.cpp
index f380a0fff001b0b20aee1d39bf4d3382e53cbf51..6a827b0145f1ee309a8bc806f064487236fb11a3 100644
--- a/src/scheme/ScalarNodalScheme.cpp
+++ b/src/scheme/ScalarNodalScheme.cpp
@@ -6,7 +6,10 @@
class ScalarNodalSchemeHandler::IScalarNodalScheme
{
public:
- virtual std::shared_ptr<const IDiscreteFunction> getSolution() const = 0;
+ 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;
IScalarNodalScheme() = default;
virtual ~IScalarNodalScheme() = default;
@@ -20,7 +23,9 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
using MeshType = Mesh<ConnectivityType>;
using MeshDataType = MeshData<Dimension>;
- std::shared_ptr<const DiscreteFunctionP0<Dimension, double>> m_solution;
+ 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
{
@@ -136,11 +141,16 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
};
public:
- std::shared_ptr<const IDiscreteFunction>
- getSolution() const final
+ std::tuple<std::shared_ptr<const IDiscreteFunction>, std::shared_ptr<const IDiscreteFunction>>
+ getSolution() const
{
- return m_solution;
+ return {m_cell_temperature, m_node_temperature};
}
+ // std::shared_ptr<const IDiscreteFunction>
+ // getSolution() const final
+ // {
+ // return m_solution;
+ // }
ScalarNodalScheme(const std::shared_ptr<const MeshType>& mesh,
const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyMatrix<Dimension>>>& cell_k,
@@ -303,6 +313,11 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
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();
@@ -313,7 +328,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
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& 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();
@@ -737,7 +752,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
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);
- CRSMatrixDescriptor<double> C(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];
@@ -747,7 +761,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
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;
- C(cell_id_j, cell_id_k) = 0;
}
}
}
@@ -835,32 +848,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
}
}
- for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
- if (not is_boundary_face[face_id]) {
- const auto& face_to_cell = face_to_cell_matrix[face_id];
- CellId cell_id_j = face_to_cell[0];
- CellId cell_id_k = face_to_cell[1];
- C(cell_id_j, cell_id_j) += dt * mes_l[face_id] / l2Norm(xj[cell_id_j] - xj[cell_id_k]);
- C(cell_id_k, cell_id_k) += dt * mes_l[face_id] / l2Norm(xj[cell_id_j] - xj[cell_id_k]);
- C(cell_id_j, cell_id_k) -= dt * mes_l[face_id] / l2Norm(xj[cell_id_j] - xj[cell_id_k]);
- C(cell_id_k, cell_id_j) -= dt * mes_l[face_id] / l2Norm(xj[cell_id_j] - xj[cell_id_k]);
- } else {
- const auto& face_to_cell = face_to_cell_matrix[face_id];
- CellId cell_id_j = face_to_cell[0];
- C(cell_id_j, cell_id_j) += dt * mes_l[face_id] / l2Norm(xj[cell_id_j] - xl[face_id]);
- }
- }
-
- const double epsilon = 0.5;
-
- for (CellId cell_id_j = 0; cell_id_j < mesh->numberOfCells(); ++cell_id_j) {
- for (CellId cell_id_k = 0; cell_id_k < mesh->numberOfCells(); ++cell_id_k) {
- if (S1(cell_id_j, cell_id_k) != 0 or C(cell_id_j, cell_id_k) != 0) {
- S1(cell_id_j, cell_id_k) = (1. - epsilon) * S1(cell_id_j, cell_id_k) + epsilon * C(cell_id_j, cell_id_k);
- }
- }
- }
-
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];
@@ -1010,15 +997,51 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
LinearSolver solver;
solver.solveLocalSystem(A1, T, B);
- m_solution = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh);
- auto& solution = *m_solution;
+ 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(
- mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { solution[cell_id] = T[cell_id]; });
+ dual_mesh->numberOfCells(),
+ PUGS_LAMBDA(CellId cell_id) { node_temperature[cell_id] = dual_node_temperature[cell_id]; });
}
}
};
-std::shared_ptr<const IDiscreteFunction>
+// std::shared_ptr<const IDiscreteFunction>
+// ScalarNodalSchemeHandler::solution() const
+// {
+// return m_scheme->getSolution();
+// }
+std::tuple<std::shared_ptr<const IDiscreteFunction>, std::shared_ptr<const IDiscreteFunction>>
ScalarNodalSchemeHandler::solution() const
{
return m_scheme->getSolution();
diff --git a/src/scheme/ScalarNodalScheme.hpp b/src/scheme/ScalarNodalScheme.hpp
index 930b15d2513a6edf8061bf9884cfc2eac1d69982..f5a0b3f7c8f2a16d4fa2a774684ec9298d1a9a8b 100644
--- a/src/scheme/ScalarNodalScheme.hpp
+++ b/src/scheme/ScalarNodalScheme.hpp
@@ -17,6 +17,7 @@
#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>
@@ -38,7 +39,8 @@ class ScalarNodalSchemeHandler
public:
std::unique_ptr<IScalarNodalScheme> m_scheme;
- std::shared_ptr<const IDiscreteFunction> solution() const;
+ // std::shared_ptr<const IDiscreteFunction> solution() const;
+ std::tuple<std::shared_ptr<const IDiscreteFunction>, std::shared_ptr<const IDiscreteFunction>> solution() const;
ScalarNodalSchemeHandler(
const std::shared_ptr<const IDiscreteFunction>& cell_k,