diff --git a/src/language/algorithms/CMakeLists.txt b/src/language/algorithms/CMakeLists.txt
index 45c9a5c489271b68e718bace2552b2e02ec56793..e8bdb537aeac576b769a107f8d1d751c4d75b957 100644
--- a/src/language/algorithms/CMakeLists.txt
+++ b/src/language/algorithms/CMakeLists.txt
@@ -3,6 +3,7 @@
add_library(PugsLanguageAlgorithms
AcousticSolverAlgorithm.cpp
Heat5PointsAlgorithm.cpp
+ HeatDiamondAlgorithm.cpp
)
diff --git a/src/language/algorithms/HeatDiamondAlgorithm.cpp b/src/language/algorithms/HeatDiamondAlgorithm.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4c2e5d69257abfe17e349ca05f1b695c0ea79f1d
--- /dev/null
+++ b/src/language/algorithms/HeatDiamondAlgorithm.cpp
@@ -0,0 +1,251 @@
+#include <language/algorithms/HeatDiamondAlgorithm.hpp>
+
+#include <algebra/CRSMatrix.hpp>
+#include <algebra/LocalRectangularMatrix.hpp>
+#include <algebra/PCG.hpp>
+#include <algebra/SparseMatrixDescriptor.hpp>
+#include <algebra/TinyVector.hpp>
+#include <algebra/Vector.hpp>
+#include <language/utils/InterpolateItemValue.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/ConnectivityToDiamondDualConnectivityDataMapper.hpp>
+#include <mesh/DiamondDualConnectivityManager.hpp>
+#include <mesh/DiamondDualMeshManager.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <output/VTKWriter.hpp>
+
+template <size_t Dimension>
+HeatDiamondScheme<Dimension>::HeatDiamondScheme(
+ std::shared_ptr<const IMesh> i_mesh,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ const FunctionSymbolId& T_id,
+ const FunctionSymbolId& kappa_id)
+{
+ using ConnectivityType = Connectivity<Dimension>;
+ using MeshType = Mesh<ConnectivityType>;
+ using MeshDataType = MeshData<Dimension>;
+
+ std::cout << "number of bc descr = " << bc_descriptor_list.size() << '\n';
+ std::shared_ptr m_mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+
+ if constexpr (Dimension > 1) {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
+
+ CellValue<double> Tj =
+ InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(T_id, mesh_data.xj());
+
+ NodeValue<double> Tr(m_mesh->connectivity());
+ const NodeValue<const TinyVector<Dimension>>& xr = m_mesh->xr();
+ const CellValue<const TinyVector<Dimension>>& xj = mesh_data.xj();
+ const auto& node_to_cell_matrix = m_mesh->connectivity().nodeToCellMatrix();
+ CellValuePerNode<double> w_rj{m_mesh->connectivity()};
+
+ for (NodeId i_node = 0; i_node < m_mesh->numberOfNodes(); ++i_node) {
+ Vector<double> b{Dimension + 1};
+ b[0] = 1;
+ for (size_t i = 1; i < Dimension + 1; i++) {
+ b[i] = xr[i_node][i - 1];
+ }
+ const auto& node_to_cell = node_to_cell_matrix[i_node];
+
+ LocalRectangularMatrix<double> A{Dimension + 1, node_to_cell.size()};
+ for (size_t j = 0; j < node_to_cell.size(); j++) {
+ A(0, j) = 1;
+ }
+ for (size_t i = 1; i < Dimension + 1; i++) {
+ for (size_t j = 0; j < node_to_cell.size(); j++) {
+ const CellId J = node_to_cell[j];
+ A(i, j) = xj[J][i - 1];
+ }
+ }
+ Vector<double> x{node_to_cell.size()};
+ x = 0;
+
+ LocalRectangularMatrix B = transpose(A) * A;
+ Vector y = transpose(A) * b;
+ PCG{y, B, B, x, 10, 1e-12};
+
+ Tr[i_node] = 0;
+ for (size_t j = 0; j < node_to_cell.size(); j++) {
+ Tr[i_node] += x[j] * Tj[node_to_cell[j]];
+ w_rj(i_node, j) = x[j];
+ }
+ }
+
+ {
+ VTKWriter vtk_writer("T_" + std::to_string(Dimension), 0.01);
+
+ vtk_writer.write(m_mesh,
+ {NamedItemValue{"Tr", Tr}, NamedItemValue{"temperature", Tj},
+ NamedItemValue{"coords", m_mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()},
+ NamedItemValue{"cell_owner", m_mesh->connectivity().cellOwner()},
+ NamedItemValue{"node_owner", m_mesh->connectivity().nodeOwner()}},
+ 0, true); // forces last output
+ }
+
+ {
+ std::shared_ptr diamond_mesh = DiamondDualMeshManager::instance().getDiamondDualMesh(m_mesh);
+
+ MeshDataType& diamond_mesh_data = MeshDataManager::instance().getMeshData(*diamond_mesh);
+
+ std::shared_ptr mapper =
+ DiamondDualConnectivityManager::instance().getConnectivityToDiamondDualConnectivityDataMapper(
+ m_mesh->connectivity());
+
+ NodeValue<double> Trd{diamond_mesh->connectivity()};
+
+ mapper->toDualNode(Tr, Tj, Trd);
+
+ CellValue<double> kappaj =
+ InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(kappa_id,
+ mesh_data.xj());
+
+ CellValue<double> dual_kappaj =
+ InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(kappa_id,
+ diamond_mesh_data
+ .xj());
+
+ VTKWriter vtk_writer("D_" + std::to_string(Dimension), 0.01);
+
+ vtk_writer.write(diamond_mesh,
+ {NamedItemValue{"kappaj", dual_kappaj}, NamedItemValue{"coords", diamond_mesh->xr()},
+ NamedItemValue{"Vj", diamond_mesh_data.Vj()}, NamedItemValue{"xj", diamond_mesh_data.xj()}},
+ 0, true); // forces last output
+
+ const CellValue<const double> dual_Vj = diamond_mesh_data.Vj();
+ const auto& face_to_node_matrix = m_mesh->connectivity().faceToNodeMatrix();
+
+ const FaceValue<const double> mes_l = [=] {
+ FaceValue<double> compute_mes_l{m_mesh->connectivity()};
+ const NodeValue<const TinyVector<Dimension>>& xr = m_mesh->xr();
+ if constexpr (Dimension == 1) {
+ compute_mes_l.fill(1);
+ } else if constexpr (Dimension == 2) {
+ parallel_for(
+ m_mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
+ const auto& face_to_node = face_to_node_matrix[face_id];
+ const NodeId node_id1 = face_to_node[0];
+ const NodeId node_id2 = face_to_node[1];
+ const TinyVector<Dimension, double> r = xr[node_id1] - xr[node_id2];
+ compute_mes_l[face_id] = l2Norm(r);
+ });
+ } else {
+ throw NotImplementedError("Not implemented in 3D");
+ }
+ return compute_mes_l;
+ }();
+
+ const CellValue<const double> dual_mes_l_j = [=] {
+ CellValue<double> compute_mes_j{diamond_mesh->connectivity()};
+ mapper->toDualCell(mes_l, compute_mes_j);
+
+ return compute_mes_j;
+ }();
+
+ FaceValue<const double> alpha_l = [&] {
+ CellValue<double> alpha_j{diamond_mesh->connectivity()};
+
+ parallel_for(
+ diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId diamond_cell_id) {
+ alpha_j[diamond_cell_id] =
+ dual_mes_l_j[diamond_cell_id] * dual_kappaj[diamond_cell_id] / dual_Vj[diamond_cell_id];
+ });
+
+ FaceValue<double> computed_alpha_l{m_mesh->connectivity()};
+ mapper->fromDualCell(alpha_j, computed_alpha_l);
+
+ VTKWriter vtk_writer("DD_" + std::to_string(Dimension), 0.01);
+
+ vtk_writer.write(diamond_mesh,
+ {NamedItemValue{"alphaj", alpha_j}, NamedItemValue{"xj", diamond_mesh_data.xj()},
+ NamedItemValue{"Vl", diamond_mesh_data.Vj()}, NamedItemValue{"|l|", dual_mes_l_j}},
+ 0,
+ true); // forces last output
+
+ return computed_alpha_l;
+ }();
+
+ SparseMatrixDescriptor S{m_mesh->numberOfCells()};
+ const auto& face_to_cell_matrix = m_mesh->connectivity().faceToCellMatrix();
+ for (FaceId face_id = 0; face_id < m_mesh->numberOfFaces(); ++face_id) {
+ const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
+ const double beta_l = 0.5 * alpha_l[face_id] * mes_l[face_id];
+ const NodeValuePerCell<const TinyVector<Dimension>>& Cjr = mesh_data.Cjr();
+
+ // CellValue<double> dual_cell_id{diamond_mesh->connectivity()};
+ // mapper->toDualCell(face_id, dual_cell_id);
+
+ for (size_t i_cell = 0; i_cell < primal_face_to_cell.size(); ++i_cell) {
+ const CellId cell_id1 = primal_face_to_cell[i_cell];
+ 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_id1, cell_id2) -= beta_l;
+ } else {
+ S(cell_id1, cell_id2) += beta_l;
+ }
+ // S(cell_id1, cell_id2)+= alpha_l[face_id]*(,Cjr[dual_cell_id,dual_node_j]);
+ }
+ }
+ }
+ const CellValue<const double> primal_Vj = mesh_data.Vj();
+ CRSMatrix A{S};
+ Vector<double> b{m_mesh->numberOfCells()};
+ double pi = 3.14159265;
+ for (CellId i_cell = 0; i_cell < m_mesh->numberOfCells(); ++i_cell) {
+ if constexpr (Dimension == 2) {
+ b[i_cell] =
+ -2 * kappaj[i_cell] * pi * pi * sin(pi * xj[i_cell][0]) * sin(pi * xj[i_cell][1]) * primal_Vj[i_cell];
+
+ } else if constexpr (Dimension == 3) {
+ b[i_cell] = -3 * kappaj[i_cell] * pi * pi * sin(pi * xj[i_cell][0]) * sin(pi * xj[i_cell][1]) *
+ sin(pi * xj[i_cell][2]) * primal_Vj[i_cell];
+ }
+ }
+
+ Vector<double> T{m_mesh->numberOfCells()};
+ T = 0;
+ PCG{b, A, A, T, 100, 1e-12};
+
+ CellValue<double> Temperature{m_mesh->connectivity()};
+
+ parallel_for(
+ m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { Temperature[cell_id] = T[cell_id]; });
+
+ Vector<double> error{m_mesh->numberOfCells()};
+ parallel_for(
+ m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { error[cell_id] = Temperature[cell_id] - Tj[cell_id]; });
+
+ std::cout << "||Error||_2 = " << std::sqrt((error, error)) << "\n";
+
+ {
+ VTKWriter vtk_writer("Temperature_" + std::to_string(Dimension), 0.01);
+
+ vtk_writer.write(m_mesh, {NamedItemValue{"T", Temperature}}, 0,
+ true); // forces last output
+ }
+ }
+ } else {
+ throw NotImplementedError("not done in 3d");
+ }
+}
+
+template HeatDiamondScheme<1>::HeatDiamondScheme(
+ std::shared_ptr<const IMesh>,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+ const FunctionSymbolId&,
+ const FunctionSymbolId&);
+
+template HeatDiamondScheme<2>::HeatDiamondScheme(
+ std::shared_ptr<const IMesh>,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+ const FunctionSymbolId&,
+ const FunctionSymbolId&);
+
+template HeatDiamondScheme<3>::HeatDiamondScheme(
+ std::shared_ptr<const IMesh>,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+ const FunctionSymbolId&,
+ const FunctionSymbolId&);
diff --git a/src/language/algorithms/HeatDiamondAlgorithm.hpp b/src/language/algorithms/HeatDiamondAlgorithm.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..9516b8277c8bacdb3550133ca4a15fcf062309c1
--- /dev/null
+++ b/src/language/algorithms/HeatDiamondAlgorithm.hpp
@@ -0,0 +1,20 @@
+#ifndef HEAT_DIAMOND_ALGORITHM_HPP
+#define HEAT_DIAMOND_ALGORITHM_HPP
+
+#include <language/utils/FunctionSymbolId.hpp>
+#include <mesh/IMesh.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+
+#include <memory>
+#include <vector>
+
+template <size_t Dimension>
+struct HeatDiamondScheme
+{
+ HeatDiamondScheme(std::shared_ptr<const IMesh> i_mesh,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ const FunctionSymbolId& T_id,
+ const FunctionSymbolId& kappa_id);
+};
+
+#endif // HEAT_DIAMOND_ALGORITHM_HPP
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 6e7bd22fe6490e180bb275c79c2f976e478d68c1..af11321e2b0111aa2e0bff7b076aaa5b003b3636 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -1,256 +1,20 @@
#include <language/modules/SchemeModule.hpp>
-#include <algebra/CRSMatrix.hpp>
-#include <algebra/LocalRectangularMatrix.hpp>
-#include <algebra/PCG.hpp>
-#include <algebra/SparseMatrixDescriptor.hpp>
-#include <algebra/TinyVector.hpp>
-#include <algebra/Vector.hpp>
#include <language/algorithms/AcousticSolverAlgorithm.hpp>
#include <language/algorithms/Heat5PointsAlgorithm.hpp>
+#include <language/algorithms/HeatDiamondAlgorithm.hpp>
#include <language/utils/BuiltinFunctionEmbedder.hpp>
-#include <language/utils/PugsFunctionAdapter.hpp>
#include <language/utils/TypeDescriptor.hpp>
-#include <mesh/ConnectivityToDiamondDualConnectivityDataMapper.hpp>
-#include <mesh/DiamondDualConnectivityManager.hpp>
-#include <mesh/DiamondDualMeshBuilder.hpp>
-#include <mesh/DiamondDualMeshManager.hpp>
-#include <mesh/Mesh.hpp>
-#include <mesh/MeshData.hpp>
-#include <mesh/SubItemValuePerItem.hpp>
-#include <output/VTKWriter.hpp>
-#include <scheme/AcousticSolver.hpp>
+#include <mesh/IMesh.hpp>
#include <scheme/IBoundaryConditionDescriptor.hpp>
#include <scheme/IBoundaryDescriptor.hpp>
#include <scheme/NamedBoundaryDescriptor.hpp>
#include <scheme/NumberedBoundaryDescriptor.hpp>
-
-#include <memory>
-
+#include <scheme/PressureBoundaryConditionDescriptor.hpp>
#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
-/////////// TEMPORARY
-
-template <size_t Dimension>
-struct HeatDiamondScheme
-{
- using ConnectivityType = Connectivity<Dimension>;
- using MeshType = Mesh<ConnectivityType>;
- using MeshDataType = MeshData<Dimension>;
- using UnknownsType = FiniteVolumesEulerUnknowns<MeshType>;
-
- std::shared_ptr<const MeshType> m_mesh;
-
- HeatDiamondScheme(std::shared_ptr<const IMesh> i_mesh,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
- const FunctionSymbolId& T_id,
- const FunctionSymbolId& kappa_id)
- : m_mesh{std::dynamic_pointer_cast<const MeshType>(i_mesh)}
- {
- std::cout << "number of bc descr = " << bc_descriptor_list.size() << '\n';
-
- if constexpr (Dimension > 1) {
- MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
-
- CellValue<double> Tj =
- InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(T_id, mesh_data.xj());
-
- NodeValue<double> Tr(m_mesh->connectivity());
- const NodeValue<const TinyVector<Dimension>>& xr = m_mesh->xr();
- const CellValue<const TinyVector<Dimension>>& xj = mesh_data.xj();
- const auto& node_to_cell_matrix = m_mesh->connectivity().nodeToCellMatrix();
- CellValuePerNode<double> w_rj{m_mesh->connectivity()};
-
- for (NodeId i_node = 0; i_node < m_mesh->numberOfNodes(); ++i_node) {
- Vector<double> b{Dimension + 1};
- b[0] = 1;
- for (size_t i = 1; i < Dimension + 1; i++) {
- b[i] = xr[i_node][i - 1];
- }
- const auto& node_to_cell = node_to_cell_matrix[i_node];
-
- LocalRectangularMatrix<double> A{Dimension + 1, node_to_cell.size()};
- for (size_t j = 0; j < node_to_cell.size(); j++) {
- A(0, j) = 1;
- }
- for (size_t i = 1; i < Dimension + 1; i++) {
- for (size_t j = 0; j < node_to_cell.size(); j++) {
- const CellId J = node_to_cell[j];
- A(i, j) = xj[J][i - 1];
- }
- }
- Vector<double> x{node_to_cell.size()};
- x = 0;
-
- LocalRectangularMatrix B = transpose(A) * A;
- Vector y = transpose(A) * b;
- PCG{y, B, B, x, 10, 1e-12};
-
- Tr[i_node] = 0;
- for (size_t j = 0; j < node_to_cell.size(); j++) {
- Tr[i_node] += x[j] * Tj[node_to_cell[j]];
- w_rj(i_node, j) = x[j];
- }
- }
-
- {
- VTKWriter vtk_writer("T_" + std::to_string(Dimension), 0.01);
-
- vtk_writer.write(m_mesh,
- {NamedItemValue{"Tr", Tr}, NamedItemValue{"temperature", Tj},
- NamedItemValue{"coords", m_mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()},
- NamedItemValue{"cell_owner", m_mesh->connectivity().cellOwner()},
- NamedItemValue{"node_owner", m_mesh->connectivity().nodeOwner()}},
- 0, true); // forces last output
- }
-
- {
- std::shared_ptr diamond_mesh = DiamondDualMeshManager::instance().getDiamondDualMesh(m_mesh);
-
- MeshDataType& diamond_mesh_data = MeshDataManager::instance().getMeshData(*diamond_mesh);
-
- std::shared_ptr mapper =
- DiamondDualConnectivityManager::instance().getConnectivityToDiamondDualConnectivityDataMapper(
- m_mesh->connectivity());
-
- NodeValue<double> Trd{diamond_mesh->connectivity()};
-
- mapper->toDualNode(Tr, Tj, Trd);
-
- CellValue<double> kappaj =
- InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(kappa_id,
- mesh_data.xj());
-
- CellValue<double> dual_kappaj =
- InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(kappa_id,
- diamond_mesh_data
- .xj());
+#include <scheme/VelocityBoundaryConditionDescriptor.hpp>
- VTKWriter vtk_writer("D_" + std::to_string(Dimension), 0.01);
-
- vtk_writer.write(diamond_mesh,
- {NamedItemValue{"kappaj", dual_kappaj}, NamedItemValue{"coords", diamond_mesh->xr()},
- NamedItemValue{"Vj", diamond_mesh_data.Vj()}, NamedItemValue{"xj", diamond_mesh_data.xj()}},
- 0, true); // forces last output
-
- const CellValue<const double> dual_Vj = diamond_mesh_data.Vj();
- const auto& face_to_node_matrix = m_mesh->connectivity().faceToNodeMatrix();
-
- const FaceValue<const double> mes_l = [=] {
- FaceValue<double> compute_mes_l{m_mesh->connectivity()};
- const NodeValue<const TinyVector<Dimension>>& xr = m_mesh->xr();
- if constexpr (Dimension == 1) {
- compute_mes_l.fill(1);
- } else if constexpr (Dimension == 2) {
- parallel_for(
- m_mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
- const auto& face_to_node = face_to_node_matrix[face_id];
- const NodeId node_id1 = face_to_node[0];
- const NodeId node_id2 = face_to_node[1];
- const TinyVector<Dimension, double> r = xr[node_id1] - xr[node_id2];
- compute_mes_l[face_id] = l2Norm(r);
- });
- } else {
- throw NotImplementedError("Not implemented in 3D");
- }
- return compute_mes_l;
- }();
-
- const CellValue<const double> dual_mes_l_j = [=] {
- CellValue<double> compute_mes_j{diamond_mesh->connectivity()};
- mapper->toDualCell(mes_l, compute_mes_j);
-
- return compute_mes_j;
- }();
-
- FaceValue<const double> alpha_l = [&] {
- CellValue<double> alpha_j{diamond_mesh->connectivity()};
-
- parallel_for(
- diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId diamond_cell_id) {
- alpha_j[diamond_cell_id] =
- dual_mes_l_j[diamond_cell_id] * dual_kappaj[diamond_cell_id] / dual_Vj[diamond_cell_id];
- });
-
- FaceValue<double> computed_alpha_l{m_mesh->connectivity()};
- mapper->fromDualCell(alpha_j, computed_alpha_l);
-
- VTKWriter vtk_writer("DD_" + std::to_string(Dimension), 0.01);
-
- vtk_writer.write(diamond_mesh,
- {NamedItemValue{"alphaj", alpha_j}, NamedItemValue{"xj", diamond_mesh_data.xj()},
- NamedItemValue{"Vl", diamond_mesh_data.Vj()}, NamedItemValue{"|l|", dual_mes_l_j}},
- 0,
- true); // forces last output
-
- return computed_alpha_l;
- }();
-
- SparseMatrixDescriptor S{m_mesh->numberOfCells()};
- const auto& face_to_cell_matrix = m_mesh->connectivity().faceToCellMatrix();
- for (FaceId face_id = 0; face_id < m_mesh->numberOfFaces(); ++face_id) {
- const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
- const double beta_l = 0.5 * alpha_l[face_id] * mes_l[face_id];
- const NodeValuePerCell<const TinyVector<Dimension>>& Cjr = mesh_data.Cjr();
-
- // CellValue<double> dual_cell_id{diamond_mesh->connectivity()};
- // mapper->toDualCell(face_id, dual_cell_id);
-
- for (size_t i_cell = 0; i_cell < primal_face_to_cell.size(); ++i_cell) {
- const CellId cell_id1 = primal_face_to_cell[i_cell];
- 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_id1, cell_id2) -= beta_l;
- } else {
- S(cell_id1, cell_id2) += beta_l;
- }
- // S(cell_id1, cell_id2)+= alpha_l[face_id]*(,Cjr[dual_cell_id,dual_node_j]);
- }
- }
- }
- const CellValue<const double> primal_Vj = mesh_data.Vj();
- CRSMatrix A{S};
- Vector<double> b{m_mesh->numberOfCells()};
- double pi = 3.14159265;
- for (CellId i_cell = 0; i_cell < m_mesh->numberOfCells(); ++i_cell) {
- if constexpr (Dimension == 2) {
- b[i_cell] =
- -2 * kappaj[i_cell] * pi * pi * sin(pi * xj[i_cell][0]) * sin(pi * xj[i_cell][1]) * primal_Vj[i_cell];
-
- } else if constexpr (Dimension == 3) {
- b[i_cell] = -3 * kappaj[i_cell] * pi * pi * sin(pi * xj[i_cell][0]) * sin(pi * xj[i_cell][1]) *
- sin(pi * xj[i_cell][2]) * primal_Vj[i_cell];
- }
- }
-
- Vector<double> T{m_mesh->numberOfCells()};
- T = 0;
- PCG{b, A, A, T, 100, 1e-12};
-
- CellValue<double> Temperature{m_mesh->connectivity()};
-
- parallel_for(
- m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { Temperature[cell_id] = T[cell_id]; });
-
- Vector<double> error{m_mesh->numberOfCells()};
- parallel_for(
- m_mesh->numberOfCells(),
- PUGS_LAMBDA(CellId cell_id) { error[cell_id] = Temperature[cell_id] - Tj[cell_id]; });
-
- std::cout << "||Error||_2 = " << std::sqrt((error, error)) << "\n";
-
- {
- VTKWriter vtk_writer("Temperature_" + std::to_string(Dimension), 0.01);
-
- vtk_writer.write(m_mesh, {NamedItemValue{"T", Temperature}}, 0,
- true); // forces last output
- }
- }
- } else {
- throw NotImplementedError("not done in 3d");
- }
- }
-};
+#include <memory>
SchemeModule::SchemeModule()
{