Displace HeatDiamondAlgorithm in its own file

src/language/algorithms/CMakeLists.txt

@@ -3,6 +3,7 @@
 add_library(PugsLanguageAlgorithms
   AcousticSolverAlgorithm.cpp
   Heat5PointsAlgorithm.cpp
+  HeatDiamondAlgorithm.cpp
 )
 
 

src/language/algorithms/HeatDiamondAlgorithm.cpp

+#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&);

src/language/algorithms/HeatDiamondAlgorithm.hpp

+#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

src/language/modules/SchemeModule.cpp

 #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());
-
-        VTKWriter vtk_writer("D_" + std::to_string(Dimension), 0.01);
+#include <scheme/VelocityBoundaryConditionDescriptor.hpp>
 
-        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()
 {