Extend time step calculation for polygonal and cubic meshes

src/geometry/LineTransformation.hpp

@@ -73,12 +73,20 @@ class LineTransformation
     return m_velocity;
   }
 
+  PUGS_INLINE
   double
   velocityNorm() const
   {
     return m_velocity_norm;
   }
 
+  PUGS_INLINE
+  double
+  velocityNorm(const TinyVector<1>&) const
+  {
+    return m_velocity_norm;
+  }
+
   PUGS_INLINE
   LineTransformation(const TinyVector<Dimension>& a, const TinyVector<Dimension>& b)
     : m_velocity{0.5 * (b - a)}, m_velocity_norm{l2Norm(m_velocity)}, m_shift{0.5 * (a + b)}

src/scheme/VariationalSolver.cpp

@@ -24,6 +24,7 @@
 
 #include <analysis/GaussQuadratureDescriptor.hpp>
 #include <analysis/QuadratureManager.hpp>
+#include <geometry/LineCubicTransformation.hpp>
 #include <geometry/LineParabolicTransformation.hpp>
 
 #include <algebra/CRSMatrixDescriptor.hpp>
@@ -38,40 +39,80 @@
 #include <analysis/GaussLegendreQuadratureDescriptor.hpp>
 #include <language/utils/IntegrateOnCells.hpp>
 
-double
-variational_acoustic_dt(const std::shared_ptr<const DiscreteFunctionVariant>& c_v)
+template <MeshConcept MeshType, size_t mesh_degree>
+struct LineTransformationAccessor
 {
-  const auto& c = c_v->get<DiscreteFunctionP0<const double>>();
+};
 
