diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index abb60970a238402f360412b380ea05092f360740..dc9502f90106ee888d7d5b5f0ae7ca687b2aab66 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -603,25 +603,39 @@ 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::vector<
- std::shared_ptr<const IBoundaryConditionDescriptor>>&,
- const std::shared_ptr<const IDiscreteFunction>&, const double&)>>(
-
- [](const std::shared_ptr<const IDiscreteFunction>& kappa,
- const std::shared_ptr<const IDiscreteFunction>& f,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
- bc_descriptor_list,
- const std::shared_ptr<const IDiscreteFunction>& P,
- const double& dt) -> const std::shared_ptr<const IDiscreteFunction> {
- return ScalarNodalSchemeHandler{kappa, f, bc_descriptor_list, P, dt}.solution();
- }
+ 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&)>>(
+
+ [](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 std::shared_ptr<const IDiscreteFunction> {
+ return ScalarNodalSchemeHandler{kappa,
+ f,
+ f_prev,
+ bc_descriptor_list,
+ bc_prev_descriptor_list,
+ P,
+ dt,
+ cn_coeff}
+ .solution();
+ }
- ));
+ ));
this->_addBuiltinFunction("unsteadyelasticity",
std::make_shared<BuiltinFunctionEmbedder<
diff --git a/src/scheme/ScalarNodalScheme.cpp b/src/scheme/ScalarNodalScheme.cpp
index 7fd220f64f08b536255d45ec89fbb6bfa46fbbd0..06fe956c09ea9e03fddb41a9abf6e5c170e5f4e6 100644
--- a/src/scheme/ScalarNodalScheme.cpp
+++ b/src/scheme/ScalarNodalScheme.cpp
@@ -145,9 +145,12 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
ScalarNodalScheme(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& dt,
+ const double& cn_coeff)
{
using BoundaryCondition =
std::variant<DirichletBoundaryCondition, FourierBoundaryCondition, NeumannBoundaryCondition>;
@@ -227,6 +230,79 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
}
}
+ 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> 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(), 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);
const NodeValue<const TinyVector<Dimension>>& xr = mesh->xr();
@@ -408,11 +484,70 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
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]) {
+ 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& 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_neumann[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];
+ } else {
+ compute_node_boundary_values[node_id] = value_list[i_face];
+ }
+ }
+ }
+ }
+ },
+ 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])) {
+ compute_node_boundary_values[node_id] /= sum_mes_l[node_id];
+ } else if ((not node_is_neumann[node_id]) && (node_is_dirichlet[node_id])) {
+ compute_node_boundary_values[node_id] /= sum_mes_l[node_id];
+ }
+ }
+ return compute_node_boundary_values;
+ }();
+
{
CellValue<const TinyMatrix<Dimension>> cell_kappaj = cell_k->cellValues();
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;
@@ -431,6 +566,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
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);
@@ -511,6 +647,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
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];
@@ -554,6 +691,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
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])) {
@@ -573,6 +711,16 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
const Array<int> non_zeros{mesh->numberOfCells()};
non_zeros.fill(Dimension); // Modif pour optimiser
CRSMatrixDescriptor<double> S(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];
+ S(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];
@@ -726,6 +874,8 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
parallel_for(
mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { B[cell_id] = dt * b[cell_id] + Ph[cell_id]; });
+ // std::cout << "g = " << node_boundary_values << "\n";
+
LinearSolver solver;
solver.solveLocalSystem(A, T, B);
@@ -746,9 +896,12 @@ ScalarNodalSchemeHandler::solution() const
ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
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& dt,
+ const double& cn_coeff)
{
const std::shared_ptr i_mesh = getCommonMesh({cell_k, f});
if (not i_mesh) {
@@ -772,8 +925,10 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
m_scheme =
std::make_unique<ScalarNodalScheme<2>>(mesh, std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k),
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
- bc_descriptor_list,
- std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(P), dt);
+ 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);
break;
}
case 3: {
diff --git a/src/scheme/ScalarNodalScheme.hpp b/src/scheme/ScalarNodalScheme.hpp
index a43ff3b4ca7b46b45948faeae252bdff62caa450..930b15d2513a6edf8061bf9884cfc2eac1d69982 100644
--- a/src/scheme/ScalarNodalScheme.hpp
+++ b/src/scheme/ScalarNodalScheme.hpp
@@ -40,11 +40,15 @@ class ScalarNodalSchemeHandler
std::shared_ptr<const IDiscreteFunction> solution() const;
- ScalarNodalSchemeHandler(const std::shared_ptr<const IDiscreteFunction>& cell_k,
- const std::shared_ptr<const IDiscreteFunction>& f,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
- const std::shared_ptr<const IDiscreteFunction>& P,
- const double& dt);
+ ScalarNodalSchemeHandler(
+ 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);
~ScalarNodalSchemeHandler();
};