Provide a scheme module prototype

It is a first step: it reproduces Pug's main. The only available
parameter is the mesh.

src/language/modules/CMakeLists.txt

@@ -5,6 +5,7 @@ add_library(PugsLanguageModules
   MathModule.cpp
   MeshModule.cpp
   ModuleRepository.cpp
+  SchemeModule.cpp
   VTKModule.cpp
 )
 

src/language/modules/ModuleRepository.cpp

@@ -3,6 +3,7 @@
 #include <language/ast/ASTNode.hpp>
 #include <language/modules/MathModule.hpp>
 #include <language/modules/MeshModule.hpp>
+#include <language/modules/SchemeModule.hpp>
 #include <language/modules/VTKModule.hpp>
 #include <language/utils/BuiltinFunctionEmbedder.hpp>
 #include <language/utils/SymbolTable.hpp>
@@ -20,6 +21,7 @@ ModuleRepository::ModuleRepository()
   this->_subscribe(std::make_unique<MathModule>());
   this->_subscribe(std::make_unique<MeshModule>());
   this->_subscribe(std::make_unique<VTKModule>());
+  this->_subscribe(std::make_unique<SchemeModule>());
 }
 
 template <typename NameEmbedderMapT, typename EmbedderTableT>

src/language/modules/SchemeModule.cpp

