diff --git a/src/algebra/SparseMatrixDescriptor.hpp b/src/algebra/SparseMatrixDescriptor.hpp
index d5c14713bb6b91235b796a56f40b9a2f5a605805..e9a8899351149734f7a65cb21b60e7fba69287ca 100644
--- a/src/algebra/SparseMatrixDescriptor.hpp
+++ b/src/algebra/SparseMatrixDescriptor.hpp
@@ -34,16 +34,14 @@ class SparseMatrixDescriptor
return m_id_value_map.size();
}
- const DataType&
- operator[](const IndexType& j) const
+ const DataType& operator[](const IndexType& j) const
{
auto i_value = m_id_value_map.find(j);
Assert(i_value != m_id_value_map.end());
return i_value->second;
}
- DataType&
- operator[](const IndexType& j)
+ DataType& operator[](const IndexType& j)
{
auto i_value = m_id_value_map.find(j);
if (i_value != m_id_value_map.end()) {
@@ -100,6 +98,18 @@ class SparseMatrixDescriptor
return m_row_array[i][j];
}
+ SparseMatrixDescriptor&
+ operator*=(const double& lambda)
+ {
+ for (IndexType i_row = 0; i_row < m_row_array.size(); ++i_row) {
+ SparseRowDescriptor& row = m_row_array[i_row];
+ for (auto& [id, value] : row.m_id_value_map) {
+ value *= lambda;
+ }
+ }
+ return *this;
+ }
+
const DataType&
operator()(const IndexType& i, const IndexType& j) const
{
diff --git a/src/language/algorithms/ParabolicHeat.cpp b/src/language/algorithms/ParabolicHeat.cpp
index a3813dd47ae4ead9804ff668f793681529cd4138..5dcfb1bfeb186fab3693ffb32697f755c067aa4b 100644
--- a/src/language/algorithms/ParabolicHeat.cpp
+++ b/src/language/algorithms/ParabolicHeat.cpp
@@ -31,7 +31,7 @@ ParabolicHeatScheme<Dimension>::ParabolicHeatScheme(
const FunctionSymbolId& T_id,
const FunctionSymbolId& kappa_id,
const FunctionSymbolId& f_id,
- const double& final_time,
+ const double& Tf,
const double& dt)
{
using ConnectivityType = Connectivity<Dimension>;
@@ -252,7 +252,15 @@ ParabolicHeatScheme<Dimension>::ParabolicHeatScheme(
NamedItemValue{"node_owner", mesh->connectivity().nodeOwner()}},
0, true); // forces last output
}
- ScalarDiamondScheme<Dimension>(i_mesh, bc_descriptor_list, kappa_id, f_id, Temperature, Temperature_face);
+ double time = 0;
+ Assert(dt > 0, "The time-step have to be positive!");
+ do {
+ double deltat = std::min(dt, Tf - time);
+ std::cout << "Current time = " << time << " time-step = " << deltat << " final time = " << Tf << "\n";
+ ScalarDiamondScheme<Dimension>(i_mesh, bc_descriptor_list, kappa_id, f_id, Temperature, Temperature_face, Tf,
+ deltat);
+ time += deltat;
+ } while (time < Tf && std::abs(time - Tf) > 1e-15);
{
Vector<double> error{mesh->numberOfCells()};
CellValue<double> cell_error{mesh->connectivity()};
diff --git a/src/scheme/ScalarDiamondScheme.hpp b/src/scheme/ScalarDiamondScheme.hpp
index 704b4c21e6ea4747534d0be558cbd63f901a8dc4..87e3ae3205d621d252155b602445da525e2aedf8 100644
--- a/src/scheme/ScalarDiamondScheme.hpp
+++ b/src/scheme/ScalarDiamondScheme.hpp
@@ -150,7 +150,9 @@ class ScalarDiamondScheme
const FunctionSymbolId& kappa_id,
const FunctionSymbolId& f_id,
CellValue<double>& Temperature,
- FaceValue<double>& Temperature_face)
+ FaceValue<double>& Temperature_face,
+ const double& Tf,
+ const double& dt)
{
using ConnectivityType = Connectivity<Dimension>;
using MeshType = Mesh<ConnectivityType>;
@@ -583,6 +585,11 @@ class ScalarDiamondScheme
return computed_alpha_l;
}();
+ double lambda = (Tf == 0) ? 0 : 1;
+ double time_factor = (lambda == 0) ? 1 : dt;
+ std::cout << " lambda = " << lambda << " time_factor =" << time_factor << "\n";
+ const CellValue<const double> primal_Vj = mesh_data.Vj();
+
SparseMatrixDescriptor S{number_of_dof};
for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
@@ -593,12 +600,13 @@ class ScalarDiamondScheme
for (size_t j_cell = 0; j_cell < primal_face_to_cell.size(); ++j_cell) {
const CellId cell_id2 = primal_face_to_cell[j_cell];
if (i_cell == j_cell) {
- S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) += beta_l;
+ S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) +=
+ (time_factor * beta_l + lambda * primal_Vj[cell_id1]);
if (primal_face_is_neumann[face_id]) {
S(face_dof_number[face_id], cell_dof_number[cell_id2]) -= beta_l;
}
} else {
- S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) -= beta_l;
+ S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) -= time_factor * beta_l;
}
}
}
@@ -636,7 +644,7 @@ class ScalarDiamondScheme
for (size_t j_cell = 0; j_cell < primal_node_to_cell_matrix[node_id].size(); ++j_cell) {
CellId j_id = primal_node_to_cell_matrix[node_id][j_cell];
- S(cell_dof_number[i_id], cell_dof_number[j_id]) -= w_rj(node_id, j_cell) * a_ir;
+ S(cell_dof_number[i_id], cell_dof_number[j_id]) -= time_factor * w_rj(node_id, j_cell) * a_ir;
if (primal_face_is_neumann[face_id]) {
S(face_dof_number[face_id], cell_dof_number[j_id]) += w_rj(node_id, j_cell) * a_ir;
}
@@ -645,7 +653,7 @@ class ScalarDiamondScheme
for (size_t l_face = 0; l_face < node_to_face_matrix[node_id].size(); ++l_face) {
FaceId l_id = node_to_face_matrix[node_id][l_face];
if (primal_face_is_on_boundary[l_id]) {
- S(cell_dof_number[i_id], face_dof_number[l_id]) -= w_rl(node_id, l_face) * a_ir;
+ S(cell_dof_number[i_id], face_dof_number[l_id]) -= time_factor * w_rl(node_id, l_face) * a_ir;
if (primal_face_is_neumann[face_id]) {
S(face_dof_number[face_id], face_dof_number[l_id]) += w_rl(node_id, l_face) * a_ir;
}
@@ -667,12 +675,12 @@ class ScalarDiamondScheme
InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(f_id,
mesh_data.xj());
- const CellValue<const double> primal_Vj = mesh_data.Vj();
CRSMatrix A{S};
Vector<double> b{number_of_dof};
for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
- b[cell_dof_number[cell_id]] = fj[cell_id] * primal_Vj[cell_id];
+ b[cell_dof_number[cell_id]] =
+ (time_factor * fj[cell_id] + lambda * Temperature[cell_id]) * primal_Vj[cell_id];
}
// Dirichlet on b^L_D
{
@@ -703,7 +711,8 @@ class ScalarDiamondScheme
const auto& value_list = bc.valueList();
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
FaceId face_id = face_list[i_face];
- b[face_dof_number[face_id]] += mes_l[face_id] * value_list[i_face]; // sign
+ b[face_dof_number[face_id]] +=
+ mes_l[face_id] * value_list[i_face]; // + lambda * Temperature_face[face_id]); // sign
}
}
},
@@ -735,7 +744,9 @@ template ScalarDiamondScheme<1>::ScalarDiamondScheme(
const FunctionSymbolId&,
const FunctionSymbolId&,
CellValue<double>&,
- FaceValue<double>&);
+ FaceValue<double>&,
+ const double&,
+ const double&);
template ScalarDiamondScheme<2>::ScalarDiamondScheme(
std::shared_ptr<const IMesh>,
@@ -743,7 +754,9 @@ template ScalarDiamondScheme<2>::ScalarDiamondScheme(
const FunctionSymbolId&,
const FunctionSymbolId&,
CellValue<double>&,
- FaceValue<double>&);
+ FaceValue<double>&,
+ const double&,
+ const double&);
template ScalarDiamondScheme<3>::ScalarDiamondScheme(
std::shared_ptr<const IMesh>,
@@ -751,4 +764,6 @@ template ScalarDiamondScheme<3>::ScalarDiamondScheme(
const FunctionSymbolId&,
const FunctionSymbolId&,
CellValue<double>&,
- FaceValue<double>&);
+ FaceValue<double>&,
+ const double&,
+ const double&);
diff --git a/tests/test_SparseMatrixDescriptor.cpp b/tests/test_SparseMatrixDescriptor.cpp
index 8616a41da0ff8a60857ef314c74e04becf6c651f..c18fd804df81f310f1d9d49008a593dc910181cf 100644
--- a/tests/test_SparseMatrixDescriptor.cpp
+++ b/tests/test_SparseMatrixDescriptor.cpp
@@ -196,4 +196,31 @@ TEST_CASE("SparseMatrixDescriptor", "[algebra]")
)";
REQUIRE(output.str() == expected_output);
}
+ SECTION("Multiplication by scalar")
+ {
+ SparseMatrixDescriptor<int, uint8_t> S{5};
+ S(0, 2) = 3;
+ S(1, 2) = 11;
+ S(1, 1) = 1;
+ S(0, 2) += 2;
+ S(3, 3) = 5;
+ S(3, 1) = -3;
+ S(4, 1) = 1;
+ S(2, 2) = 4;
+ S(4, 4) = -2;
+ S *= 3;
+ const auto value_array = S.valueArray();
+
+ REQUIRE(value_array.size() == 9);
+
+ REQUIRE(value_array[0] == 0);
+ REQUIRE(value_array[1] == 15);
+ REQUIRE(value_array[2] == 3);
+ REQUIRE(value_array[3] == 33);
+ REQUIRE(value_array[4] == 12);
+ REQUIRE(value_array[5] == -9);
+ REQUIRE(value_array[6] == 15);
+ REQUIRE(value_array[7] == 3);
+ REQUIRE(value_array[8] == -6);
+ }
}