diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 57aea488cd78fc4bd8b2d3f13a2b501eca8cc9ce..1f51634718f8f5bda9b510ddb83ba78bf0e8ec71 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -496,7 +496,7 @@ SchemeModule::SchemeModule()
));
- this->_addBuiltinFunction("local_dt_eucclhyd_solver_order2",
+ this->_addBuiltinFunction("local_dt_eucclhyd_solver_order2",
std::function(
[](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
@@ -514,17 +514,17 @@ SchemeModule::SchemeModule()
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list2)
-> std::tuple<
- std::shared_ptr<const IMesh>, std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const IMesh>,
+ std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const MeshVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>> {
return LocalDtAcousticSolverHandler{getCommonMesh({rho1, u1, E1, c1, p1}),
getCommonMesh({rho2, u2, E2, c2, p2})}
.solver()
- .applyOrder2(LocalDtAcousticSolverHandler::SolverType::Eucclhyd, dt1, rho1, rho2, u1, u2,
- E1, E2, c1, c2, p1, p2, bc_descriptor_list1, bc_descriptor_list2);
+ .applyOrder2(LocalDtAcousticSolverHandler::SolverType::Eucclhyd, dt1, rho1, rho2, u1,
+ u2, E1, E2, c1, c2, p1, p2, bc_descriptor_list1, bc_descriptor_list2);
}
));
@@ -548,9 +548,9 @@ SchemeModule::SchemeModule()
bc_descriptor_list2,
const size_t& q)
-> std::tuple<
- std::shared_ptr<const IMesh>, std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const IMesh>,
+ std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const MeshVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>> {
@@ -581,17 +581,17 @@ SchemeModule::SchemeModule()
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list2)
-> std::tuple<
- std::shared_ptr<const IMesh>, std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const IMesh>,
+ std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const MeshVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>> {
return LocalDtAcousticSolverHandler{getCommonMesh({rho1, u1, E1, c1, p1}),
getCommonMesh({rho2, u2, E2, c2, p2})}
.solver()
- .apply(LocalDtAcousticSolverHandler::SolverType::Glace, dt1,rho1, rho2, u1, u2,
- E1, E2, c1, c2, p1, p2, bc_descriptor_list1, bc_descriptor_list2);
+ .apply(LocalDtAcousticSolverHandler::SolverType::Glace, dt1, rho1, rho2, u1, u2, E1,
+ E2, c1, c2, p1, p2, bc_descriptor_list1, bc_descriptor_list2);
}
));
@@ -615,9 +615,9 @@ SchemeModule::SchemeModule()
bc_descriptor_list2,
const size_t& q)
-> std::tuple<
- std::shared_ptr<const IMesh>, std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const IMesh>,
+ std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const MeshVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>> {
@@ -630,7 +630,7 @@ SchemeModule::SchemeModule()
));
- this->_addBuiltinFunction("local_dt_eucclhyd_solver",
+ this->_addBuiltinFunction("local_dt_eucclhyd_solver",
std::function(
[](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
@@ -648,16 +648,16 @@ SchemeModule::SchemeModule()
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list2)
-> std::tuple<
- std::shared_ptr<const IMesh>, std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const IMesh>,
+ std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const MeshVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>> {
return LocalDtAcousticSolverHandler{getCommonMesh({rho1, u1, E1, c1, p1}),
getCommonMesh({rho2, u2, E2, c2, p2})}
.solver()
- .apply(LocalDtAcousticSolverHandler::SolverType::Eucclhyd, dt1,rho1, rho2, u1, u2,
+ .apply(LocalDtAcousticSolverHandler::SolverType::Eucclhyd, dt1, rho1, rho2, u1, u2,
E1, E2, c1, c2, p1, p2, bc_descriptor_list1, bc_descriptor_list2);
}
@@ -784,7 +784,7 @@ SchemeModule::SchemeModule()
const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list1,
- const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
const std::shared_ptr<const DiscreteFunctionVariant>& u2,
const std::shared_ptr<const DiscreteFunctionVariant>& E2,
const std::shared_ptr<const DiscreteFunctionVariant>& CG2,
@@ -793,24 +793,24 @@ SchemeModule::SchemeModule()
const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list2,
- const double& mu,
- const double& lambda,
- const double& dt1,
- const size_t& q) -> std::tuple<std::shared_ptr<const IMesh>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const IMesh>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>> {
- return LocalDtHyperelasticSolverHandler{getCommonMesh({rho1, u1, E1, CG1, aL1, aT1, sigma1}),
- getCommonMesh({rho2, u2, E2, CG2, aL2, aT2, sigma2})}
+ const double& mu, const double& lambda, const double& dt1, const size_t& q)
+ -> 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>, 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>> {
+ return LocalDtHyperelasticSolverHandler{getCommonMesh(
+ {rho1, u1, E1, CG1, aL1, aT1, sigma1}),
+ getCommonMesh(
+ {rho2, u2, E2, CG2, aL2, aT2, sigma2})}
.solver()
- .apply(LocalDtHyperelasticSolverHandler::SolverType::Eucclhyd, dt1, q, rho1, rho2, u1, u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2,
- sigma1, sigma2, bc_descriptor_list1, bc_descriptor_list2, mu, lambda);
+ .apply(LocalDtHyperelasticSolverHandler::SolverType::Eucclhyd, dt1, q, rho1, rho2, u1,
+ u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, bc_descriptor_list1,
+ bc_descriptor_list2, mu, lambda);
}
));
@@ -827,7 +827,7 @@ SchemeModule::SchemeModule()
const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list1,
- const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
const std::shared_ptr<const DiscreteFunctionVariant>& u2,
const std::shared_ptr<const DiscreteFunctionVariant>& E2,
const std::shared_ptr<const DiscreteFunctionVariant>& CG2,
@@ -836,24 +836,24 @@ SchemeModule::SchemeModule()
const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list2,
- const double& mu,
- const double& lambda,
- const double& dt1,
- const size_t& q) -> std::tuple<std::shared_ptr<const IMesh>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const IMesh>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>> {
- return LocalDtHyperelasticSolverHandler{getCommonMesh({rho1, u1, E1, CG1, aL1, aT1, sigma1}),
- getCommonMesh({rho2, u2, E2, CG2, aL2, aT2, sigma2})}
+ const double& mu, const double& lambda, const double& dt1, const size_t& q)
+ -> 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>, 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>> {
+ return LocalDtHyperelasticSolverHandler{getCommonMesh(
+ {rho1, u1, E1, CG1, aL1, aT1, sigma1}),
+ getCommonMesh(
+ {rho2, u2, E2, CG2, aL2, aT2, sigma2})}
.solver()
- .apply(LocalDtHyperelasticSolverHandler::SolverType::Glace, dt1, q, rho1, rho2, u1, u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2,
- sigma1, sigma2, bc_descriptor_list1, bc_descriptor_list2, mu, lambda);
+ .apply(LocalDtHyperelasticSolverHandler::SolverType::Glace, dt1, q, rho1, rho2, u1,
+ u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, bc_descriptor_list1,
+ bc_descriptor_list2, mu, lambda);
}
));
diff --git a/src/scheme/LocalDtAcousticSolver.cpp b/src/scheme/LocalDtAcousticSolver.cpp
index 62091415cda41c2aac061bb041320a236d35197b..b653baa4830736e4ed11a8ed8e4cf48a51f6089d 100644
--- a/src/scheme/LocalDtAcousticSolver.cpp
+++ b/src/scheme/LocalDtAcousticSolver.cpp
@@ -6,7 +6,11 @@
#include <mesh/MeshFaceBoundary.hpp>
#include <mesh/MeshFlatNodeBoundary.hpp>
#include <mesh/MeshNodeBoundary.hpp>
+#include <mesh/MeshTraits.hpp>
+#include <mesh/MeshVariant.hpp>
#include <mesh/SubItemValuePerItemVariant.hpp>
+#include <scheme/AcousticSolver.hpp>
+#include <scheme/CouplingBoundaryConditionDescriptor.hpp>
#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
#include <scheme/DiscreteFunctionP0.hpp>
#include <scheme/DiscreteFunctionUtils.hpp>
@@ -15,26 +19,25 @@
#include <scheme/IBoundaryConditionDescriptor.hpp>
#include <scheme/IDiscreteFunctionDescriptor.hpp>
#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
-#include <scheme/CouplingBoundaryConditionDescriptor.hpp>
-#include <scheme/AcousticSolver.hpp>
#include <utils/Socket.hpp>
#include <variant>
#include <vector>
-template <size_t Dimension>
+template <MeshConcept MeshType>
class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
: public LocalDtAcousticSolverHandler::ILocalDtAcousticSolver
{
private:
+ static constexpr size_t Dimension = MeshType::Dimension;
+
using Rdxd = TinyMatrix<Dimension>;
using Rd = TinyVector<Dimension>;
- using MeshType = Mesh<Connectivity<Dimension>>;
- using MeshDataType = MeshData<Dimension>;
+ using MeshDataType = MeshData<MeshType>;
- using DiscreteScalarFunction = DiscreteFunctionP0<Dimension, const double>;
- using DiscreteVectorFunction = DiscreteFunctionP0<Dimension, const Rd>;
+ using DiscreteScalarFunction = DiscreteFunctionP0<const double>;
+ using DiscreteVectorFunction = DiscreteFunctionP0<const Rd>;
class FixedBoundaryCondition;
class PressureBoundaryCondition;
@@ -48,7 +51,7 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
SymmetryBoundaryCondition,
VelocityBoundaryCondition,
ExternalFSIVelocityBoundaryCondition,
- CouplingBoundaryCondition>;
+ CouplingBoundaryCondition>;
using BoundaryConditionList = std::vector<BoundaryCondition>;
@@ -168,7 +171,7 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
}
NodeValue<Rd>
- _computeBr(const Mesh<Connectivity<Dimension>>& mesh,
+ _computeBr(const Mesh<Dimension>& mesh,
const NodeValuePerCell<const Rdxd>& Ajr,
const DiscreteVectorFunction& u,
const DiscreteScalarFunction& p) const
@@ -235,8 +238,7 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
if (dirichlet_bc_descriptor.name() == "velocity") {
- MeshNodeBoundary<Dimension> mesh_node_boundary =
- getMeshNodeBoundary(mesh, dirichlet_bc_descriptor.boundaryDescriptor());
+ MeshNodeBoundary mesh_node_boundary = getMeshNodeBoundary(mesh, dirichlet_bc_descriptor.boundaryDescriptor());
Array<const Rd> value_list =
InterpolateItemValue<Rd(Rd)>::template interpolate<ItemType::node>(dirichlet_bc_descriptor.rhsSymbolId(),
@@ -248,17 +250,15 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
const FunctionSymbolId pressure_id = dirichlet_bc_descriptor.rhsSymbolId();
if constexpr (Dimension == 1) {
- MeshNodeBoundary<Dimension> mesh_node_boundary =
- getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());
+ MeshNodeBoundary mesh_node_boundary = getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());
Array<const double> node_values =
InterpolateItemValue<double(Rd)>::template interpolate<ItemType::node>(pressure_id, mesh.xr(),
mesh_node_boundary.nodeList());
-
- bc_list.emplace_back(PressureBoundaryCondition{mesh_node_boundary, node_values});
+ PressureBoundaryCondition p(mesh_node_boundary, node_values);
+ bc_list.emplace_back(p);
} else {
- MeshFaceBoundary<Dimension> mesh_face_boundary =
- getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor());
+ MeshFaceBoundary mesh_face_boundary = getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor());
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
Array<const double> face_values =
@@ -273,10 +273,10 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
break;
}
case IBoundaryConditionDescriptor::Type::coupling: {
- const CouplingBoundaryConditionDescriptor& coupling_bc_descriptor =
- dynamic_cast<const CouplingBoundaryConditionDescriptor&>(*bc_descriptor);
- bc_list.emplace_back(
- CouplingBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),coupling_bc_descriptor.label()));
+ const CouplingBoundaryConditionDescriptor& coupling_bc_descriptor =
+ dynamic_cast<const CouplingBoundaryConditionDescriptor&>(*bc_descriptor);
+ bc_list.emplace_back(CouplingBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+ coupling_bc_descriptor.label()));
break;
}
default: {
@@ -301,23 +301,23 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
NodeValue<Rdxd>& Ar,
NodeValue<Rd>& br) const;
void _applyCouplingBC(NodeValue<Rdxd>& Ar1,
- NodeValue<Rdxd>& Ar2,
- NodeValue<Rd>& br1,
- NodeValue<Rd>& br2,
- const std::vector<NodeId>& map1,
- const std::vector<NodeId>& map2) const;
+ NodeValue<Rdxd>& Ar2,
+ NodeValue<Rd>& br1,
+ NodeValue<Rd>& br2,
+ const std::vector<NodeId>& map1,
+ const std::vector<NodeId>& map2) const;
void _applyCouplingBC(const MeshType& mesh,
- NodeValue<Rd>& ur,
- NodeValue<Rd>& CR_ur,
- NodeValuePerCell<Rd>& Fjr,
- NodeValuePerCell<Rd>& CR_Fjr,
- const std::vector<NodeId>& map2) const;
+ NodeValue<Rd>& ur,
+ NodeValue<Rd>& CR_ur,
+ NodeValuePerCell<Rd>& Fjr,
+ NodeValuePerCell<Rd>& CR_Fjr,
+ const std::vector<NodeId>& map2) const;
void _applyCouplingBC2(NodeValue<Rdxd>& Ar1,
- NodeValue<Rdxd>& Ar2,
- NodeValue<Rd>& ur1,
- NodeValue<Rd>& ur2,
- const std::vector<NodeId>& map1,
- const std::vector<NodeId>& map2) const;
+ NodeValue<Rdxd>& Ar2,
+ NodeValue<Rd>& ur1,
+ NodeValue<Rd>& ur2,
+ const std::vector<NodeId>& map1,
+ const std::vector<NodeId>& map2) const;
void
_applyBoundaryConditions(const BoundaryConditionList& bc_list,
@@ -366,58 +366,59 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
return F;
}
- std::tuple<std::vector<NodeId>,
- std::vector<NodeId>>
- _computeMapping(std::shared_ptr<const IMesh>& i_mesh1,
- std::shared_ptr<const IMesh>& i_mesh2,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2) const
+ std::tuple<std::vector<NodeId>, std::vector<NodeId>>
+ _computeMapping(std::shared_ptr<const MeshVariant>& i_mesh1,
+ std::shared_ptr<const MeshVariant>& i_mesh2,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2) const
{
- const MeshType& mesh1 = dynamic_cast<const MeshType&>(*i_mesh1);
- const MeshType& mesh2 = dynamic_cast<const MeshType&>(*i_mesh2);
-
+ const MeshType& mesh1 = *i_mesh1->get<MeshType>();
+ const MeshType& mesh2 = *i_mesh2->get<MeshType>();
+
const BoundaryConditionList bc_list1 = this->_getBCList(mesh1, bc_descriptor_list1);
const BoundaryConditionList bc_list2 = this->_getBCList(mesh2, bc_descriptor_list2);
-
+
std::vector<NodeId> map1;
std::vector<NodeId> map2;
-
+
NodeValue<Rd> xr1 = copy(mesh1.xr());
NodeValue<Rd> xr2 = copy(mesh2.xr());
- for(const auto& boundary_condition1 : bc_list1){
- std::visit([&](auto&& bc1) {
- using T1 = std::decay_t<decltype(bc1)>;
- if constexpr (std::is_same_v<CouplingBoundaryCondition,T1>){
- const auto& node_list1 = bc1.nodeList();
- for(const auto& boundary_condition2 : bc_list2){
- std::visit([&](auto&& bc2) {
- using T2 = std::decay_t<decltype(bc2)>;
- if constexpr (std::is_same_v<CouplingBoundaryCondition,T2>){
- if(bc1.label() == bc2.label()){
- const auto& node_list2 = bc2.nodeList();
- for(size_t i = 0; i<node_list1.size(); i++){
- NodeId node_id1 = node_list1[i];
- NodeId node_id2;
- for(size_t j = 0; j<node_list2.size(); j++){
- node_id2 = node_list2[j];
- double err = 0;
- for(size_t j = 0; j<Dimension; j++){
- err += (xr1[node_id1][j] - xr2[node_id2][j])*(xr1[node_id1][j] - xr2[node_id2][j]);
- }
- if(sqrt(err) < 1e-14){
- map1.push_back(node_id1);
- map2.push_back(node_id2);
- }
- }
- }
- }
- }
- },boundary_condition2);
- }
- }
- },boundary_condition1);
+ for (const auto& boundary_condition1 : bc_list1) {
+ std::visit(
+ [&](auto&& bc1) {
+ using T1 = std::decay_t<decltype(bc1)>;
+ if constexpr (std::is_same_v<CouplingBoundaryCondition, T1>) {
+ const auto& node_list1 = bc1.nodeList();
+ for (const auto& boundary_condition2 : bc_list2) {
+ std::visit(
+ [&](auto&& bc2) {
+ using T2 = std::decay_t<decltype(bc2)>;
+ if constexpr (std::is_same_v<CouplingBoundaryCondition, T2>) {
+ if (bc1.label() == bc2.label()) {
+ const auto& node_list2 = bc2.nodeList();
+ for (size_t i = 0; i < node_list1.size(); i++) {
+ NodeId node_id1 = node_list1[i];
+ NodeId node_id2;
+ for (size_t j = 0; j < node_list2.size(); j++) {
+ node_id2 = node_list2[j];
+ double err = 0;
+ err = dot((xr1[node_id1] - xr2[node_id2]), (xr1[node_id1] - xr2[node_id2]));
+ if (sqrt(err) < 1e-14) {
+ map1.push_back(node_id1);
+ map2.push_back(node_id2);
+ }
+ }
+ }
+ }
+ }
+ },
+ boundary_condition2);
+ }
+ }
+ },
+ boundary_condition1);
}
- return std::make_tuple(map1,map2);
+ return std::make_tuple(map1, map2);
}
public:
@@ -438,7 +439,7 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
throw NormalError("acoustic solver expects P0 functions");
}
- const MeshType& mesh = dynamic_cast<const MeshType&>(*i_mesh);
+ const MeshType& mesh = *i_mesh->get<MeshType>();
const DiscreteScalarFunction& rho = rho_v->get<DiscreteScalarFunction>();
const DiscreteScalarFunction& c = c_v->get<DiscreteScalarFunction>();
const DiscreteVectorFunction& u = u_v->get<DiscreteVectorFunction>();
@@ -464,11 +465,11 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
}
std::tuple<const std::shared_ptr<const ItemValueVariant>,
- const std::shared_ptr<const SubItemValuePerItemVariant>,
- const std::shared_ptr<const ItemValueVariant>,
- const std::shared_ptr<const SubItemValuePerItemVariant>,
- NodeValue<Rd>,
- NodeValuePerCell<Rd>>
+ const std::shared_ptr<const SubItemValuePerItemVariant>,
+ const std::shared_ptr<const ItemValueVariant>,
+ const std::shared_ptr<const SubItemValuePerItemVariant>,
+ NodeValue<Rd>,
+ NodeValuePerCell<Rd>>
compute_fluxes(const SolverType& solver_type,
const std::shared_ptr<const DiscreteFunctionVariant>& rho1_v,
const std::shared_ptr<const DiscreteFunctionVariant>& c1_v,
@@ -480,8 +481,8 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
const std::shared_ptr<const DiscreteFunctionVariant>& u2_v,
const std::shared_ptr<const DiscreteFunctionVariant>& p2_v,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2,
- const std::vector<NodeId>& map1,
- const std::vector<NodeId>& map2) const
+ const std::vector<NodeId>& map1,
+ const std::vector<NodeId>& map2) const
{
std::shared_ptr i_mesh1 = getCommonMesh({rho1_v, c1_v, u1_v, p1_v});
std::shared_ptr i_mesh2 = getCommonMesh({rho2_v, c2_v, u2_v, p2_v});
@@ -501,14 +502,13 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
throw NormalError("acoustic solver expects P0 functions");
}
-
- const MeshType& mesh1 = dynamic_cast<const MeshType&>(*i_mesh1);
+ const MeshType& mesh1 = *i_mesh1->get<MeshType>();
const DiscreteScalarFunction& rho1 = rho1_v->get<DiscreteScalarFunction>();
const DiscreteScalarFunction& c1 = c1_v->get<DiscreteScalarFunction>();
const DiscreteVectorFunction& u1 = u1_v->get<DiscreteVectorFunction>();
const DiscreteScalarFunction& p1 = p1_v->get<DiscreteScalarFunction>();
- const MeshType& mesh2 = dynamic_cast<const MeshType&>(*i_mesh2);
+ const MeshType& mesh2 = *i_mesh2->get<MeshType>();
const DiscreteScalarFunction& rho2 = rho2_v->get<DiscreteScalarFunction>();
const DiscreteScalarFunction& c2 = c2_v->get<DiscreteScalarFunction>();
const DiscreteVectorFunction& u2 = u2_v->get<DiscreteVectorFunction>();
@@ -528,26 +528,24 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
this->_applyBoundaryConditions(bc_list2, mesh2, Ar2, br2);
this->_applyExternalVelocityBC(bc_list1, p1, Ar1, br1);
this->_applyExternalVelocityBC(bc_list2, p2, Ar2, br2);
- this->_applyCouplingBC(Ar1,Ar2,br1,br2,map1,map2);
+ this->_applyCouplingBC(Ar1, Ar2, br1, br2, map1, map2);
synchronize(Ar1);
synchronize(br1);
synchronize(Ar2);
synchronize(br2);
- NodeValue<Rd> ur1 = this->_computeUr(mesh1, Ar1, br1);
- NodeValuePerCell<Rd> Fjr1 = this->_computeFjr(mesh1, Ajr1, ur1, u1, p1);
- NodeValue<Rd> ur2 = this->_computeUr(mesh2, Ar2, br2);
- NodeValuePerCell<Rd> Fjr2 = this->_computeFjr(mesh2, Ajr2, ur2, u2, p2);
- NodeValue<Rd> CR_ur = this->_computeUr(mesh2, Ar2, br2);
+ NodeValue<Rd> ur1 = this->_computeUr(mesh1, Ar1, br1);
+ NodeValuePerCell<Rd> Fjr1 = this->_computeFjr(mesh1, Ajr1, ur1, u1, p1);
+ NodeValue<Rd> ur2 = this->_computeUr(mesh2, Ar2, br2);
+ NodeValuePerCell<Rd> Fjr2 = this->_computeFjr(mesh2, Ajr2, ur2, u2, p2);
+ NodeValue<Rd> CR_ur = this->_computeUr(mesh2, Ar2, br2);
NodeValuePerCell<Rd> CR_Fjr = this->_computeFjr(mesh2, Ajr2, ur2, u2, p2);
return std::make_tuple(std::make_shared<const ItemValueVariant>(ur1),
std::make_shared<const SubItemValuePerItemVariant>(Fjr1),
- std::make_shared<const ItemValueVariant>(ur2),
- std::make_shared<const SubItemValuePerItemVariant>(Fjr2),
- CR_ur,
- CR_Fjr);
+ std::make_shared<const ItemValueVariant>(ur2),
+ std::make_shared<const SubItemValuePerItemVariant>(Fjr2), CR_ur, CR_Fjr);
}
std::tuple<const std::shared_ptr<const ItemValueVariant>, const std::shared_ptr<const SubItemValuePerItemVariant>>
@@ -557,9 +555,9 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
const std::shared_ptr<const DiscreteFunctionVariant>& p_v,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
- NodeValue<Rd> CR_ur,
- NodeValuePerCell<Rd> CR_Fjr,
- const std::vector<NodeId> map2) const
+ NodeValue<Rd> CR_ur,
+ NodeValuePerCell<Rd> CR_Fjr,
+ const std::vector<NodeId> map2) const
{
std::shared_ptr i_mesh = getCommonMesh({rho_v, c_v, u_v, p_v});
if (not i_mesh) {
@@ -570,7 +568,7 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
throw NormalError("acoustic solver expects P0 functions");
}
- const MeshType& mesh = dynamic_cast<const MeshType&>(*i_mesh);
+ const MeshType& mesh = *i_mesh->get<MeshType>();
const DiscreteScalarFunction& rho = rho_v->get<DiscreteScalarFunction>();
const DiscreteScalarFunction& c = c_v->get<DiscreteScalarFunction>();
const DiscreteVectorFunction& u = u_v->get<DiscreteVectorFunction>();
@@ -580,7 +578,7 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
NodeValue<Rd> br = this->_computeBr(mesh, Ajr, u, p);
-
+
const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
this->_applyBoundaryConditions(bc_list, mesh, Ar, br);
this->_applyExternalVelocityBC(bc_list, p, Ar, br);
@@ -591,31 +589,31 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
NodeValue<Rd> ur = this->_computeUr(mesh, Ar, br);
NodeValuePerCell<Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u, p);
- this->_applyCouplingBC(mesh,ur,CR_ur,Fjr,CR_Fjr,map2);
+ this->_applyCouplingBC(mesh, ur, CR_ur, Fjr, CR_Fjr, map2);
return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
std::make_shared<const SubItemValuePerItemVariant>(Fjr));
}
- std::tuple<const std::shared_ptr<const ItemValueVariant>,
- const std::shared_ptr<const SubItemValuePerItemVariant>,
- const std::shared_ptr<const ItemValueVariant>,
- const std::shared_ptr<const SubItemValuePerItemVariant>,
- NodeValue<Rd>,
- NodeValuePerCell<Rd>>
+ std::tuple<const std::shared_ptr<const ItemValueVariant>,
+ const std::shared_ptr<const SubItemValuePerItemVariant>,
+ const std::shared_ptr<const ItemValueVariant>,
+ const std::shared_ptr<const SubItemValuePerItemVariant>,
+ NodeValue<Rd>,
+ NodeValuePerCell<Rd>>
compute_fluxes2(const SolverType& solver_type,
- const std::shared_ptr<const DiscreteFunctionVariant>& rho1_v,
- const std::shared_ptr<const DiscreteFunctionVariant>& c1_v,
- const std::shared_ptr<const DiscreteFunctionVariant>& u1_v,
- const std::shared_ptr<const DiscreteFunctionVariant>& p1_v,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
- const std::shared_ptr<const DiscreteFunctionVariant>& rho2_v,
- const std::shared_ptr<const DiscreteFunctionVariant>& c2_v,
- const std::shared_ptr<const DiscreteFunctionVariant>& u2_v,
- const std::shared_ptr<const DiscreteFunctionVariant>& p2_v,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2,
- const std::vector<NodeId>& map1,
- const std::vector<NodeId>& map2) const
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho1_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c1_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u1_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p1_v,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho2_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c2_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u2_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p2_v,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2,
+ const std::vector<NodeId>& map1,
+ const std::vector<NodeId>& map2) const
{
std::shared_ptr i_mesh1 = getCommonMesh({rho1_v, c1_v, u1_v, p1_v});
std::shared_ptr i_mesh2 = getCommonMesh({rho2_v, c2_v, u2_v, p2_v});
@@ -635,14 +633,13 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
throw NormalError("acoustic solver expects P0 functions");
}
-
- const MeshType& mesh1 = dynamic_cast<const MeshType&>(*i_mesh1);
+ const MeshType& mesh1 = *i_mesh1->get<MeshType>();
const DiscreteScalarFunction& rho1 = rho1_v->get<DiscreteScalarFunction>();
const DiscreteScalarFunction& c1 = c1_v->get<DiscreteScalarFunction>();
const DiscreteVectorFunction& u1 = u1_v->get<DiscreteVectorFunction>();
const DiscreteScalarFunction& p1 = p1_v->get<DiscreteScalarFunction>();
- const MeshType& mesh2 = dynamic_cast<const MeshType&>(*i_mesh2);
+ const MeshType& mesh2 = *i_mesh2->get<MeshType>();
const DiscreteScalarFunction& rho2 = rho2_v->get<DiscreteScalarFunction>();
const DiscreteScalarFunction& c2 = c2_v->get<DiscreteScalarFunction>();
const DiscreteVectorFunction& u2 = u2_v->get<DiscreteVectorFunction>();
@@ -668,24 +665,22 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
synchronize(Ar2);
synchronize(br2);
- NodeValue<Rd> ur1 = this->_computeUr(mesh1, Ar1, br1);
- NodeValue<Rd> ur2 = this->_computeUr(mesh2, Ar2, br2);
- this->_applyCouplingBC2(Ar1,Ar2,ur1,ur2,map1,map2);
+ NodeValue<Rd> ur1 = this->_computeUr(mesh1, Ar1, br1);
+ NodeValue<Rd> ur2 = this->_computeUr(mesh2, Ar2, br2);
+ this->_applyCouplingBC2(Ar1, Ar2, ur1, ur2, map1, map2);
NodeValuePerCell<Rd> Fjr1 = this->_computeFjr(mesh1, Ajr1, ur1, u1, p1);
NodeValuePerCell<Rd> Fjr2 = this->_computeFjr(mesh2, Ajr2, ur2, u2, p2);
return std::make_tuple(std::make_shared<const ItemValueVariant>(ur1),
std::make_shared<const SubItemValuePerItemVariant>(Fjr1),
- std::make_shared<const ItemValueVariant>(ur2),
- std::make_shared<const SubItemValuePerItemVariant>(Fjr2),
- ur2,
- Fjr2);
+ std::make_shared<const ItemValueVariant>(ur2),
+ std::make_shared<const SubItemValuePerItemVariant>(Fjr2), ur2, Fjr2);
}
- std::tuple<std::shared_ptr<const IMesh>,
+ 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>>
apply_fluxes(const double& dt,
const MeshType& mesh,
const DiscreteScalarFunction& rho,
@@ -734,12 +729,13 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { new_rho[j] *= Vj[j] / new_Vj[j]; });
- return {new_mesh, std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_rho)),
+ 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::tuple<std::shared_ptr<const IMesh>,
+ std::tuple<std::shared_ptr<const MeshVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
@@ -750,8 +746,8 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
const std::shared_ptr<const ItemValueVariant>& ur,
const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr) const
{
- std::shared_ptr i_mesh = getCommonMesh({rho_v, u_v, E_v});
- if (not i_mesh) {
+ std::shared_ptr mesh_v = getCommonMesh({rho_v, u_v, E_v});
+ if (not mesh_v) {
throw NormalError("discrete functions are not defined on the same mesh");
}
@@ -759,36 +755,35 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
throw NormalError("acoustic solver expects P0 functions");
}
-
- return this->apply_fluxes(dt, //
- dynamic_cast<const MeshType&>(*i_mesh), //
- rho_v->get<DiscreteScalarFunction>(), //
- u_v->get<DiscreteVectorFunction>(), //
- E_v->get<DiscreteScalarFunction>(), //
- ur->get<NodeValue<const Rd>>(), //
+ return this->apply_fluxes(dt, //
+ *mesh_v->get<MeshType>(), //
+ rho_v->get<DiscreteScalarFunction>(), //
+ u_v->get<DiscreteVectorFunction>(), //
+ E_v->get<DiscreteScalarFunction>(), //
+ ur->get<NodeValue<const Rd>>(), //
Fjr->get<NodeValuePerCell<const Rd>>());
}
- std::tuple<std::shared_ptr<const IMesh>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>>
+ std::tuple<std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
mesh_correction(const MeshType& mesh1,
- const MeshType& mesh2,
- const DiscreteScalarFunction& rho,
- const DiscreteVectorFunction& u,
- const DiscreteScalarFunction& E,
- std::vector<NodeId>& map1,
- std::vector<NodeId>& map2) const
- {
+ const MeshType& mesh2,
+ const DiscreteScalarFunction& rho,
+ const DiscreteVectorFunction& u,
+ const DiscreteScalarFunction& E,
+ std::vector<NodeId>& map1,
+ std::vector<NodeId>& map2) const
+ {
NodeValue<Rd> xr1 = copy(mesh1.xr());
NodeValue<Rd> new_xr2 = copy(mesh2.xr());
size_t n = map1.size();
- for(size_t i = 0; i<n; i++){
- for(size_t j = 0; j<Dimension; j++){
- new_xr2[map2[i]][j] = xr1[map1[i]][j];
+ for (size_t i = 0; i < n; i++) {
+ for (size_t j = 0; j < Dimension; j++) {
+ new_xr2[map2[i]][j] = xr1[map1[i]][j];
}
}
@@ -798,72 +793,69 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
CellValue<Rd> new_u = copy(u.cellValues());
CellValue<double> new_E = copy(E.cellValues());
- return {new_mesh2,std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh2, new_rho)),
+ return {std::make_shared<MeshVariant>(new_mesh2),
+ std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh2, new_rho)),
std::make_shared<DiscreteFunctionVariant>(DiscreteVectorFunction(new_mesh2, new_u)),
std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh2, new_E))};
}
- std::tuple<std::shared_ptr<const IMesh>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>>
- mesh_correction(const std::shared_ptr<const IMesh>& i_mesh1,
- const std::shared_ptr<const IMesh>& i_mesh2,
- const std::shared_ptr<const DiscreteFunctionVariant>& rho,
- const std::shared_ptr<const DiscreteFunctionVariant>& u,
- const std::shared_ptr<const DiscreteFunctionVariant>& E,
- std::vector<NodeId>& map1,
- std::vector<NodeId>& map2) const
+ std::tuple<std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+ mesh_correction(const std::shared_ptr<const MeshVariant>& i_mesh1,
+ const std::shared_ptr<const MeshVariant>& i_mesh2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E,
+ std::vector<NodeId>& map1,
+ std::vector<NodeId>& map2) const
{
-
- return this->mesh_correction(dynamic_cast<const MeshType&>(*i_mesh1),
- dynamic_cast<const MeshType&>(*i_mesh2),
- rho->get<DiscreteScalarFunction>(),
- u->get<DiscreteVectorFunction>(),
- E->get<DiscreteScalarFunction>(),
- map1,
- map2);
+ return this->mesh_correction(*i_mesh1->get<MeshType>(), *i_mesh2->get<MeshType>(),
+ rho->get<DiscreteScalarFunction>(), u->get<DiscreteVectorFunction>(),
+ E->get<DiscreteScalarFunction>(), map1, map2);
}
- std::tuple<std::shared_ptr<const IMesh>,
+ 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 MeshVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const IMesh>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>>
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
apply(const SolverType& solver_type,
const double& dt1,
- const size_t& q,
+ const size_t& q,
const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
- const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
const std::shared_ptr<const DiscreteFunctionVariant>& u1,
- const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u2,
const std::shared_ptr<const DiscreteFunctionVariant>& E1,
- const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E2,
const std::shared_ptr<const DiscreteFunctionVariant>& c1,
- const std::shared_ptr<const DiscreteFunctionVariant>& c2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c2,
const std::shared_ptr<const DiscreteFunctionVariant>& p1,
- const std::shared_ptr<const DiscreteFunctionVariant>& p2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p2,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2) const
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2) const
{
- std::shared_ptr<const IMesh> new_m2 = getCommonMesh({rho2, c2, u2, p2});
- std::shared_ptr<const IMesh> m1 = getCommonMesh({rho1, c1, u1, p1});
-
+ std::shared_ptr new_m2 = getCommonMesh({rho2, c2, u2, p2});
+ std::shared_ptr m1 = getCommonMesh({rho1, c1, u1, p1});
+
std::shared_ptr<const DiscreteFunctionVariant> new_rho2 = rho2;
std::shared_ptr<const DiscreteFunctionVariant> new_u2 = u2;
std::shared_ptr<const DiscreteFunctionVariant> new_E2 = E2;
- auto [map1, map2] = _computeMapping(m1,new_m2,bc_descriptor_list1,bc_descriptor_list2);
-
- auto [ur1, Fjr1, ur2, Fjr2,CR_ur,CR_Fjr] = compute_fluxes(solver_type, rho1, c1, u1, p1, bc_descriptor_list1, rho2, c2, u2, p2, bc_descriptor_list2,map1,map2);
- const auto [new_m1,new_rho1,new_u1,new_E1] = apply_fluxes(dt1, rho1, u1, E1, ur1, Fjr1);
-
- double dt2 = dt1/q;
+ auto [map1, map2] = _computeMapping(m1, new_m2, bc_descriptor_list1, bc_descriptor_list2);
+
+ auto [ur1, Fjr1, ur2, Fjr2, CR_ur, CR_Fjr] = compute_fluxes(solver_type, rho1, c1, u1, p1, bc_descriptor_list1,
+ rho2, c2, u2, p2, bc_descriptor_list2, map1, map2);
+ const auto [new_m1, new_rho1, new_u1, new_E1] = apply_fluxes(dt1, rho1, u1, E1, ur1, Fjr1);
+
+ double dt2 = dt1 / q;
const double& gamma = 1.4;
- double sum_dt = 0;
+ double sum_dt = 0;
// do{
@@ -878,69 +870,71 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
// const DiscreteScalarFunction& p_d = (gamma-1) * rho_d * eps;
// const DiscreteScalarFunction& c_d = sqrt(gamma * p_d / rho_d);
- // const std::shared_ptr<const DiscreteFunctionVariant>& new_p2 = std::make_shared<const DiscreteFunctionVariant>(p_d);
- // const std::shared_ptr<const DiscreteFunctionVariant>& new_c2 = std::make_shared<const DiscreteFunctionVariant>(c_d);
-
- // auto [ur2,Fjr2] = compute_fluxes(solver_type,new_rho2,new_c2,new_u2,new_p2,bc_descriptor_list2,CR_ur,CR_Fjr,map2);
+ // const std::shared_ptr<const DiscreteFunctionVariant>& new_p2 = std::make_shared<const
+ // DiscreteFunctionVariant>(p_d); const std::shared_ptr<const DiscreteFunctionVariant>& new_c2 =
+ // std::make_shared<const DiscreteFunctionVariant>(c_d);
+
+ // auto [ur2,Fjr2] =
+ // compute_fluxes(solver_type,new_rho2,new_c2,new_u2,new_p2,bc_descriptor_list2,CR_ur,CR_Fjr,map2);
// std::tie(new_m2,new_rho2,new_u2,new_E2) = apply_fluxes(dt2,new_rho2,new_u2,new_E2,ur2,Fjr2);
// sum_dt += dt2;
// }while(sum_dt < dt1);
- for(size_t i = 0; i<q; i++){
-
- if(i == q-1){
- dt2 = dt1 - sum_dt;
+ for (size_t i = 0; i < q; i++) {
+ if (i == q - 1) {
+ dt2 = dt1 - sum_dt;
}
const DiscreteScalarFunction& rho_d = new_rho2->get<DiscreteScalarFunction>();
const DiscreteVectorFunction& u_d = new_u2->get<DiscreteVectorFunction>();
const DiscreteScalarFunction& E_d = new_E2->get<DiscreteScalarFunction>();
- const DiscreteScalarFunction& eps = E_d - 0.5 * dot(u_d,u_d);
- const DiscreteScalarFunction& p_d = (gamma-1) * rho_d * eps;
+ const DiscreteScalarFunction& eps = E_d - 0.5 * dot(u_d, u_d);
+ const DiscreteScalarFunction& p_d = (gamma - 1) * rho_d * eps;
const DiscreteScalarFunction& c_d = sqrt(gamma * p_d / rho_d);
- const std::shared_ptr<const DiscreteFunctionVariant>& new_p2 = std::make_shared<const DiscreteFunctionVariant>(p_d);
- const std::shared_ptr<const DiscreteFunctionVariant>& new_c2 = std::make_shared<const DiscreteFunctionVariant>(c_d);
-
- auto [ur2,Fjr2] = compute_fluxes(solver_type,new_rho2,new_c2,new_u2,new_p2,bc_descriptor_list2,CR_ur,CR_Fjr,map2);
+ const std::shared_ptr<const DiscreteFunctionVariant>& new_p2 =
+ std::make_shared<const DiscreteFunctionVariant>(p_d);
+ const std::shared_ptr<const DiscreteFunctionVariant>& new_c2 =
+ std::make_shared<const DiscreteFunctionVariant>(c_d);
- std::tie(new_m2,new_rho2,new_u2,new_E2) = apply_fluxes(dt2,new_rho2,new_u2,new_E2,ur2,Fjr2);
+ auto [ur2, Fjr2] =
+ compute_fluxes(solver_type, new_rho2, new_c2, new_u2, new_p2, bc_descriptor_list2, CR_ur, CR_Fjr, map2);
+
+ std::tie(new_m2, new_rho2, new_u2, new_E2) = apply_fluxes(dt2, new_rho2, new_u2, new_E2, ur2, Fjr2);
sum_dt += dt2;
}
- std::tie(new_m2,new_rho2,new_u2,new_E2) = mesh_correction(new_m1,new_m2,new_rho2,new_u2,new_E2,map1,map2);
-
- return {new_m1,new_rho1,new_u1,new_E1,new_m2,new_rho2,new_u2,new_E2};
+ std::tie(new_m2, new_rho2, new_u2, new_E2) = mesh_correction(new_m1, new_m2, new_rho2, new_u2, new_E2, map1, map2);
+
+ return {new_m1, new_rho1, new_u1, new_E1, new_m2, new_rho2, new_u2, new_E2};
}
- std::tuple<std::shared_ptr<const IMesh>,
+ 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 IMesh>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>>
+ std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
apply(const SolverType& solver_type,
- const double& dt1,
+ const double& dt1,
const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
- const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
const std::shared_ptr<const DiscreteFunctionVariant>& u1,
- const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u2,
const std::shared_ptr<const DiscreteFunctionVariant>& E1,
- const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E2,
const std::shared_ptr<const DiscreteFunctionVariant>& c1,
- const std::shared_ptr<const DiscreteFunctionVariant>& c2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c2,
const std::shared_ptr<const DiscreteFunctionVariant>& p1,
- const std::shared_ptr<const DiscreteFunctionVariant>& p2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p2,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2) const
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2) const
{
-
-
const double& gamma = 1.4;
std::shared_ptr<const DiscreteFunctionVariant> new_rho2 = rho2;
@@ -948,56 +942,58 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
std::shared_ptr<const DiscreteFunctionVariant> new_E2 = E2;
std::shared_ptr<const DiscreteFunctionVariant> new_c2 = c2;
std::shared_ptr<const DiscreteFunctionVariant> new_p2 = p2;
- std::shared_ptr<const IMesh> new_m2 = getCommonMesh({new_rho2, new_c2, new_u2, new_p2});
- std::shared_ptr<const IMesh> m1 = getCommonMesh({rho1, c1, u1, p1});
+ std::shared_ptr new_m2 = getCommonMesh({new_rho2, new_c2, new_u2, new_p2});
+ std::shared_ptr m1 = getCommonMesh({rho1, c1, u1, p1});
double dt2 = 0.4 * acoustic_dt(new_c2);
- auto [map1, map2] = _computeMapping(m1,new_m2,bc_descriptor_list1,bc_descriptor_list2);
+ auto [map1, map2] = _computeMapping(m1, new_m2, bc_descriptor_list1, bc_descriptor_list2);
- auto [ur1, Fjr1, ur2, Fjr2,CR_ur,CR_Fjr] = compute_fluxes2(solver_type, rho1, c1, u1, p1, bc_descriptor_list1, new_rho2, new_c2, new_u2, new_p2, bc_descriptor_list2,map1,map2);
+ auto [ur1, Fjr1, ur2, Fjr2, CR_ur, CR_Fjr] =
+ compute_fluxes2(solver_type, rho1, c1, u1, p1, bc_descriptor_list1, new_rho2, new_c2, new_u2, new_p2,
+ bc_descriptor_list2, map1, map2);
- const auto [new_m1,new_rho1,new_u1,new_E1] = apply_fluxes(dt1, rho1, u1, E1, ur1, Fjr1);
- std::tie(new_m2,new_rho2,new_u2,new_E2) = apply_fluxes(dt2,new_rho2,new_u2,new_E2,ur2,Fjr2);
-
- double sum_dt = dt2;
+ const auto [new_m1, new_rho1, new_u1, new_E1] = apply_fluxes(dt1, rho1, u1, E1, ur1, Fjr1);
+ std::tie(new_m2, new_rho2, new_u2, new_E2) = apply_fluxes(dt2, new_rho2, new_u2, new_E2, ur2, Fjr2);
- while(sum_dt < dt1){
+ double sum_dt = dt2;
+ while (sum_dt < dt1) {
const DiscreteScalarFunction& rho_d = new_rho2->get<DiscreteScalarFunction>();
const DiscreteVectorFunction& u_d = new_u2->get<DiscreteVectorFunction>();
const DiscreteScalarFunction& E_d = new_E2->get<DiscreteScalarFunction>();
- const DiscreteScalarFunction& eps = E_d - 0.5 * dot(u_d,u_d);
- const DiscreteScalarFunction& p_d = (gamma-1) * rho_d * eps;
+ const DiscreteScalarFunction& eps = E_d - 0.5 * dot(u_d, u_d);
+ const DiscreteScalarFunction& p_d = (gamma - 1) * rho_d * eps;
const DiscreteScalarFunction& c_d = sqrt(gamma * p_d / rho_d);
new_p2 = std::make_shared<DiscreteFunctionVariant>(p_d);
new_c2 = std::make_shared<DiscreteFunctionVariant>(c_d);
double dt2 = 0.4 * acoustic_dt(new_c2);
- //std::cout << "dt2 = " << dt2 << std::endl;
- if(sum_dt + dt2 > dt1){
- dt2 = dt1 - sum_dt;
+ // std::cout << "dt2 = " << dt2 << std::endl;
+ if (sum_dt + dt2 > dt1) {
+ dt2 = dt1 - sum_dt;
}
- auto [ur2,Fjr2] = compute_fluxes(solver_type,new_rho2,new_c2,new_u2,new_p2,bc_descriptor_list2,CR_ur,CR_Fjr,map2);
- std::tie(new_m2,new_rho2,new_u2,new_E2) = apply_fluxes(dt2,new_rho2,new_u2,new_E2,ur2,Fjr2);
+ auto [ur2, Fjr2] =
+ compute_fluxes(solver_type, new_rho2, new_c2, new_u2, new_p2, bc_descriptor_list2, CR_ur, CR_Fjr, map2);
+ std::tie(new_m2, new_rho2, new_u2, new_E2) = apply_fluxes(dt2, new_rho2, new_u2, new_E2, ur2, Fjr2);
sum_dt += dt2;
}
- std::tie(new_m2,new_rho2,new_u2,new_E2) = mesh_correction(new_m1,new_m2,new_rho2,new_u2,new_E2,map1,map2);
- return {new_m1,new_rho1,new_u1,new_E1,new_m2,new_rho2,new_u2,new_E2};
+ std::tie(new_m2, new_rho2, new_u2, new_E2) = mesh_correction(new_m1, new_m2, new_rho2, new_u2, new_E2, map1, map2);
+ return {new_m1, new_rho1, new_u1, new_E1, new_m2, new_rho2, new_u2, new_E2};
}
//************************************************
//**** Nouvelles fcts pour l'ordre 2 en temps ****
//************************************************
- std::tuple<std::shared_ptr<const IMesh>,
+ 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>>
apply_fluxes(const double& dt,
- const double& tk,
+ const double& tk,
const MeshType& mesh,
const DiscreteScalarFunction& rho,
const DiscreteVectorFunction& u,
@@ -1011,31 +1007,28 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
if ((mesh.shared_connectivity() != ur_n.connectivity_ptr()) or
(mesh.shared_connectivity() != Fjr_n.connectivity_ptr())) {
- // throw NormalError("fluxes are not defined on the same connectivity than the mesh");
+ // throw NormalError("fluxes are not defined on the same connectivity than the mesh");
}
NodeValue<Rd> ur{mesh.connectivity()};
parallel_for(
- mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r){
-
- for(size_t i = 0; i<Dimension; ++i){
- ur[r][i] = (1-tk)*ur_n[r][i] + tk*ur_app[r][i];
+ mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+ for (size_t i = 0; i < Dimension; ++i) {
+ ur[r][i] = (1 - tk) * ur_n[r][i] + tk * ur_app[r][i];
}
- }
- );
+ });
NodeValuePerCell<Rd> Fjr{mesh.connectivity()};
parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(CellId j){
+ mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
const auto& cell_nodes = cell_to_node_matrix[j];
- for (size_t r = 0; r < cell_nodes.size(); ++r){
- for(size_t i = 0; i<Dimension; i++){
- Fjr(j,r)[i] = (1-tk)*Fjr_n(j,r)[i] + tk*Fjr_app(j,r)[i];
+ for (size_t r = 0; r < cell_nodes.size(); ++r) {
+ for (size_t i = 0; i < Dimension; i++) {
+ Fjr(j, r)[i] = (1 - tk) * Fjr_n(j, r)[i] + tk * Fjr_app(j, r)[i];
}
}
- }
- );
+ });
NodeValue<Rd> new_xr = copy(mesh.xr());
parallel_for(
@@ -1054,18 +1047,18 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
const auto& cell_nodes = cell_to_node_matrix[j];
Rd momentum_fluxes = zero;
- double energy_fluxes = 0 ;
- //double energy_fluxes_n = 0;
- //double energy_fluxes_app = 0;
+ double energy_fluxes = 0;
+ // double energy_fluxes_n = 0;
+ // double energy_fluxes_app = 0;
for (size_t R = 0; R < cell_nodes.size(); ++R) {
const NodeId r = cell_nodes[R];
momentum_fluxes += Fjr(j, R);
- energy_fluxes += ((1-tk)*dot(Fjr_n(j, R), ur_n[r]) + tk*dot(Fjr_app(j,R),ur_app[r]));
- //energy_fluxes_n += dot(Fjr_n(j, R), ur_n[r]);
- //energy_fluxes_app += dot(Fjr_app(j, R), ur_app[r]);
+ energy_fluxes += ((1 - tk) * dot(Fjr_n(j, R), ur_n[r]) + tk * dot(Fjr_app(j, R), ur_app[r]));
+ // energy_fluxes_n += dot(Fjr_n(j, R), ur_n[r]);
+ // energy_fluxes_app += dot(Fjr_app(j, R), ur_app[r]);
}
- //double energy_fluxes = (1-tk)*energy_fluxes_n + tk*energy_fluxes_app;
- const double dt_over_Mj = dt / (rho[j] * Vj[j]);
+ // double energy_fluxes = (1-tk)*energy_fluxes_n + tk*energy_fluxes_app;
+ 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;
});
@@ -1075,15 +1068,16 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { new_rho[j] *= Vj[j] / new_Vj[j]; });
- return {new_mesh, std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_rho)),
+ 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::tuple<std::shared_ptr<const IMesh>,
+ 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>>
apply_fluxes(const double& dt,
const MeshType& mesh,
const DiscreteScalarFunction& rho,
@@ -1098,26 +1092,22 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
if ((mesh.shared_connectivity() != ur_n.connectivity_ptr()) or
(mesh.shared_connectivity() != Fjr_n.connectivity_ptr())) {
- // throw NormalError("fluxes are not defined on the same connectivity than the mesh");
+ // throw NormalError("fluxes are not defined on the same connectivity than the mesh");
}
NodeValue<Rd> ur{mesh.connectivity()};
parallel_for(
- mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r){
- ur[r] = 0.5*(ur_n[r] + ur_app[r]);
- }
- );
+ mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { ur[r] = 0.5 * (ur_n[r] + ur_app[r]); });
NodeValuePerCell<Rd> Fjr{mesh.connectivity()};
parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(CellId j){
+ mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
const auto& cell_nodes = cell_to_node_matrix[j];
- for (size_t r = 0; r < cell_nodes.size(); ++r){
- Fjr(j,r) = 0.5*(Fjr_n(j,r) + Fjr_app(j,r));
+ for (size_t r = 0; r < cell_nodes.size(); ++r) {
+ Fjr(j, r) = 0.5 * (Fjr_n(j, r) + Fjr_app(j, r));
}
- }
- );
+ });
NodeValue<Rd> new_xr = copy(mesh.xr());
parallel_for(
@@ -1136,13 +1126,13 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
const auto& cell_nodes = cell_to_node_matrix[j];
Rd momentum_fluxes = zero;
- double energy_fluxes = 0 ;
+ double energy_fluxes = 0;
for (size_t R = 0; R < cell_nodes.size(); ++R) {
const NodeId r = cell_nodes[R];
momentum_fluxes += Fjr(j, R);
- energy_fluxes += 0.5*dot(Fjr_n(j,R),ur_n[r]) + 0.5*dot(Fjr_app(j,R),ur_app[r]);
+ energy_fluxes += 0.5 * dot(Fjr_n(j, R), ur_n[r]) + 0.5 * dot(Fjr_app(j, R), ur_app[r]);
}
- const double dt_over_Mj = dt / (rho[j] * Vj[j]);
+ 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;
});
@@ -1152,12 +1142,13 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { new_rho[j] *= Vj[j] / new_Vj[j]; });
- return {new_mesh, std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_rho)),
+ 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::tuple<std::shared_ptr<const IMesh>,
+ std::tuple<std::shared_ptr<const MeshVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
@@ -1180,20 +1171,19 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
throw NormalError("acoustic solver expects P0 functions");
}
-
return this->apply_fluxes(dt,
- tk, //
- dynamic_cast<const MeshType&>(*i_mesh), //
- rho_v->get<DiscreteScalarFunction>(), //
- u_v->get<DiscreteVectorFunction>(), //
- E_v->get<DiscreteScalarFunction>(), //
- ur_n->get<NodeValue<const Rd>>(), //
+ tk, //
+ *i_mesh->get<MeshType>(), //
+ rho_v->get<DiscreteScalarFunction>(), //
+ u_v->get<DiscreteVectorFunction>(), //
+ E_v->get<DiscreteScalarFunction>(), //
+ ur_n->get<NodeValue<const Rd>>(), //
Fjr_n->get<NodeValuePerCell<const Rd>>(),
- ur_app->get<NodeValue<const Rd>>(), //
+ ur_app->get<NodeValue<const Rd>>(), //
Fjr_app->get<NodeValuePerCell<const Rd>>());
}
- std::tuple<std::shared_ptr<const IMesh>,
+ std::tuple<std::shared_ptr<const MeshVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
@@ -1215,33 +1205,32 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
throw NormalError("acoustic solver expects P0 functions");
}
-
- return this->apply_fluxes(dt, //
- dynamic_cast<const MeshType&>(*i_mesh), //
- rho_v->get<DiscreteScalarFunction>(), //
- u_v->get<DiscreteVectorFunction>(), //
- E_v->get<DiscreteScalarFunction>(), //
- ur_n->get<NodeValue<const Rd>>(), //
+ return this->apply_fluxes(dt, //
+ *i_mesh->get<MeshType>(), //
+ rho_v->get<DiscreteScalarFunction>(), //
+ u_v->get<DiscreteVectorFunction>(), //
+ E_v->get<DiscreteScalarFunction>(), //
+ ur_n->get<NodeValue<const Rd>>(), //
Fjr_n->get<NodeValuePerCell<const Rd>>(),
- ur_app->get<NodeValue<const Rd>>(), //
+ ur_app->get<NodeValue<const Rd>>(), //
Fjr_app->get<NodeValuePerCell<const Rd>>());
}
std::tuple<const std::shared_ptr<const ItemValueVariant>,
- const std::shared_ptr<const SubItemValuePerItemVariant>,
- const std::shared_ptr<const ItemValueVariant>,
- const std::shared_ptr<const SubItemValuePerItemVariant>>
+ const std::shared_ptr<const SubItemValuePerItemVariant>,
+ const std::shared_ptr<const ItemValueVariant>,
+ const std::shared_ptr<const SubItemValuePerItemVariant>>
compute_fluxes2(const SolverType& solver_type,
- const std::shared_ptr<const DiscreteFunctionVariant>& rho1_v,
- const std::shared_ptr<const DiscreteFunctionVariant>& c1_v,
- const std::shared_ptr<const DiscreteFunctionVariant>& u1_v,
- const std::shared_ptr<const DiscreteFunctionVariant>& p1_v,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
- const std::shared_ptr<const DiscreteFunctionVariant>& rho2_v,
- const std::shared_ptr<const DiscreteFunctionVariant>& c2_v,
- const std::shared_ptr<const DiscreteFunctionVariant>& u2_v,
- const std::shared_ptr<const DiscreteFunctionVariant>& p2_v,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2) const
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho1_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c1_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u1_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p1_v,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho2_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c2_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u2_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p2_v,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2) const
{
std::shared_ptr i_mesh1 = getCommonMesh({rho1_v, c1_v, u1_v, p1_v});
std::shared_ptr i_mesh2 = getCommonMesh({rho2_v, c2_v, u2_v, p2_v});
@@ -1261,7 +1250,7 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
throw NormalError("acoustic solver expects P0 functions");
}
- auto [map1,map2] = _computeMapping(i_mesh1,i_mesh2,bc_descriptor_list1,bc_descriptor_list2);
+ auto [map1, map2] = _computeMapping(i_mesh1, i_mesh2, bc_descriptor_list1, bc_descriptor_list2);
const MeshType& mesh1 = dynamic_cast<const MeshType&>(*i_mesh1);
const DiscreteScalarFunction& rho1 = rho1_v->get<DiscreteScalarFunction>();
@@ -1295,15 +1284,15 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
synchronize(Ar2);
synchronize(br2);
- NodeValue<Rd> ur1 = this->_computeUr(mesh1, Ar1, br1);
- NodeValue<Rd> ur2 = this->_computeUr(mesh2, Ar2, br2);
- this->_applyCouplingBC2(Ar1,Ar2,ur1,ur2,map1,map2);
+ NodeValue<Rd> ur1 = this->_computeUr(mesh1, Ar1, br1);
+ NodeValue<Rd> ur2 = this->_computeUr(mesh2, Ar2, br2);
+ this->_applyCouplingBC2(Ar1, Ar2, ur1, ur2, map1, map2);
NodeValuePerCell<Rd> Fjr1 = this->_computeFjr(mesh1, Ajr1, ur1, u1, p1);
NodeValuePerCell<Rd> Fjr2 = this->_computeFjr(mesh2, Ajr2, ur2, u2, p2);
return std::make_tuple(std::make_shared<const ItemValueVariant>(ur1),
std::make_shared<const SubItemValuePerItemVariant>(Fjr1),
- std::make_shared<const ItemValueVariant>(ur2),
+ std::make_shared<const ItemValueVariant>(ur2),
std::make_shared<const SubItemValuePerItemVariant>(Fjr2));
}
@@ -1314,8 +1303,8 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
const std::shared_ptr<const DiscreteFunctionVariant>& p_v,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
- const std::shared_ptr<const ItemValueVariant>& CR_ur,
- const std::shared_ptr<const SubItemValuePerItemVariant>& CR_Fjr) const
+ const std::shared_ptr<const ItemValueVariant>& CR_ur,
+ const std::shared_ptr<const SubItemValuePerItemVariant>& CR_Fjr) const
{
std::shared_ptr i_mesh = getCommonMesh({rho_v, c_v, u_v, p_v});
if (not i_mesh) {
@@ -1333,25 +1322,27 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
const DiscreteScalarFunction& p = p_v->get<DiscreteScalarFunction>();
std::vector<NodeId> map2;
- const BoundaryConditionList bc_list2 = this->_getBCList(mesh,bc_descriptor_list);
-
- for(const auto& boundary_condition2 : bc_list2){
- std::visit([&](auto&& bc2){
- using T2 = std::decay_t<decltype(bc2)>;
- if constexpr (std::is_same_v<CouplingBoundaryCondition,T2>){
- const auto& node_list2 = bc2.nodeList();
- for(size_t i = 0; i<node_list2.size(); i++){
- map2.push_back(node_list2[i]);
+ const BoundaryConditionList bc_list2 = this->_getBCList(mesh, bc_descriptor_list);
+
+ for (const auto& boundary_condition2 : bc_list2) {
+ std::visit(
+ [&](auto&& bc2) {
+ using T2 = std::decay_t<decltype(bc2)>;
+ if constexpr (std::is_same_v<CouplingBoundaryCondition, T2>) {
+ const auto& node_list2 = bc2.nodeList();
+ for (size_t i = 0; i < node_list2.size(); i++) {
+ map2.push_back(node_list2[i]);
+ }
}
- }
- },boundary_condition2);
+ },
+ boundary_condition2);
}
NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * c);
NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
NodeValue<Rd> br = this->_computeBr(mesh, Ajr, u, p);
-
+
const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
this->_applyBoundaryConditions(bc_list, mesh, Ar, br);
this->_applyExternalVelocityBC(bc_list, p, Ar, br);
@@ -1362,60 +1353,57 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
NodeValue<Rd> ur = this->_computeUr(mesh, Ar, br);
NodeValuePerCell<Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u, p);
- NodeValue<const Rd> ur_c = CR_ur->get<NodeValue<const Rd>>();
- NodeValuePerCell<const Rd> Fjr_c = CR_Fjr->get<NodeValuePerCell<const Rd>>();
+ NodeValue<const Rd> ur_c = CR_ur->get<NodeValue<const Rd>>();
+ NodeValuePerCell<const Rd> Fjr_c = CR_Fjr->get<NodeValuePerCell<const Rd>>();
const size_t& n = map2.size();
- const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
- const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+ const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
+ const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
- for(size_t i = 0; i<n; i++){
-
- const NodeId node_id = map2[i];
- const auto& node_to_cell = node_to_cell_matrix[node_id];
- const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
+ for (size_t i = 0; i < n; i++) {
+ const NodeId node_id = map2[i];
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+ const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
- for(size_t j = 0; j < node_to_cell.size(); j++){
- const CellId J = node_to_cell[j];
- const unsigned int R = node_local_number_in_its_cells[j];
- Fjr(J,R) = Fjr_c(J,R);
+ for (size_t j = 0; j < node_to_cell.size(); j++) {
+ const CellId J = node_to_cell[j];
+ const unsigned int R = node_local_number_in_its_cells[j];
+ Fjr(J, R) = Fjr_c(J, R);
+ }
+ ur[node_id] = ur_c[node_id];
}
- ur[node_id] = ur_c[node_id];
- }
-
- //this->_applyCouplingBC(mesh,ur,ur_c,Fjr,Fjr_c,map2);
- //this->_applyCouplingBC(mesh,ur,CR_ur,Fjr,CR_Fjr,map2);
+
+ // this->_applyCouplingBC(mesh,ur,ur_c,Fjr,Fjr_c,map2);
+ // this->_applyCouplingBC(mesh,ur,CR_ur,Fjr,CR_Fjr,map2);
return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
std::make_shared<const SubItemValuePerItemVariant>(Fjr));
}
- std::tuple<std::shared_ptr<const IMesh>,
+ 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 MeshVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const IMesh>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>>
+ std::shared_ptr<const DiscreteFunctionVariant>>
applyOrder2(const SolverType& solver_type,
- const double& dt1,
- const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
- const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
- const std::shared_ptr<const DiscreteFunctionVariant>& u1,
- const std::shared_ptr<const DiscreteFunctionVariant>& u2,
- const std::shared_ptr<const DiscreteFunctionVariant>& E1,
- const std::shared_ptr<const DiscreteFunctionVariant>& E2,
- const std::shared_ptr<const DiscreteFunctionVariant>& c1,
- const std::shared_ptr<const DiscreteFunctionVariant>& c2,
- const std::shared_ptr<const DiscreteFunctionVariant>& p1,
- const std::shared_ptr<const DiscreteFunctionVariant>& p2,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2) const
+ const double& dt1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p2,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2) const
{
-
-
const double& gamma = 1.4;
std::shared_ptr<const DiscreteFunctionVariant> new_rho2 = rho2;
@@ -1423,20 +1411,21 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
std::shared_ptr<const DiscreteFunctionVariant> new_E2 = E2;
std::shared_ptr<const DiscreteFunctionVariant> new_c2 = c2;
std::shared_ptr<const DiscreteFunctionVariant> new_p2 = p2;
- std::shared_ptr<const IMesh> new_m2 = getCommonMesh({new_rho2, new_c2, new_u2, new_p2});
- std::shared_ptr<const IMesh> m1 = getCommonMesh({rho1, c1, u1, p1});
-
- const MeshType& mesh1 = dynamic_cast<const MeshType&>(*m1);
- const MeshType& mesh2 = dynamic_cast<const MeshType&>(*new_m2);
+ std::shared_ptr new_m2 = getCommonMesh({new_rho2, new_c2, new_u2, new_p2});
+ std::shared_ptr m1 = getCommonMesh({rho1, c1, u1, p1});
- //double dt2 = 0.4 * acoustic_dt(new_c2);
- //double dt2_k = dt1/2;
- double dt2_k = dt1/2;
+ const MeshType& mesh1 = *m1->get<MeshType>();
+ const MeshType& mesh2 = *new_m2->get<MeshType>();
- auto [map1, map2] = _computeMapping(m1,new_m2,bc_descriptor_list1,bc_descriptor_list2);
+ // double dt2 = 0.4 * acoustic_dt(new_c2);
+ // double dt2_k = dt1/2;
+ double dt2_k = dt1 / 2;
- const auto [ur1_n, Fjr1_n, ur2_n, Fjr2_n,CR_ur,CR_Fjr] = compute_fluxes2(solver_type, rho1, c1, u1, p1, bc_descriptor_list1, new_rho2, new_c2, new_u2, new_p2, bc_descriptor_list2,map1,map2);
+ auto [map1, map2] = _computeMapping(m1, new_m2, bc_descriptor_list1, bc_descriptor_list2);
+ const auto [ur1_n, Fjr1_n, ur2_n, Fjr2_n, CR_ur, CR_Fjr] =
+ compute_fluxes2(solver_type, rho1, c1, u1, p1, bc_descriptor_list1, new_rho2, new_c2, new_u2, new_p2,
+ bc_descriptor_list2, map1, map2);
std::shared_ptr<const DiscreteFunctionVariant> rho1_app;
std::shared_ptr<const DiscreteFunctionVariant> u1_app;
@@ -1444,49 +1433,49 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
std::shared_ptr<const DiscreteFunctionVariant> rho2_app;
std::shared_ptr<const DiscreteFunctionVariant> u2_app;
std::shared_ptr<const DiscreteFunctionVariant> E2_app;
- std::shared_ptr<const IMesh> m1_app;
- std::shared_ptr<const IMesh> m2_app;
+ std::shared_ptr<const MeshVariant> m1_app;
+ std::shared_ptr<const MeshVariant> m2_app;
- std::tie(m1_app,rho1_app,u1_app,E1_app) = apply_fluxes(dt1, rho1, u1, E1, ur1_n, Fjr1_n);
- std::tie(m2_app,rho2_app,u2_app,E2_app) = apply_fluxes(dt2_k, rho2, u2, E2, ur2_n, Fjr2_n);
+ std::tie(m1_app, rho1_app, u1_app, E1_app) = apply_fluxes(dt1, rho1, u1, E1, ur1_n, Fjr1_n);
+ std::tie(m2_app, rho2_app, u2_app, E2_app) = apply_fluxes(dt2_k, rho2, u2, E2, ur2_n, Fjr2_n);
double sum_dt = dt2_k;
- while(sum_dt < dt1){
-
+ while (sum_dt < dt1) {
const DiscreteScalarFunction& rho_d = rho2_app->get<DiscreteScalarFunction>();
const DiscreteVectorFunction& u_d = u2_app->get<DiscreteVectorFunction>();
const DiscreteScalarFunction& E_d = E2_app->get<DiscreteScalarFunction>();
- const DiscreteScalarFunction& eps = E_d - 0.5 * dot(u_d,u_d);
- const DiscreteScalarFunction& p_d = (gamma-1) * rho_d * eps;
+ const DiscreteScalarFunction& eps = E_d - 0.5 * dot(u_d, u_d);
+ const DiscreteScalarFunction& p_d = (gamma - 1) * rho_d * eps;
const DiscreteScalarFunction& c_d = sqrt(gamma * p_d / rho_d);
const std::shared_ptr<const DiscreteFunctionVariant> p2_k = std::make_shared<DiscreteFunctionVariant>(p_d);
const std::shared_ptr<const DiscreteFunctionVariant> c2_k = std::make_shared<DiscreteFunctionVariant>(c_d);
- //double dt2 = 0.4 * acoustic_dt(new_c2);
- dt2_k = dt1/2;
- //std::cout << "dt2 = " << dt2 << std::endl;
- if(sum_dt + dt2_k > dt1){
- dt2_k = dt1 - sum_dt;
+ // double dt2 = 0.4 * acoustic_dt(new_c2);
+ dt2_k = dt1 / 2;
+ // std::cout << "dt2 = " << dt2 << std::endl;
+ if (sum_dt + dt2_k > dt1) {
+ dt2_k = dt1 - sum_dt;
}
- auto [ur_k,Fjr_k] = compute_fluxes(solver_type,rho2_app,c2_k,u2_app,p2_k,bc_descriptor_list2,CR_ur,CR_Fjr,map2);
- std::tie(m2_app,rho2_app,u2_app,E2_app) = apply_fluxes(dt2_k,rho2_app,u2_app,E2_app,ur_k,Fjr_k);
+ auto [ur_k, Fjr_k] =
+ compute_fluxes(solver_type, rho2_app, c2_k, u2_app, p2_k, bc_descriptor_list2, CR_ur, CR_Fjr, map2);
+ std::tie(m2_app, rho2_app, u2_app, E2_app) = apply_fluxes(dt2_k, rho2_app, u2_app, E2_app, ur_k, Fjr_k);
sum_dt += dt2_k;
}
const DiscreteScalarFunction& rho1_d = rho1_app->get<DiscreteScalarFunction>();
const DiscreteVectorFunction& u1_d = u1_app->get<DiscreteVectorFunction>();
const DiscreteScalarFunction& E1_d = E1_app->get<DiscreteScalarFunction>();
- const DiscreteScalarFunction& eps1 = E1_d - 0.5 * dot(u1_d,u1_d);
- const DiscreteScalarFunction& p1_d = (gamma-1) * rho1_d * eps1;
+ const DiscreteScalarFunction& eps1 = E1_d - 0.5 * dot(u1_d, u1_d);
+ const DiscreteScalarFunction& p1_d = (gamma - 1) * rho1_d * eps1;
const DiscreteScalarFunction& c1_d = sqrt(gamma * p1_d / rho1_d);
const DiscreteScalarFunction& rho2_d = rho2_app->get<DiscreteScalarFunction>();
const DiscreteVectorFunction& u2_d = u2_app->get<DiscreteVectorFunction>();
const DiscreteScalarFunction& E2_d = E2_app->get<DiscreteScalarFunction>();
- const DiscreteScalarFunction& eps2 = E2_d - 0.5 * dot(u2_d,u2_d);
- const DiscreteScalarFunction& p2_d = (gamma-1) * rho2_d * eps2;
+ const DiscreteScalarFunction& eps2 = E2_d - 0.5 * dot(u2_d, u2_d);
+ const DiscreteScalarFunction& p2_d = (gamma - 1) * rho2_d * eps2;
const DiscreteScalarFunction& c2_d = sqrt(gamma * p2_d / rho2_d);
const std::shared_ptr<const DiscreteFunctionVariant> c1_app = std::make_shared<const DiscreteFunctionVariant>(c1_d);
@@ -1494,89 +1483,88 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
const std::shared_ptr<const DiscreteFunctionVariant> c2_app = std::make_shared<const DiscreteFunctionVariant>(c2_d);
const std::shared_ptr<const DiscreteFunctionVariant> p2_app = std::make_shared<const DiscreteFunctionVariant>(p2_d);
+ const auto [ur1_app, Fjr1_app, ur2_app, Fjr2_app, CRapp1, CRapp2] =
+ compute_fluxes2(solver_type, rho1_app, c1_app, u1_app, p1_app, bc_descriptor_list1, rho2_app, c2_app, u2_app,
+ p2_app, bc_descriptor_list2, map1, map2);
+ // double dt2 = dt1/2;
+ double dt2 = dt1 / 2;
+ double tk = (dt2 * 0.5) / dt1;
+ // double tk = 1;
+
+ // std::shared_ptr<const ItemValueVariant> ur_inter =_InterpolateUr(mesh2,ur2_n,ur2_app,tk);
+ // std::shared_ptr<const SubItemValuePerItemVariant> Fjr_inter= _InterpolateFjr(mesh2,Fjr2_n,Fjr2_app,tk);
+ // std::shared_ptr<const ItemValueVariant> ur2 = _computeOrder2Ur(mesh2,ur_inter,ur2_n);
+ // std::shared_ptr<const SubItemValuePerItemVariant> Fjr2 = _computeOrder2Fjr(mesh2,Fjr_inter,Fjr2_n);
+ // std::shared_ptr<const ItemValueVariant> ur1 = _computeOrder2Ur(mesh1,ur1_app,ur1_n);
+ // std::shared_ptr<const SubItemValuePerItemVariant> Fjr1 = _computeOrder2Fjr(mesh1,Fjr1_app,Fjr1_n);
+
+ const auto [new_m1, new_rho1, new_u1, new_E1] = apply_fluxes(dt1, rho1, u1, E1, ur1_n, Fjr1_n, ur1_app, Fjr1_app);
+ // const auto [new_m1,new_rho1,new_u1,new_E1] = apply_fluxes(dt1,rho1,u1,E1,ur1,Fjr1);
+ // std::tie(new_m2,new_rho2,new_u2,new_E2) = apply_fluxes(dt2,rho2,u2,E2,ur2,Fjr2);
+ std::tie(new_m2, new_rho2, new_u2, new_E2) = apply_fluxes(dt2, tk, rho2, u2, E2, ur2_n, Fjr2_n, ur2_app, Fjr2_app);
+ // std::tie(new_m2,new_rho2,new_u2,new_E2) =
+ // apply_fluxes(dt2,tk,rho2,u2,E2,ur2_n,Fjr2_n,ur2_app,Fjr2_app,ur2_n,Fjr2_n);
- const auto [ur1_app, Fjr1_app, ur2_app, Fjr2_app, CRapp1, CRapp2] = compute_fluxes2(solver_type, rho1_app, c1_app, u1_app, p1_app, bc_descriptor_list1, rho2_app, c2_app, u2_app, p2_app, bc_descriptor_list2,map1,map2);
- //double dt2 = dt1/2;
- double dt2 = dt1/2;
- double tk = (dt2*0.5)/dt1;
- //double tk = 1;
-
- //std::shared_ptr<const ItemValueVariant> ur_inter =_InterpolateUr(mesh2,ur2_n,ur2_app,tk);
- //std::shared_ptr<const SubItemValuePerItemVariant> Fjr_inter= _InterpolateFjr(mesh2,Fjr2_n,Fjr2_app,tk);
- //std::shared_ptr<const ItemValueVariant> ur2 = _computeOrder2Ur(mesh2,ur_inter,ur2_n);
- //std::shared_ptr<const SubItemValuePerItemVariant> Fjr2 = _computeOrder2Fjr(mesh2,Fjr_inter,Fjr2_n);
- //std::shared_ptr<const ItemValueVariant> ur1 = _computeOrder2Ur(mesh1,ur1_app,ur1_n);
- //std::shared_ptr<const SubItemValuePerItemVariant> Fjr1 = _computeOrder2Fjr(mesh1,Fjr1_app,Fjr1_n);
-
- const auto [new_m1,new_rho1,new_u1,new_E1] = apply_fluxes(dt1,rho1,u1,E1,ur1_n,Fjr1_n,ur1_app,Fjr1_app);
- //const auto [new_m1,new_rho1,new_u1,new_E1] = apply_fluxes(dt1,rho1,u1,E1,ur1,Fjr1);
- //std::tie(new_m2,new_rho2,new_u2,new_E2) = apply_fluxes(dt2,rho2,u2,E2,ur2,Fjr2);
- std::tie(new_m2,new_rho2,new_u2,new_E2) = apply_fluxes(dt2,tk,rho2,u2,E2,ur2_n,Fjr2_n,ur2_app,Fjr2_app);
- //std::tie(new_m2,new_rho2,new_u2,new_E2) = apply_fluxes(dt2,tk,rho2,u2,E2,ur2_n,Fjr2_n,ur2_app,Fjr2_app,ur2_n,Fjr2_n);
-
sum_dt = dt2;
- while(sum_dt < dt1){
-
+ while (sum_dt < dt1) {
const DiscreteScalarFunction& rho_d = new_rho2->get<DiscreteScalarFunction>();
const DiscreteVectorFunction& u_d = new_u2->get<DiscreteVectorFunction>();
const DiscreteScalarFunction& E_d = new_E2->get<DiscreteScalarFunction>();
- const DiscreteScalarFunction& eps = E_d - 0.5 * dot(u_d,u_d);
- const DiscreteScalarFunction& p_d = (gamma-1) * rho_d * eps;
+ const DiscreteScalarFunction& eps = E_d - 0.5 * dot(u_d, u_d);
+ const DiscreteScalarFunction& p_d = (gamma - 1) * rho_d * eps;
const DiscreteScalarFunction& c_d = sqrt(gamma * p_d / rho_d);
new_p2 = std::make_shared<DiscreteFunctionVariant>(p_d);
new_c2 = std::make_shared<DiscreteFunctionVariant>(c_d);
- //double dt2 = 0.4 * acoustic_dt(new_c2);
- dt2 = dt1/2;
- //dt2 = dt1;
- //std::cout << "boucle" << std::endl;
- if(sum_dt + dt2 > dt1){
- dt2 = dt1 - sum_dt;
+ // double dt2 = 0.4 * acoustic_dt(new_c2);
+ dt2 = dt1 / 2;
+ // dt2 = dt1;
+ // std::cout << "boucle" << std::endl;
+ if (sum_dt + dt2 > dt1) {
+ dt2 = dt1 - sum_dt;
}
- tk = (sum_dt + dt2*0.5)/dt1;
- //tk = 0.5*3;
-
- //const MeshType& new_mesh2 = dynamic_cast<const MeshType&>(*new_m2);
-//
-//
- //auto ur_inter = _InterpolateUr(new_mesh2,ur2_n,ur2_app,tk);
- //auto Fjr_inter = _InterpolateFjr(new_mesh2,Fjr2_n,Fjr2_app,tk);
- //auto [ur2_k,Fjr2_k] = compute_fluxes(solver_type,new_rho2,new_c2,new_u2,new_p2,bc_descriptor_list2,CR_ur,CR_Fjr,map2);
- //auto new_ur2 = _computeOrder2Ur(new_mesh2,ur_inter,ur2_k);
- //auto new_Fjr2 = _computeOrder2Fjr(new_mesh2,Fjr_inter,Fjr2_k);
-
- //std::tie(new_m2,new_rho2,new_u2,new_E2) = apply_fluxes(dt2,new_rho2,new_u2,new_E2,new_ur2,new_Fjr2);
- std::tie(new_m2,new_rho2,new_u2,new_E2) = apply_fluxes(dt2,tk,new_rho2,new_u2,new_E2,ur2_n,Fjr2_n,ur2_app,Fjr2_app);
+ tk = (sum_dt + dt2 * 0.5) / dt1;
+ // tk = 0.5*3;
+
+ // const MeshType& new_mesh2 = dynamic_cast<const MeshType&>(*new_m2);
+ //
+ //
+ // auto ur_inter = _InterpolateUr(new_mesh2,ur2_n,ur2_app,tk);
+ // auto Fjr_inter = _InterpolateFjr(new_mesh2,Fjr2_n,Fjr2_app,tk);
+ // auto [ur2_k,Fjr2_k] =
+ // compute_fluxes(solver_type,new_rho2,new_c2,new_u2,new_p2,bc_descriptor_list2,CR_ur,CR_Fjr,map2); auto new_ur2
+ // = _computeOrder2Ur(new_mesh2,ur_inter,ur2_k); auto new_Fjr2 =
+ // _computeOrder2Fjr(new_mesh2,Fjr_inter,Fjr2_k);
+
+ // std::tie(new_m2,new_rho2,new_u2,new_E2) = apply_fluxes(dt2,new_rho2,new_u2,new_E2,new_ur2,new_Fjr2);
+ std::tie(new_m2, new_rho2, new_u2, new_E2) =
+ apply_fluxes(dt2, tk, new_rho2, new_u2, new_E2, ur2_n, Fjr2_n, ur2_app, Fjr2_app);
sum_dt += dt2;
}
- std::tie(new_m2,new_rho2,new_u2,new_E2) = mesh_correction(new_m1,new_m2,new_rho2,new_u2,new_E2,map1,map2);
- return {new_m1,new_rho1,new_u1,new_E1,new_m2,new_rho2,new_u2,new_E2};
+ std::tie(new_m2, new_rho2, new_u2, new_E2) = mesh_correction(new_m1, new_m2, new_rho2, new_u2, new_E2, map1, map2);
+ return {new_m1, new_rho1, new_u1, new_E1, new_m2, new_rho2, new_u2, new_E2};
}
- //************************************************
- //************************************************ Fin
- //************************************************
-
LocalDtAcousticSolver() = default;
LocalDtAcousticSolver(LocalDtAcousticSolver&&) = default;
~LocalDtAcousticSolver() = default;
};
-template <size_t Dimension>
+template <MeshConcept MeshType>
void
-LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applyPressureBC(const BoundaryConditionList& bc_list,
- const MeshType& mesh,
- NodeValue<Rd>& br) const
+LocalDtAcousticSolverHandler::LocalDtAcousticSolver<MeshType>::_applyPressureBC(const BoundaryConditionList& bc_list,
+ const MeshType& mesh,
+ NodeValue<Rd>& br) const
{
for (const auto& boundary_condition : bc_list) {
std::visit(
[&](auto&& bc) {
using T = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<PressureBoundaryCondition, T>) {
- MeshData<Dimension>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+ MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
if constexpr (Dimension == 1) {
const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
@@ -1630,11 +1618,11 @@ LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applyPressureBC
}
}
-template <size_t Dimension>
+template <MeshConcept MeshType>
void
-LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applySymmetryBC(const BoundaryConditionList& bc_list,
- NodeValue<Rdxd>& Ar,
- NodeValue<Rd>& br) const
+LocalDtAcousticSolverHandler::LocalDtAcousticSolver<MeshType>::_applySymmetryBC(const BoundaryConditionList& bc_list,
+ NodeValue<Rdxd>& Ar,
+ NodeValue<Rd>& br) const
{
for (const auto& boundary_condition : bc_list) {
std::visit(
@@ -1661,11 +1649,11 @@ LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applySymmetryBC
}
}
-template <size_t Dimension>
+template <MeshConcept MeshType>
void
-LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applyVelocityBC(const BoundaryConditionList& bc_list,
- NodeValue<Rdxd>& Ar,
- NodeValue<Rd>& br) const
+LocalDtAcousticSolverHandler::LocalDtAcousticSolver<MeshType>::_applyVelocityBC(const BoundaryConditionList& bc_list,
+ NodeValue<Rdxd>& Ar,
+ NodeValue<Rd>& br) const
{
for (const auto& boundary_condition : bc_list) {
std::visit(
@@ -1698,9 +1686,9 @@ LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applyVelocityBC
}
}
-template <size_t Dimension>
+template <MeshConcept MeshType>
void
-LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applyExternalVelocityBC(
+LocalDtAcousticSolverHandler::LocalDtAcousticSolver<MeshType>::_applyExternalVelocityBC(
const BoundaryConditionList& bc_list,
const DiscreteScalarFunction& p,
NodeValue<Rdxd>& Ar,
@@ -1728,18 +1716,18 @@ LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applyExternalVe
}
}
-template <size_t Dimension>
-void LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applyCouplingBC(NodeValue<Rdxd>& Ar1,
- NodeValue<Rdxd>& Ar2,
- NodeValue<Rd>& br1,
- NodeValue<Rd>& br2,
- const std::vector<NodeId>& map1,
- const std::vector<NodeId>& map2) const
+template <MeshConcept MeshType>
+void
+LocalDtAcousticSolverHandler::LocalDtAcousticSolver<MeshType>::_applyCouplingBC(NodeValue<Rdxd>& Ar1,
+ NodeValue<Rdxd>& Ar2,
+ NodeValue<Rd>& br1,
+ NodeValue<Rd>& br2,
+ const std::vector<NodeId>& map1,
+ const std::vector<NodeId>& map2) const
{
size_t n = map1.size();
- for(size_t i = 0; i<n; i++){
-
+ for (size_t i = 0; i < n; i++) {
NodeId node_id1 = map1[i];
NodeId node_id2 = map2[i];
@@ -1748,66 +1736,65 @@ void LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applyCoupl
br1[node_id1] += br2[node_id2];
br2[node_id2] = br1[node_id1];
-
}
}
-template <size_t Dimension>
-void LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applyCouplingBC2(NodeValue<Rdxd>& Ar1,
- NodeValue<Rdxd>& Ar2,
- NodeValue<Rd>& ur1,
- NodeValue<Rd>& ur2,
- const std::vector<NodeId>& map1,
- const std::vector<NodeId>& map2) const
+template <MeshConcept MeshType>
+void
+LocalDtAcousticSolverHandler::LocalDtAcousticSolver<MeshType>::_applyCouplingBC2(NodeValue<Rdxd>& Ar1,
+ NodeValue<Rdxd>& Ar2,
+ NodeValue<Rd>& ur1,
+ NodeValue<Rd>& ur2,
+ const std::vector<NodeId>& map1,
+ const std::vector<NodeId>& map2) const
{
size_t n = map1.size();
Rd lambda;
- for(size_t i = 0; i<n; i++){
-
+ for (size_t i = 0; i < n; i++) {
NodeId node_id1 = map1[i];
NodeId node_id2 = map2[i];
-
- lambda = inverse(inverse(Ar2[node_id2]) + inverse(Ar1[node_id1]))*(ur2[node_id2]-ur1[node_id1]);
- ur1[node_id1] = ur1[node_id1] + inverse(Ar1[node_id1])*lambda;
- ur2[node_id2] = ur2[node_id2] - inverse(Ar2[node_id2])*lambda;
+ lambda = inverse(inverse(Ar2[node_id2]) + inverse(Ar1[node_id1])) * (ur2[node_id2] - ur1[node_id1]);
+
+ ur1[node_id1] = ur1[node_id1] + inverse(Ar1[node_id1]) * lambda;
+ ur2[node_id2] = ur2[node_id2] - inverse(Ar2[node_id2]) * lambda;
}
}
-template <size_t Dimension>
-void LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applyCouplingBC(const MeshType& mesh,
- NodeValue<Rd>& ur,
- NodeValue<Rd>& CR_ur,
- NodeValuePerCell<Rd>& Fjr,
- NodeValuePerCell<Rd>& CR_Fjr,
- const std::vector<NodeId>& map2) const
+template <MeshConcept MeshType>
+void
+LocalDtAcousticSolverHandler::LocalDtAcousticSolver<MeshType>::_applyCouplingBC(const MeshType& mesh,
+ NodeValue<Rd>& ur,
+ NodeValue<Rd>& CR_ur,
+ NodeValuePerCell<Rd>& Fjr,
+ NodeValuePerCell<Rd>& CR_Fjr,
+ const std::vector<NodeId>& map2) const
{
const size_t& n = map2.size();
const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
- for(size_t i = 0; i<n; i++){
-
+ for (size_t i = 0; i < n; i++) {
const NodeId node_id = map2[i];
const auto& node_to_cell = node_to_cell_matrix[node_id];
const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
- for(size_t j = 0; j < node_to_cell.size(); j++){
+ for (size_t j = 0; j < node_to_cell.size(); j++) {
const CellId J = node_to_cell[j];
const unsigned int R = node_local_number_in_its_cells[j];
- Fjr(J,R) = CR_Fjr(J,R);
+ Fjr(J, R) = CR_Fjr(J, R);
}
ur[node_id] = CR_ur[node_id];
}
}
-template <size_t Dimension>
-class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::FixedBoundaryCondition
+template <MeshConcept MeshType>
+class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<MeshType>::FixedBoundaryCondition
{
private:
- const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
+ const MeshNodeBoundary m_mesh_node_boundary;
public:
const Array<const NodeId>&
@@ -1816,18 +1803,16 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::FixedBound
return m_mesh_node_boundary.nodeList();
}
- FixedBoundaryCondition(const MeshNodeBoundary<Dimension> mesh_node_boundary)
- : m_mesh_node_boundary{mesh_node_boundary}
- {}
+ FixedBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary) : m_mesh_node_boundary{mesh_node_boundary} {}
~FixedBoundaryCondition() = default;
};
-template <size_t Dimension>
-class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::PressureBoundaryCondition
+template <MeshConcept MeshType>
+class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<MeshType>::PressureBoundaryCondition
{
private:
- const MeshFaceBoundary<Dimension> m_mesh_face_boundary;
+ const MeshFaceBoundary m_mesh_face_boundary;
const Array<const double> m_value_list;
public:
@@ -1843,8 +1828,7 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::PressureBo
return m_value_list;
}
- PressureBoundaryCondition(const MeshFaceBoundary<Dimension>& mesh_face_boundary,
- const Array<const double>& value_list)
+ PressureBoundaryCondition(const MeshFaceBoundary& mesh_face_boundary, const Array<const double>& value_list)
: m_mesh_face_boundary{mesh_face_boundary}, m_value_list{value_list}
{}
@@ -1852,10 +1836,10 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::PressureBo
};
template <>
-class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<1>::PressureBoundaryCondition
+class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Mesh<1>>::PressureBoundaryCondition
{
private:
- const MeshNodeBoundary<1> m_mesh_node_boundary;
+ const MeshNodeBoundary m_mesh_node_boundary;
const Array<const double> m_value_list;
public:
@@ -1871,19 +1855,19 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<1>::PressureBoundaryCo
return m_value_list;
}
- PressureBoundaryCondition(const MeshNodeBoundary<1>& mesh_node_boundary, const Array<const double>& value_list)
+ PressureBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const Array<const double>& value_list)
: m_mesh_node_boundary{mesh_node_boundary}, m_value_list{value_list}
{}
~PressureBoundaryCondition() = default;
};
-template <size_t Dimension>
-class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::ExternalFSIVelocityBoundaryCondition
+template <MeshConcept MeshType>
+class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<MeshType>::ExternalFSIVelocityBoundaryCondition
{
private:
const ItemToItemMatrix<ItemType::node, ItemType::cell> m_node_to_cell_matrix;
- const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
+ const MeshNodeBoundary m_mesh_node_boundary;
const Array<TinyVector<Dimension>> m_value_list;
const std::shared_ptr<const Socket> m_socket;
@@ -1911,8 +1895,8 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::ExternalFS
return m_value_list;
}
- ExternalFSIVelocityBoundaryCondition(const Mesh<Connectivity<Dimension>>& mesh,
- const MeshNodeBoundary<Dimension>& mesh_node_boundary,
+ ExternalFSIVelocityBoundaryCondition(const Mesh<Dimension>& mesh,
+ const MeshNodeBoundary& mesh_node_boundary,
const std::shared_ptr<const Socket>& socket)
: m_node_to_cell_matrix{mesh.connectivity().nodeToCellMatrix()},
m_mesh_node_boundary{mesh_node_boundary},
@@ -1927,11 +1911,11 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::ExternalFS
~ExternalFSIVelocityBoundaryCondition() = default;
};
-template <size_t Dimension>
-class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::VelocityBoundaryCondition
+template <MeshConcept MeshType>
+class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<MeshType>::VelocityBoundaryCondition
{
private:
- const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
+ const MeshNodeBoundary m_mesh_node_boundary;
const Array<const TinyVector<Dimension>> m_value_list;
@@ -1948,7 +1932,7 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::VelocityBo
return m_value_list;
}
- VelocityBoundaryCondition(const MeshNodeBoundary<Dimension>& mesh_node_boundary,
+ VelocityBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary,
const Array<const TinyVector<Dimension>>& value_list)
: m_mesh_node_boundary{mesh_node_boundary}, m_value_list{value_list}
{}
@@ -1956,14 +1940,14 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::VelocityBo
~VelocityBoundaryCondition() = default;
};
-template <size_t Dimension>
-class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::SymmetryBoundaryCondition
+template <MeshConcept MeshType>
+class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<MeshType>::SymmetryBoundaryCondition
{
public:
using Rd = TinyVector<Dimension, double>;
private:
- const MeshFlatNodeBoundary<Dimension> m_mesh_flat_node_boundary;
+ const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_boundary;
public:
const Rd&
@@ -1984,7 +1968,7 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::SymmetryBo
return m_mesh_flat_node_boundary.nodeList();
}
- SymmetryBoundaryCondition(const MeshFlatNodeBoundary<Dimension>& mesh_flat_node_boundary)
+ SymmetryBoundaryCondition(const MeshFlatNodeBoundary<MeshType>& mesh_flat_node_boundary)
: m_mesh_flat_node_boundary(mesh_flat_node_boundary)
{
;
@@ -1993,11 +1977,11 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::SymmetryBo
~SymmetryBoundaryCondition() = default;
};
-template <size_t Dimension>
-class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::CouplingBoundaryCondition
+template <MeshConcept MeshType>
+class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<MeshType>::CouplingBoundaryCondition
{
- private:
- const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
+ private:
+ const MeshNodeBoundary m_mesh_node_boundary;
const size_t m_label;
public:
@@ -2008,21 +1992,22 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::CouplingBo
}
size_t
- label() const {
+ label() const
+ {
return m_label;
}
-
- CouplingBoundaryCondition(const MeshNodeBoundary<Dimension>& mesh_node_boundary,const size_t label)
+
+ CouplingBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const size_t label)
: m_mesh_node_boundary{mesh_node_boundary}, m_label{label}
{
;
}
~CouplingBoundaryCondition() = default;
-
-};
+};
-LocalDtAcousticSolverHandler::LocalDtAcousticSolverHandler(const std::shared_ptr<const IMesh>& i_mesh1, const std::shared_ptr<const IMesh>& i_mesh2)
+LocalDtAcousticSolverHandler::LocalDtAcousticSolverHandler(const std::shared_ptr<const MeshVariant>& i_mesh1,
+ const std::shared_ptr<const MeshVariant>& i_mesh2)
{
if (not i_mesh1) {
throw NormalError("mesh1 not defined");
@@ -2031,27 +2016,26 @@ LocalDtAcousticSolverHandler::LocalDtAcousticSolverHandler(const std::shared_ptr
if (not i_mesh2) {
throw NormalError("mesh2 not defined");
}
-
- if (not (i_mesh1->dimension()==i_mesh2->dimension())) {
- throw NormalError("mesh1 and mesh2 must have the same dimension");
- }
-
-
- switch (i_mesh1->dimension()) {
- case 1: {
- m_acoustic_solver = std::make_unique<LocalDtAcousticSolver<1>>();
- break;
- }
- case 2: {
- m_acoustic_solver = std::make_unique<LocalDtAcousticSolver<2>>();
- break;
- }
- case 3: {
- m_acoustic_solver = std::make_unique<LocalDtAcousticSolver<3>>();
- break;
- }
- default: {
- throw UnexpectedError("invalid mesh dimension");
- }
- }
+ std::visit(
+ [&](auto&& first_mesh, auto&& second_mesh) {
+ using FirstMeshType = mesh_type_t<decltype(first_mesh)>;
+ using SecondMeshType = mesh_type_t<decltype(second_mesh)>;
+
+ if constexpr (!std::is_same_v<FirstMeshType,
+ SecondMeshType>) { // <- ICI POUR LE TEST SUR LES TYPES DE MAILLAGES
+ throw NormalError("incompatible mesh types");
+ }
+ },
+ i_mesh1->variant(), i_mesh2->variant());
+
+ std::visit(
+ [&](auto&& mesh) {
+ using MeshType = mesh_type_t<decltype(mesh)>;
+ if constexpr (is_polygonal_mesh_v<MeshType>) {
+ m_acoustic_solver = std::make_unique<LocalDtAcousticSolver<MeshType>>();
+ } else {
+ throw NormalError("unexpected mesh type");
+ }
+ },
+ i_mesh1->variant());
}
diff --git a/src/scheme/LocalDtAcousticSolver.hpp b/src/scheme/LocalDtAcousticSolver.hpp
index c8d23daedde06a12c26a4a1086d26060c4cd15bb..9abb72d75e3c9ae21ac229b5dc7987e1da316b4c 100644
--- a/src/scheme/LocalDtAcousticSolver.hpp
+++ b/src/scheme/LocalDtAcousticSolver.hpp
@@ -1,13 +1,15 @@
#ifndef LOCAL_DT_ACOUSTIC_SOLVER_HPP
#define LOCAL_DT_ACOUSTIC_SOLVER_HPP
+#include <mesh/MeshTraits.hpp>
+
#include <memory>
#include <tuple>
#include <vector>
class DiscreteFunctionVariant;
class IBoundaryConditionDescriptor;
-class IMesh;
+class MeshVariant;
class ItemValueVariant;
class SubItemValuePerItemVariant;
@@ -33,30 +35,30 @@ class LocalDtAcousticSolverHandler
const std::shared_ptr<const DiscreteFunctionVariant>& p,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const = 0;
- virtual std::tuple<std::shared_ptr<const IMesh>,
+ virtual 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 MeshVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const IMesh>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>>
+ std::shared_ptr<const DiscreteFunctionVariant>>
applyOrder2(const SolverType& solver_type,
- const double& dt1,
- const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
- const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
- const std::shared_ptr<const DiscreteFunctionVariant>& u1,
- const std::shared_ptr<const DiscreteFunctionVariant>& u2,
- const std::shared_ptr<const DiscreteFunctionVariant>& E1,
- const std::shared_ptr<const DiscreteFunctionVariant>& E2,
- const std::shared_ptr<const DiscreteFunctionVariant>& c1,
- const std::shared_ptr<const DiscreteFunctionVariant>& c2,
- const std::shared_ptr<const DiscreteFunctionVariant>& p1,
- const std::shared_ptr<const DiscreteFunctionVariant>& p2,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2) const = 0;
-
- virtual std::tuple<std::shared_ptr<const IMesh>,
+ const double& dt1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p2,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2) const = 0;
+
+ virtual std::tuple<std::shared_ptr<const MeshVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
@@ -67,61 +69,61 @@ class LocalDtAcousticSolverHandler
const std::shared_ptr<const ItemValueVariant>& ur,
const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr) const = 0;
- virtual std::tuple<std::shared_ptr<const IMesh>,
+ virtual 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 MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+ apply(const SolverType& solver_type,
+ const double& dt1,
+ const size_t& q,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p2,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2) const = 0;
+
+ virtual 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 IMesh>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>>
+ std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
apply(const SolverType& solver_type,
const double& dt1,
- const size_t& q,
- const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
- const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
- const std::shared_ptr<const DiscreteFunctionVariant>& u1,
- const std::shared_ptr<const DiscreteFunctionVariant>& u2,
- const std::shared_ptr<const DiscreteFunctionVariant>& E1,
- const std::shared_ptr<const DiscreteFunctionVariant>& E2,
- const std::shared_ptr<const DiscreteFunctionVariant>& c1,
- const std::shared_ptr<const DiscreteFunctionVariant>& c2,
- const std::shared_ptr<const DiscreteFunctionVariant>& p1,
- const std::shared_ptr<const DiscreteFunctionVariant>& p2,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2) const = 0;
-
- virtual std::tuple<std::shared_ptr<const IMesh>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const IMesh>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>>
- apply(const SolverType& solver_type,
- const double& dt1,
- const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
- const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
- const std::shared_ptr<const DiscreteFunctionVariant>& u1,
- const std::shared_ptr<const DiscreteFunctionVariant>& u2,
- const std::shared_ptr<const DiscreteFunctionVariant>& E1,
- const std::shared_ptr<const DiscreteFunctionVariant>& E2,
- const std::shared_ptr<const DiscreteFunctionVariant>& c1,
- const std::shared_ptr<const DiscreteFunctionVariant>& c2,
- const std::shared_ptr<const DiscreteFunctionVariant>& p1,
- const std::shared_ptr<const DiscreteFunctionVariant>& p2,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2) const = 0;
-
- ILocalDtAcousticSolver() = default;
- ILocalDtAcousticSolver(ILocalDtAcousticSolver&&) = default;
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p2,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2) const = 0;
+
+ ILocalDtAcousticSolver() = default;
+ ILocalDtAcousticSolver(ILocalDtAcousticSolver&&) = default;
ILocalDtAcousticSolver& operator=(ILocalDtAcousticSolver&&) = default;
virtual ~ILocalDtAcousticSolver() = default;
};
- template <size_t Dimension>
+ template <MeshConcept MeshType>
class LocalDtAcousticSolver;
std::unique_ptr<ILocalDtAcousticSolver> m_acoustic_solver;
@@ -133,7 +135,8 @@ class LocalDtAcousticSolverHandler
return *m_acoustic_solver;
}
- LocalDtAcousticSolverHandler(const std::shared_ptr<const IMesh>& mesh1,const std::shared_ptr<const IMesh>& mesh2);
+ LocalDtAcousticSolverHandler(const std::shared_ptr<const MeshVariant>& mesh1,
+ const std::shared_ptr<const MeshVariant>& mesh2);
};
#endif // LOCAL_DT_ACOUSTIC_SOLVER_HPP