+#include <language/modules/SchemeModule.hpp>
+
+#include <language/utils/BuiltinFunctionEmbedder.hpp>
+#include <mesh/Mesh.hpp>
+#include <scheme/AcousticSolver.hpp>
+
+#include <memory>
+
+/////////// TEMPORARY
+
+#include <output/VTKWriter.hpp>
+#include <scheme/BoundaryConditionDescriptor.hpp>
+
+template <size_t Dimension>
+struct MeshDimensionAdapter
+{
+  using ConnectivityType = Connectivity<Dimension>;
+  using MeshType         = Mesh<ConnectivityType>;
+  using MeshDataType     = MeshData<MeshType>;
+  using UnknownsType     = FiniteVolumesEulerUnknowns<MeshDataType>;
+
+  const MeshType& m_mesh;
+
+  MeshDimensionAdapter(const IMesh& mesh) : m_mesh{dynamic_cast<const MeshType&>(mesh)}
+  {
+    std::vector sym_boundary_name_list = [&]() -> std::vector<std::string> {
+      if constexpr (Dimension == 3) {
+        return {"XMIN", "XMAX", "YMIN", "YMAX", "ZMIN", "ZMAX"};
+      } else if constexpr (Dimension == 2) {
+        return {"XMIN", "XMAX", "YMIN", "YMAX"};
+      } else {
+        return {"XMIN", "XMAX"};
+      }
+    }();
+
+    MeshDataType mesh_data(m_mesh);
+    std::vector<std::shared_ptr<BoundaryConditionDescriptor>> bc_descriptor_list;
+    for (const auto& sym_boundary_name : sym_boundary_name_list) {
+      std::shared_ptr<BoundaryDescriptor> boudary_descriptor =
+        std::shared_ptr<BoundaryDescriptor>(new NamedBoundaryDescriptor(sym_boundary_name));
+      SymmetryBoundaryConditionDescriptor* sym_bc_descriptor =
+        new SymmetryBoundaryConditionDescriptor(boudary_descriptor);
+
+      bc_descriptor_list.push_back(std::shared_ptr<BoundaryConditionDescriptor>(sym_bc_descriptor));
+    }
+
+    std::vector<BoundaryConditionHandler> bc_list;
+    {
+      constexpr ItemType FaceType = [] {
+        if constexpr (Dimension > 1) {
+          return ItemType::face;
+        } else {
+          return ItemType::node;
+        }
+      }();
+
+      for (const auto& bc_descriptor : bc_descriptor_list) {
+        switch (bc_descriptor->type()) {
+        case BoundaryConditionDescriptor::Type::symmetry: {
+          const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
+            dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
+          for (size_t i_ref_face_list = 0;
+               i_ref_face_list < m_mesh.connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
+            const auto& ref_face_list = m_mesh.connectivity().template refItemList<FaceType>(i_ref_face_list);
+            const RefId& ref          = ref_face_list.refId();
+            if (ref == sym_bc_descriptor.boundaryDescriptor()) {
+              SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc =
+                new SymmetryBoundaryCondition<MeshType::Dimension>(
+                  MeshFlatNodeBoundary<MeshType::Dimension>(m_mesh, ref_face_list));
+              std::shared_ptr<SymmetryBoundaryCondition<MeshType::Dimension>> bc(sym_bc);
+              bc_list.push_back(BoundaryConditionHandler(bc));
+            }
+          }
+          break;
+        }
+        default: {
+          throw UnexpectedError("Unknown BCDescription\n");
+        }
+        }
+      }
+    }
+
+    UnknownsType unknowns(mesh_data);
+    unknowns.initializeSod();
+
+    AcousticSolver<MeshDataType> acoustic_solver(mesh_data, bc_list);
+
+    const CellValue<const double>& Vj = mesh_data.Vj();
+
+    const double tmax = 0.2;
+    double t          = 0;
+
+    int itermax   = std::numeric_limits<int>::max();
+    int iteration = 0;
+
+    CellValue<double>& rhoj   = unknowns.rhoj();
+    CellValue<double>& ej     = unknowns.ej();
+    CellValue<double>& pj     = unknowns.pj();
+    CellValue<double>& gammaj = unknowns.gammaj();
+    CellValue<double>& cj     = unknowns.cj();
+
+    BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
+
+    VTKWriter vtk_writer("mesh_" + std::to_string(Dimension), 0.01);
+
+    while ((t < tmax) and (iteration < itermax)) {
+      vtk_writer.write(m_mesh,
+                       {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
+                        NamedItemValue{"coords", m_mesh.xr()},
+                        NamedItemValue{"cell_owner", m_mesh.connectivity().cellOwner()},
+                        NamedItemValue{"node_owner", m_mesh.connectivity().nodeOwner()}},
+                       t);
+      double dt = 0.4 * acoustic_solver.acoustic_dt(Vj, cj);
+      if (t + dt > tmax) {
+        dt = tmax - t;
+      }
+      acoustic_solver.computeNextStep(t, dt, unknowns);
+
+      block_eos.updatePandCFromRhoE();
+
+      t += dt;
+      ++iteration;
+    }
+    vtk_writer.write(m_mesh,
+                     {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
+                      NamedItemValue{"coords", m_mesh.xr()},
+                      NamedItemValue{"cell_owner", m_mesh.connectivity().cellOwner()},
+                      NamedItemValue{"node_owner", m_mesh.connectivity().nodeOwner()}},
+                     t, true);   // forces last output
+  }
+};
+
+SchemeModule::SchemeModule()
+{
+  this->_addBuiltinFunction("glace", std::make_shared<BuiltinFunctionEmbedder<void, std::shared_ptr<IMesh>>>(
+                                       std::function<void(std::shared_ptr<IMesh>)>{
+
+                                         [](std::shared_ptr<const IMesh> p_mesh) -> void {
+                                           switch (p_mesh->dimension()) {
+                                           case 1: {
+                                             MeshDimensionAdapter<1>{*p_mesh};
+                                             break;
+                                           }
+                                           case 2: {
+                                             MeshDimensionAdapter<2>{*p_mesh};
+                                             break;
+                                           }
+                                           case 3: {
+                                             MeshDimensionAdapter<3>{*p_mesh};
+                                             break;
+                                           }
+                                           default: {
+                                             throw UnexpectedError("invalid mesh dimension");
+                                           }
+                                           }
+                                         }}
+
+                                       ));
+}

src/language/modules/SchemeModule.hpp

+#ifndef SCHEME_MODULE_HPP
+#define SCHEME_MODULE_HPP
+
+#include <language/modules/BuiltinModule.hpp>
+#include <utils/PugsMacros.hpp>
+
+class SchemeModule : public BuiltinModule
+{
+ public:
+  std::string_view
+  name() const final
+  {
+    return "scheme";
+  }
+
+  SchemeModule();
+
+  ~SchemeModule() = default;
+};
+
+#endif   // SCHEME_MODULE_HPP