diff --git a/src/scheme/LocalDtHyperelasticSolver.cpp b/src/scheme/LocalDtHyperelasticSolver.cpp
index 1b31007949c9d9b4561b37ede9318d57a9cd91a6..f66ea554d4a44032352f5b792844429d2f830a09 100644
--- a/src/scheme/LocalDtHyperelasticSolver.cpp
+++ b/src/scheme/LocalDtHyperelasticSolver.cpp
@@ -6,6 +6,8 @@
#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/CouplingBoundaryConditionDescriptor.hpp>
#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
@@ -22,20 +24,21 @@
#include <variant>
#include <vector>
-template <size_t Dimension>
+template <MeshConcept MeshType>
class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final
: public LocalDtHyperelasticSolverHandler::ILocalDtHyperelasticSolver
{
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 DiscreteTensorFunction = DiscreteFunctionP0<Dimension, const Rdxd>;
+ using DiscreteScalarFunction = DiscreteFunctionP0<const double>;
+ using DiscreteVectorFunction = DiscreteFunctionP0<const Rd>;
+ using DiscreteTensorFunction = DiscreteFunctionP0<const Rdxd>;
class FixedBoundaryCondition;
class PressureBoundaryCondition;
@@ -178,7 +181,7 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final
}
NodeValue<Rd>
- _computeBr(const Mesh<Connectivity<Dimension>>& mesh,
+ _computeBr(const Mesh<Dimension>& mesh,
const NodeValuePerCell<const Rdxd>& Ajr,
const DiscreteVectorFunction& u,
const DiscreteTensorFunction& sigma) const
@@ -233,8 +236,7 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver 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(),
@@ -246,8 +248,7 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver 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(),
@@ -255,8 +256,7 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final
bc_list.emplace_back(PressureBoundaryCondition{mesh_node_boundary, node_values});
} 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 =
@@ -269,8 +269,7 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final
const FunctionSymbolId normal_stress_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 Rdxd> node_values =
InterpolateItemValue<Rdxd(Rd)>::template interpolate<ItemType::node>(normal_stress_id, mesh.xr(),
@@ -278,8 +277,7 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final
bc_list.emplace_back(NormalStressBoundaryCondition{mesh_node_boundary, node_values});
} 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 Rdxd> face_values =
@@ -379,13 +377,13 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final
}
std::tuple<std::vector<NodeId>, std::vector<NodeId>>
- _computeMapping(std::shared_ptr<const IMesh>& i_mesh1,
- std::shared_ptr<const IMesh>& i_mesh2,
+ _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);
@@ -414,9 +412,7 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final
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]);
- }
+ 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);
@@ -454,7 +450,7 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final
throw NormalError("hyperelastic 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 DiscreteVectorFunction& u = u_v->get<DiscreteVectorFunction>();
const DiscreteScalarFunction& aL = aL_v->get<DiscreteScalarFunction>();
@@ -518,14 +514,14 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver 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 DiscreteVectorFunction& u1 = u1_v->get<DiscreteVectorFunction>();
const DiscreteScalarFunction& aL1 = aL1_v->get<DiscreteScalarFunction>();
const DiscreteScalarFunction& aT1 = aT1_v->get<DiscreteScalarFunction>();
const DiscreteTensorFunction& sigma1 = sigma1_v->get<DiscreteTensorFunction>();
- const MeshType& mesh2 = dynamic_cast<const MeshType&>(*i_mesh2);
+ const MeshType& mesh2 = *i_mesh2->get<MeshType>();
const DiscreteScalarFunction& rho2 = rho2_v->get<DiscreteScalarFunction>();
const DiscreteVectorFunction& u2 = u2_v->get<DiscreteVectorFunction>();
const DiscreteScalarFunction& aL2 = aL2_v->get<DiscreteScalarFunction>();
@@ -585,7 +581,7 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final
throw NormalError("hyperelastic 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 DiscreteVectorFunction& u = u_v->get<DiscreteVectorFunction>();
const DiscreteScalarFunction& aL = aL_v->get<DiscreteScalarFunction>();
@@ -612,7 +608,7 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final
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>,
@@ -675,13 +671,14 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver 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::make_shared<DiscreteFunctionVariant>(DiscreteTensorFunction(new_mesh, new_CG))};
}
- 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>,
@@ -703,17 +700,17 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final
throw NormalError("hyperelastic solver expects P0 functions");
}
- return this->apply_fluxes(dt, //
- dynamic_cast<const MeshType&>(*i_mesh), //
- rho->get<DiscreteScalarFunction>(), //
- u->get<DiscreteVectorFunction>(), //
- E->get<DiscreteScalarFunction>(), //
- CG->get<DiscreteTensorFunction>(), //
- ur->get<NodeValue<const Rd>>(), //
+ return this->apply_fluxes(dt, //
+ *i_mesh->get<MeshType>(), //
+ rho->get<DiscreteScalarFunction>(), //
+ u->get<DiscreteVectorFunction>(), //
+ E->get<DiscreteScalarFunction>(), //
+ CG->get<DiscreteTensorFunction>(), //
+ ur->get<NodeValue<const Rd>>(), //
Fjr->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>,
@@ -745,19 +742,20 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final
CellValue<double> new_E = copy(E.cellValues());
CellValue<Rdxd> new_CG = copy(CG.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::make_shared<DiscreteFunctionVariant>(DiscreteTensorFunction(new_mesh2, new_CG))};
}
- 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>>
- mesh_correction(const std::shared_ptr<const IMesh>& i_mesh1,
- const std::shared_ptr<const IMesh>& i_mesh2,
+ 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,
@@ -765,17 +763,17 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final
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),
+ return this->mesh_correction(*i_mesh1->get<MeshType>(), *i_mesh2->get<MeshType>(),
rho->get<DiscreteScalarFunction>(), u->get<DiscreteVectorFunction>(),
E->get<DiscreteScalarFunction>(), CG->get<DiscreteTensorFunction>(), 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 DiscreteFunctionVariant>,
- std::shared_ptr<const IMesh>,
+ std::shared_ptr<const MeshVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
@@ -802,8 +800,8 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final
const double& mu,
const double& lambda) const
{
- std::shared_ptr<const IMesh> new_m2 = getCommonMesh({rho2, aL2, aT2, u2, sigma2});
- std::shared_ptr<const IMesh> m1 = getCommonMesh({rho1, aL1, aT1, u1, sigma1});
+ std::shared_ptr new_m2 = getCommonMesh({rho2, aL2, aT2, u2, sigma2});
+ std::shared_ptr m1 = getCommonMesh({rho1, aL1, aT1, u1, sigma1});
std::shared_ptr<const DiscreteFunctionVariant> new_rho2 = rho2;
std::shared_ptr<const DiscreteFunctionVariant> new_u2 = u2;
@@ -837,7 +835,7 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final
const DiscreteScalarFunction& E_d = new_E2->get<DiscreteScalarFunction>();
const DiscreteScalarFunction& eps = E_d - 0.5 * dot(u_d, u_d);
const DiscreteTensorFunction& CG_d = new_CG2->get<DiscreteTensorFunction>();
- const DiscreteScalarFunction& p = fluid * (gamma - 1) * rho_d * eps;
+ const DiscreteScalarFunction& p = fluid * (gamma - 1) * rho_d * eps;
const DiscreteTensorFunction& sigma_d =
mu / sqrt(det(CG_d)) * (CG_d - I) + lambda / sqrt(det(CG_d)) * log(sqrt(det(CG_d))) * I - p * I;
const DiscreteScalarFunction& aL_d = sqrt((lambda + 2 * mu) / rho_d + gamma * p / rho_d);
@@ -871,9 +869,9 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final
~LocalDtHyperelasticSolver() = default;
};
-template <size_t Dimension>
+template <MeshConcept MeshType>
void
-LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::_applyPressureBC(
+LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<MeshType>::_applyPressureBC(
const BoundaryConditionList& bc_list,
const MeshType& mesh,
NodeValue<Rd>& br) const
@@ -883,7 +881,7 @@ LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::_applyPr
[&](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);
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
if constexpr (Dimension == 1) {
const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
@@ -937,9 +935,9 @@ LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::_applyPr
}
}
-template <size_t Dimension>
+template <MeshConcept MeshType>
void
-LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::_applyNormalStressBC(
+LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<MeshType>::_applyNormalStressBC(
const BoundaryConditionList& bc_list,
const MeshType& mesh,
NodeValue<Rd>& br) const
@@ -949,7 +947,7 @@ LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::_applyNo
[&](auto&& bc) {
using T = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<NormalStressBoundaryCondition, T>) {
- MeshData<Dimension>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
if constexpr (Dimension == 1) {
const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
@@ -1003,9 +1001,9 @@ LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::_applyNo
}
}
-template <size_t Dimension>
+template <MeshConcept MeshType>
void
-LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::_applySymmetryBC(
+LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<MeshType>::_applySymmetryBC(
const BoundaryConditionList& bc_list,
NodeValue<Rdxd>& Ar,
NodeValue<Rd>& br) const
@@ -1035,9 +1033,9 @@ LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::_applySy
}
}
-template <size_t Dimension>
+template <MeshConcept MeshType>
void
-LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::_applyVelocityBC(
+LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<MeshType>::_applyVelocityBC(
const BoundaryConditionList& bc_list,
NodeValue<Rdxd>& Ar,
NodeValue<Rd>& br) const
@@ -1073,9 +1071,9 @@ LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::_applyVe
}
}
-template <size_t Dimension>
+template <MeshConcept MeshType>
void
-LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::_applyCouplingBC(
+LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<MeshType>::_applyCouplingBC(
NodeValue<Rdxd>& Ar1,
NodeValue<Rdxd>& Ar2,
NodeValue<Rd>& br1,
@@ -1097,9 +1095,9 @@ LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::_applyCo
}
}
-template <size_t Dimension>
+template <MeshConcept MeshType>
void
-LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::_applyCouplingBC(
+LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<MeshType>::_applyCouplingBC(
const MeshType& mesh,
NodeValue<Rd>& ur,
NodeValue<Rd>& CR_ur,
@@ -1126,11 +1124,11 @@ LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::_applyCo
}
}
-template <size_t Dimension>
-class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::FixedBoundaryCondition
+template <MeshConcept MeshType>
+class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<MeshType>::FixedBoundaryCondition
{
private:
- const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
+ const MeshNodeBoundary m_mesh_node_boundary;
public:
const Array<const NodeId>&
@@ -1139,18 +1137,16 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::Fi
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 LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::PressureBoundaryCondition
+template <MeshConcept MeshType>
+class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<MeshType>::PressureBoundaryCondition
{
private:
- const MeshFaceBoundary<Dimension> m_mesh_face_boundary;
+ const MeshFaceBoundary m_mesh_face_boundary;
const Array<const double> m_value_list;
public:
@@ -1166,8 +1162,7 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::Pr
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}
{}
@@ -1175,10 +1170,10 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::Pr
};
template <>
-class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<1>::PressureBoundaryCondition
+class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<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:
@@ -1194,18 +1189,18 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<1>::PressureBo
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 LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::NormalStressBoundaryCondition
+template <MeshConcept MeshType>
+class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<MeshType>::NormalStressBoundaryCondition
{
private:
- const MeshFaceBoundary<Dimension> m_mesh_face_boundary;
+ const MeshFaceBoundary m_mesh_face_boundary;
const Array<const Rdxd> m_value_list;
public:
@@ -1221,8 +1216,7 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::No
return m_value_list;
}
- NormalStressBoundaryCondition(const MeshFaceBoundary<Dimension>& mesh_face_boundary,
- const Array<const Rdxd>& value_list)
+ NormalStressBoundaryCondition(const MeshFaceBoundary& mesh_face_boundary, const Array<const Rdxd>& value_list)
: m_mesh_face_boundary{mesh_face_boundary}, m_value_list{value_list}
{}
@@ -1230,10 +1224,10 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::No
};
template <>
-class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<1>::NormalStressBoundaryCondition
+class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Mesh<1>>::NormalStressBoundaryCondition
{
private:
- const MeshNodeBoundary<1> m_mesh_node_boundary;
+ const MeshNodeBoundary m_mesh_node_boundary;
const Array<const Rdxd> m_value_list;
public:
@@ -1249,18 +1243,18 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<1>::NormalStre
return m_value_list;
}
- NormalStressBoundaryCondition(const MeshNodeBoundary<1>& mesh_node_boundary, const Array<const Rdxd>& value_list)
+ NormalStressBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const Array<const Rdxd>& value_list)
: m_mesh_node_boundary{mesh_node_boundary}, m_value_list{value_list}
{}
~NormalStressBoundaryCondition() = default;
};
-template <size_t Dimension>
-class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::VelocityBoundaryCondition
+template <MeshConcept MeshType>
+class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<MeshType>::VelocityBoundaryCondition
{
private:
- const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
+ const MeshNodeBoundary m_mesh_node_boundary;
const Array<const TinyVector<Dimension>> m_value_list;
@@ -1277,7 +1271,7 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::Ve
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}
{}
@@ -1285,14 +1279,14 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::Ve
~VelocityBoundaryCondition() = default;
};
-template <size_t Dimension>
-class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::SymmetryBoundaryCondition
+template <MeshConcept MeshType>
+class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<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&
@@ -1313,7 +1307,7 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::Sy
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)
{
;
@@ -1322,11 +1316,11 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::Sy
~SymmetryBoundaryCondition() = default;
};
-template <size_t Dimension>
-class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::CouplingBoundaryCondition
+template <MeshConcept MeshType>
+class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<MeshType>::CouplingBoundaryCondition
{
private:
- const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
+ const MeshNodeBoundary m_mesh_node_boundary;
const size_t m_label;
public:
@@ -1342,7 +1336,7 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::Co
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}
{
;
@@ -1351,8 +1345,8 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::Co
~CouplingBoundaryCondition() = default;
};
-LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolverHandler(const std::shared_ptr<const IMesh>& i_mesh1,
- const std::shared_ptr<const IMesh>& i_mesh2)
+LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolverHandler(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");
@@ -1361,26 +1355,26 @@ LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolverHandler(const std::sh
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_hyperelastic_solver = std::make_unique<LocalDtHyperelasticSolver<1>>();
- break;
- }
- case 2: {
- m_hyperelastic_solver = std::make_unique<LocalDtHyperelasticSolver<2>>();
- break;
- }
- case 3: {
- m_hyperelastic_solver = std::make_unique<LocalDtHyperelasticSolver<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_hyperelastic_solver = std::make_unique<LocalDtHyperelasticSolver<MeshType>>();
+ } else {
+ throw NormalError("unexpected mesh type");
+ }
+ },
+ i_mesh1->variant());
}
diff --git a/src/scheme/LocalDtHyperelasticSolver.hpp b/src/scheme/LocalDtHyperelasticSolver.hpp
index 8cc83e5d468f913817facc76d210d4e003b44ffc..731e0810a4d11a233df73fc3b2ee24af429e9699 100644
--- a/src/scheme/LocalDtHyperelasticSolver.hpp
+++ b/src/scheme/LocalDtHyperelasticSolver.hpp
@@ -1,12 +1,15 @@
-#ifndef LOCAL_DT_HYPERELASTIC_SOLVER_HPP
-#define LOCAL_DT_HYPERELASTIC_SOLVER_HPP
+#ifndef LOCAL_DT_HYPERELASTIC_SOLVER_HPP
+#define LOCAL_DT_HYPERELASTIC_SOLVER_HPP
+
+#include <mesh/MeshTraits.hpp>
#include <memory>
#include <tuple>
#include <vector>
+class DiscreteFunctionVariant;
class IBoundaryConditionDescriptor;
-class IMesh;
+class MeshVariant;
class ItemValueVariant;
class SubItemValuePerItemVariant;
class DiscreteFunctionVariant;
@@ -34,7 +37,7 @@ class LocalDtHyperelasticSolverHandler
const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
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>,
@@ -47,46 +50,46 @@ class LocalDtHyperelasticSolverHandler
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 DiscreteFunctionVariant>,
- std::shared_ptr<const IMesh>,
+ 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(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>& CG1,
- const std::shared_ptr<const DiscreteFunctionVariant>& CG2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& CG2,
const std::shared_ptr<const DiscreteFunctionVariant>& aL1,
- const std::shared_ptr<const DiscreteFunctionVariant>& aL2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aL2,
const std::shared_ptr<const DiscreteFunctionVariant>& aT1,
- const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
- const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2,
- const double& mu,
- const double& lambda) const = 0;
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2,
+ const double& mu,
+ const double& lambda) const = 0;
- ILocalDtHyperelasticSolver() = default;
+ ILocalDtHyperelasticSolver() = default;
ILocalDtHyperelasticSolver(ILocalDtHyperelasticSolver&&) = default;
ILocalDtHyperelasticSolver& operator=(ILocalDtHyperelasticSolver&&) = default;
virtual ~ILocalDtHyperelasticSolver() = default;
};
- template <size_t Dimension>
+ template <MeshConcept MeshType>
class LocalDtHyperelasticSolver;
std::unique_ptr<ILocalDtHyperelasticSolver> m_hyperelastic_solver;
@@ -98,7 +101,8 @@ class LocalDtHyperelasticSolverHandler
return *m_hyperelastic_solver;
}
- LocalDtHyperelasticSolverHandler(const std::shared_ptr<const IMesh>& mesh1, const std::shared_ptr<const IMesh>& mesh2);
+ LocalDtHyperelasticSolverHandler(const std::shared_ptr<const MeshVariant>& mesh1,
+ const std::shared_ptr<const MeshVariant>& mesh2);
};
#endif // HYPERELASTIC_SOLVER_HPP