-  return std::visit(
-    [&](auto&& p_mesh) -> double {
-      const auto& mesh = *p_mesh;
+template <>
+struct LineTransformationAccessor<Mesh<2>, 1>
+{
+  const NodeValue<const TinyVector<2>> m_xr;
 
-      using MeshType = mesh_type_t<decltype(mesh)>;
-      if constexpr (is_polynomial_mesh_v<MeshType>) {
-        const auto Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
-        if (mesh.degree() != 2) {
-          throw NotImplementedError("mesh degree != 2");
+  const ItemToItemMatrix<ItemType::face, ItemType::node> m_face_to_node_matrix;
+
+  LineTransformationAccessor(const Mesh<2>& mesh)
+    : m_xr{mesh.xr()}, m_face_to_node_matrix(mesh.connectivity().faceToNodeMatrix())
+  {}
+
+  auto
+  getTransformation(const FaceId face_id) const
+  {
+    auto node_list = m_face_to_node_matrix[face_id];
+    return LineTransformation<2>{m_xr[node_list[0]], m_xr[node_list[1]]};
+  }
+};
+
+template <size_t mesh_degree>
+struct LineTransformationAccessor<PolynomialMesh<2>, mesh_degree>
+{
+  const NodeValue<const TinyVector<2>> m_xr;
+  const FaceArray<const TinyVector<2>> m_xl;
+
+  const ItemToItemMatrix<ItemType::face, ItemType::node> m_face_to_node_matrix;
+
+  LineTransformationAccessor(const PolynomialMesh<2>& mesh)
+    : m_xr{mesh.xr()}, m_xl{mesh.xl()}, m_face_to_node_matrix(mesh.connectivity().faceToNodeMatrix())
+  {}
+
+  auto
+  getTransformation(const FaceId face_id) const
+  {
+    auto node_list = m_face_to_node_matrix[face_id];
+    if constexpr (mesh_degree == 1) {
+      return LineTransformation<2>{m_xr[node_list[0]], m_xr[node_list[1]]};
+    } else if constexpr (mesh_degree == 2) {
+      return LineParabolicTransformation<2>{m_xr[node_list[0]], m_xl[face_id][0], m_xr[node_list[1]]};
+    } else if constexpr (mesh_degree == 3) {
+      return LineCubicTransformation<2>{m_xr[node_list[0]], m_xl[face_id][0], m_xl[face_id][1], m_xr[node_list[1]]};
+    } else {
+      static_assert(mesh_degree == 1, "mesh degree is not supported");
     }
+  }
+};
 
-        auto xr = mesh.xr();
-        auto xl = mesh.xl();
+template <size_t mesh_degree, MeshConcept MeshType>
+double
+variational_acoustic_dt(const DiscreteFunctionP0<const double>& c, const std::shared_ptr<const MeshType>& p_mesh)
+{
+  const auto& mesh = *p_mesh;
+  const auto Vj    = MeshDataManager::instance().getMeshData(mesh).Vj();
 
   auto cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
-        auto face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
 
   auto qf = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor(2));
 
   CellValue<double> Sj{mesh.connectivity()};
   parallel_for(
     mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+      LineTransformationAccessor<MeshType, mesh_degree> transform(*p_mesh);
+
       double sum     = 0;
       auto face_list = cell_to_face_matrix[cell_id];
       for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
         const FaceId face_id = face_list[i_face];
-              auto node_list       = face_to_node_matrix[face_id];
 
-              LineParabolicTransformation<MeshType::Dimension> t(xr[node_list[0]], xl[face_id][0], xr[node_list[1]]);
+        auto t = transform.getTransformation(face_id);
 
         for (size_t iq = 0; iq < qf.numberOfPoints(); ++iq) {
           sum += qf.weight(iq) * t.velocityNorm(qf.point(iq));
@@ -85,6 +126,39 @@ variational_acoustic_dt(const std::shared_ptr<const DiscreteFunctionVariant>& c_
     mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { local_dt[j] = 2 * Vj[j] / (Sj[j] * c[j]); });
 
   return min(local_dt);
+}
+
+double
+variational_acoustic_dt(const std::shared_ptr<const DiscreteFunctionVariant>& c_v)
+{
+  const auto& c = c_v->get<DiscreteFunctionP0<const double>>();
+
+  return std::visit(
+    [&](auto&& p_mesh) -> double {
+      const auto& mesh = *p_mesh;
+
+      using MeshType = mesh_type_t<decltype(mesh)>;
+      if constexpr (is_polygonal_mesh_v<MeshType>) {
+        if constexpr (MeshType::Dimension == 2) {
+          return variational_acoustic_dt<1>(c, p_mesh);
+        } else {
+          throw NotImplementedError("not implemented in dimension d != 2");
+        }
+      } else if constexpr (is_polynomial_mesh_v<MeshType>) {
+        switch (mesh.degree()) {
+        case 1: {
+          return variational_acoustic_dt<1>(c, p_mesh);
+        }
+        case 2: {
+          return variational_acoustic_dt<2>(c, p_mesh);
+        }
+        case 3: {
+          return variational_acoustic_dt<3>(c, p_mesh);
+        }
+        default: {
+          throw NotImplementedError("not implemented for mesh degree > 3");
+        }
+        }
       } else {
         throw NormalError("unexpected mesh type");
       }
@@ -206,23 +280,13 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
 
           Table<TinyVector<Dimension>> quadrature_points(mesh_face_boundary.faceList().size(), qf.numberOfPoints());
 
+          LineTransformationAccessor<MeshType, 2> transformation{mesh};
+
           for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
             const FaceId face_id = face_list[i_face];
 
-            const Rd& x0 = mesh.xr()[face_to_node_matrix[face_id][0]];
-            const Rd& x2 = mesh.xr()[face_to_node_matrix[face_id][1]];
-
-            const Rd x1 = [&]() {
-              if (mesh.degree() == 2) {
-                return mesh.xl()[face_id][0];
-              } else if (mesh.degree() == 1) {
-                return 0.5 * (x0 + x2);
-              } else {
-                throw NotImplementedError("degree>2");
-              }
-            }();
+            auto t = transformation.getTransformation(face_id);
 
-            LineParabolicTransformation<Dimension> t(x0, x1, x2);
             auto face_quadrature_points = quadrature_points[i_face];
             for (size_t i_xi = 0; i_xi < qf.numberOfPoints(); ++i_xi) {
               face_quadrature_points[i_xi] = t(qf.point(i_xi));
@@ -731,7 +795,7 @@ class VariationalSolverHandler::VariationalSolver final : public VariationalSolv
     const CellValue<const double> rhoc =
       (rho_v->get<DiscreteFunctionP0<const double>>() * c_v->get<DiscreteFunctionP0<const double>>()).cellValues();
 
-#if 1
+#if 0
     // Heun's RK3 method
 
     auto [ur1, ul1, momentum_fluxes1, energy_fluxes1] =