diff --git a/src/language/algorithms/HeatDiamondAlgorithm.cpp b/src/language/algorithms/HeatDiamondAlgorithm.cpp
index 6a6ff407608598fb0be6cb2dad283235cfa0d92e..3a2424705a4bb35dbac11dced0005bf8d2cf4056 100644
--- a/src/language/algorithms/HeatDiamondAlgorithm.cpp
+++ b/src/language/algorithms/HeatDiamondAlgorithm.cpp
@@ -1,5 +1,6 @@
#include <language/algorithms/HeatDiamondAlgorithm.hpp>
+#include <algebra/BiCGStab.hpp>
#include <algebra/CRSMatrix.hpp>
#include <algebra/LocalRectangularMatrix.hpp>
#include <algebra/PCG.hpp>
@@ -22,7 +23,8 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const FunctionSymbolId& T_id,
const FunctionSymbolId& kappa_id,
- const FunctionSymbolId& f_id)
+ const FunctionSymbolId& f_id,
+ const FunctionSymbolId& g_id)
{
using ConnectivityType = Connectivity<Dimension>;
using MeshType = Mesh<ConnectivityType>;
@@ -271,9 +273,9 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
S(i_id, j_id) += weight_rj * alpha_face_id * (nil, Clr);
}
}
- } else {
- std::cout << "*** ignore node " << face_node_id << '\n';
- }
+ } // else {
+ // std::cout << "*** ignore node " << face_node_id << '\n';
+ // }
}
}
}
@@ -283,15 +285,60 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
CellValue<double> fj =
InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(f_id, mesh_data.xj());
+ NodeValue<double> gr =
+ InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::node>(g_id, m_mesh->xr());
+
const CellValue<const double> primal_Vj = mesh_data.Vj();
CRSMatrix A{S};
Vector<double> b{m_mesh->numberOfCells()};
- parallel_for(
- m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { b[cell_id] = fj[cell_id] * primal_Vj[cell_id]; });
+ // parallel_for(
+ // m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { b[cell_id] = fj[cell_id] * primal_Vj[cell_id]; });
+
+ const auto& primal_cell_to_face_matrix = m_mesh->connectivity().cellToFaceMatrix();
+ for (CellId cell_id = 0; cell_id < m_mesh->numberOfCells(); ++cell_id) {
+ b[cell_id] = fj[cell_id] * primal_Vj[cell_id];
+ const auto& primal_cell_to_face = primal_cell_to_face_matrix[cell_id];
+
+ for (size_t face_id = 0; face_id < primal_cell_to_face.size(); ++face_id) { // Num local
+ FaceId i_cell_face = primal_cell_to_face_matrix[cell_id][face_id]; // Num global
+
+ const bool is_face_reversed_for_cell_i = [=] {
+ for (size_t i_face = 0; i_face < primal_cell_to_face.size(); ++i_face) { // Num local
+ FaceId primal_face_id = primal_cell_to_face[i_face]; // Num global
+ if (primal_face_id == i_cell_face) {
+ return primal_face_cell_is_reversed(cell_id, i_face);
+ }
+ }
+ }();
+ const auto& primal_face_to_node = primal_face_to_node_matrix[i_cell_face];
+ for (size_t i_node = 0; i_node < primal_face_to_node.size(); ++i_node) { // Num local
+ const TinyVector<Dimension> nil = [=] {
+ if (is_face_reversed_for_cell_i) {
+ return -primal_nlr(i_cell_face, i_node);
+ } else {
+ return primal_nlr(i_cell_face, i_node);
+ }
+ }();
+
+ CellId dual_cell_id = face_dual_cell_id[i_cell_face]; // Num global dual de la face
+
+ NodeId face_node_id = primal_face_to_node[i_node]; // Num global
+ if (primal_node_is_on_boundary[face_node_id]) {
+ for (size_t i_dual_node = 0; i_dual_node < dual_Cjr.numberOfSubValues(dual_cell_id); ++i_dual_node) {
+ const NodeId dual_node_id = dual_cell_to_node_matrix[dual_cell_id][i_dual_node];
+ if (dual_node_id == face_node_id) {
+ const TinyVector<Dimension> Clr = dual_Cjr(dual_cell_id, i_dual_node);
+ b[cell_id] -= alpha_l[i_cell_face] * gr[face_node_id] * (nil, Clr);
+ }
+ }
+ }
+ }
+ }
+ }
Vector<double> T{m_mesh->numberOfCells()};
T = 0;
- PCG{b, A, A, T, 100, 1e-12};
+ BiCGStab{b, A, T, 1000, 1e-12};
CellValue<double> Temperature{m_mesh->connectivity()};
@@ -322,6 +369,7 @@ template HeatDiamondScheme<1>::HeatDiamondScheme(
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
const FunctionSymbolId&,
const FunctionSymbolId&,
+ const FunctionSymbolId&,
const FunctionSymbolId&);
template HeatDiamondScheme<2>::HeatDiamondScheme(
@@ -329,6 +377,7 @@ template HeatDiamondScheme<2>::HeatDiamondScheme(
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
const FunctionSymbolId&,
const FunctionSymbolId&,
+ const FunctionSymbolId&,
const FunctionSymbolId&);
template HeatDiamondScheme<3>::HeatDiamondScheme(
@@ -336,4 +385,5 @@ template HeatDiamondScheme<3>::HeatDiamondScheme(
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
const FunctionSymbolId&,
const FunctionSymbolId&,
+ const FunctionSymbolId&,
const FunctionSymbolId&);
diff --git a/src/language/algorithms/HeatDiamondAlgorithm.hpp b/src/language/algorithms/HeatDiamondAlgorithm.hpp
index b1753895140fb2d655a530973bd22f934a187c42..fc8885d618f4aa5071ab74e640cad689e9314d71 100644
--- a/src/language/algorithms/HeatDiamondAlgorithm.hpp
+++ b/src/language/algorithms/HeatDiamondAlgorithm.hpp
@@ -15,7 +15,8 @@ struct HeatDiamondScheme
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const FunctionSymbolId& T_id,
const FunctionSymbolId& kappa_id,
- const FunctionSymbolId& f_id);
+ const FunctionSymbolId& f_id,
+ const FunctionSymbolId& g_id);
};
#endif // HEAT_DIAMOND_ALGORITHM_HPP
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 64c8fa03b4bd395b3a6d3f985467f506afef0093..cd5d2e42a21e6728d19ca02505d0f23c501d6541 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -176,24 +176,25 @@ SchemeModule::SchemeModule()
this->_addBuiltinFunction("heat", std::make_shared<BuiltinFunctionEmbedder<
void(std::shared_ptr<const IMesh>,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
- const FunctionSymbolId&, const FunctionSymbolId&, const FunctionSymbolId&)>>(
+ const FunctionSymbolId&, const FunctionSymbolId&, const FunctionSymbolId&,
+ const FunctionSymbolId&)>>(
[](std::shared_ptr<const IMesh> p_mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list,
const FunctionSymbolId& T_id, const FunctionSymbolId& kappa_id,
- const FunctionSymbolId& f_id) -> void {
+ const FunctionSymbolId& f_id, const FunctionSymbolId& g_id) -> void {
switch (p_mesh->dimension()) {
case 1: {
- HeatDiamondScheme<1>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
+ HeatDiamondScheme<1>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id, g_id};
break;
}
case 2: {
- HeatDiamondScheme<2>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
+ HeatDiamondScheme<2>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id, g_id};
break;
}
case 3: {
- HeatDiamondScheme<3>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
+ HeatDiamondScheme<3>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id, g_id};
break;
}
default: {