Add order 2 apply_fluxes on LocalDt

src/scheme/Order2LocalDtHyperelasticSolver.cpp

@@ -1136,6 +1136,113 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi
                               Fjr->get<NodeValuePerCell<const Rd>>());
   }
 
+  std::tuple<std::shared_ptr<const MeshVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>>
+  order2_apply_fluxes(const double& dt,
+               const MeshType& mesh,
+               const DiscreteScalarFunction& rho,
+               const DiscreteVectorFunction& u,
+               const DiscreteScalarFunction& E,
+               const DiscreteTensorFunction& CG,
+               const DiscreteScalarFunction& rho_app,
+               const DiscreteTensorFunction& CG_app,
+               const NodeValue<const Rd>& ur,
+               const NodeValuePerCell<const Rd>& Fjr) const
+  {
+    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+    if ((mesh.shared_connectivity() != ur.connectivity_ptr()) or
+        (mesh.shared_connectivity() != Fjr.connectivity_ptr())) {
+      throw NormalError("fluxes are not defined on the same connectivity than the mesh");
+    }
+
+    NodeValue<Rd> new_xr = copy(mesh.xr());
+    parallel_for(
+      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });
+
+    std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr);
+
+    CellValue<const double> Vj           = MeshDataManager::instance().getMeshData(mesh).Vj();
+    const NodeValuePerCell<const Rd> Cjr = MeshDataManager::instance().getMeshData(mesh).Cjr();
+
+    CellValue<double> new_rho = copy(rho.cellValues());
+    CellValue<Rd> new_u       = copy(u.cellValues());
+    CellValue<double> new_E   = copy(E.cellValues());
+    CellValue<Rdxd> new_CG    = copy(CG.cellValues());
+
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+        const auto& cell_nodes = cell_to_node_matrix[j];
+
+        Rd momentum_fluxes   = zero;
+        double energy_fluxes = 0;
+        Rdxd gradv           = zero;
+        for (size_t R = 0; R < cell_nodes.size(); ++R) {
+          const NodeId r = cell_nodes[R];
+          gradv += tensorProduct(ur[r], Cjr(j, R));
+          momentum_fluxes += Fjr(j, R);
+          energy_fluxes += dot(Fjr(j, R), ur[r]);
+        }
+        const Rdxd cauchy_green_fluxes = rho_app[j] * (gradv * CG_app[j] + CG_app[j] * transpose(gradv));
+        const double dt_over_Mj        = dt / (rho[j] * Vj[j]);
+        new_u[j] += dt_over_Mj * momentum_fluxes;
+        new_E[j] += dt_over_Mj * energy_fluxes;
+        new_CG[j] += dt_over_Mj * cauchy_green_fluxes;
+        new_CG[j] += transpose(new_CG[j]);
+        new_CG[j] *= 0.5;
+      });
+
+    CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
+
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { new_rho[j] *= Vj[j] / new_Vj[j]; });
+
+    return {std::make_shared<MeshVariant>(new_mesh),
+            std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_rho)),
+            std::make_shared<DiscreteFunctionVariant>(DiscreteVectorFunction(new_mesh, new_u)),
+            std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_E)),
+            std::make_shared<DiscreteFunctionVariant>(DiscreteTensorFunction(new_mesh, new_CG))};
+  }
+
+  std::tuple<std::shared_ptr<const MeshVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>>
+  order2_apply_fluxes(const double& dt,
+               const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+               const std::shared_ptr<const DiscreteFunctionVariant>& u,
+               const std::shared_ptr<const DiscreteFunctionVariant>& E,
+               const std::shared_ptr<const DiscreteFunctionVariant>& CG,
+               const std::shared_ptr<const DiscreteFunctionVariant>& rho_app,
+               const std::shared_ptr<const DiscreteFunctionVariant>& CG_app,
+               const std::shared_ptr<const ItemValueVariant>& ur,
+               const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr) const
+  {
+    std::shared_ptr mesh_v = getCommonMesh({rho, u, E});
+    if (not mesh_v) {
+      throw NormalError("discrete functions are not defined on the same mesh");
+    }
+
+    if (not checkDiscretizationType({rho, u, E}, DiscreteFunctionType::P0)) {
+      throw NormalError("hyperelastic solver expects P0 functions");
+    }
+
+    return this->order2_apply_fluxes(dt,                                   //
+                              *mesh_v->get<MeshType>(),             //
+                              rho->get<DiscreteScalarFunction>(),   //
+                              u->get<DiscreteVectorFunction>(),     //
+                              E->get<DiscreteScalarFunction>(),     //
+                              CG->get<DiscreteTensorFunction>(),    //
+                              rho_app->get<DiscreteScalarFunction>(),   //
+                              CG_app->get<DiscreteTensorFunction>(),    //
+                              ur->get<NodeValue<const Rd>>(),       //
+                              Fjr->get<NodeValuePerCell<const Rd>>());
+  }
+
   std::tuple<std::shared_ptr<const MeshVariant>,
              std::shared_ptr<const DiscreteFunctionVariant>,
              std::shared_ptr<const DiscreteFunctionVariant>,