diff --git a/src/scheme/LocalDtHyperelasticSolver.cpp b/src/scheme/LocalDtHyperelasticSolver.cpp
index 69f3ed208be443d8763475af58284187ad6d048c..1b31007949c9d9b4561b37ede9318d57a9cd91a6 100644
--- a/src/scheme/LocalDtHyperelasticSolver.cpp
+++ b/src/scheme/LocalDtHyperelasticSolver.cpp
@@ -7,6 +7,7 @@
#include <mesh/MeshFlatNodeBoundary.hpp>
#include <mesh/MeshNodeBoundary.hpp>
#include <mesh/SubItemValuePerItemVariant.hpp>
+#include <scheme/CouplingBoundaryConditionDescriptor.hpp>
#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
#include <scheme/DiscreteFunctionP0.hpp>
#include <scheme/DiscreteFunctionUtils.hpp>
@@ -16,14 +17,14 @@
#include <scheme/IBoundaryConditionDescriptor.hpp>
#include <scheme/IDiscreteFunctionDescriptor.hpp>
#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
-#include <scheme/CouplingBoundaryConditionDescriptor.hpp>
#include <utils/Socket.hpp>
#include <variant>
#include <vector>
template <size_t Dimension>
-class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final : public LocalDtHyperelasticSolverHandler::ILocalDtHyperelasticSolver
+class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final
+ : public LocalDtHyperelasticSolverHandler::ILocalDtHyperelasticSolver
{
private:
using Rdxd = TinyMatrix<Dimension>;
@@ -48,7 +49,7 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final : public
NormalStressBoundaryCondition,
SymmetryBoundaryCondition,
VelocityBoundaryCondition,
- CouplingBoundaryCondition>;
+ CouplingBoundaryCondition>;
using BoundaryConditionList = std::vector<BoundaryCondition>;
@@ -293,10 +294,11 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final : public
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()));
- break;
+ 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: {
is_valid_boundary_condition = false;
@@ -317,17 +319,17 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final : public
void _applySymmetryBC(const BoundaryConditionList& bc_list, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br) const;
void _applyVelocityBC(const BoundaryConditionList& bc_list, 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
_applyBoundaryConditions(const BoundaryConditionList& bc_list,
const MeshType& mesh,
@@ -376,58 +378,61 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final : public
return F;
}
- std::tuple<std::vector<NodeId>,
- std::vector<NodeId>>
+ 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::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
{
const MeshType& mesh1 = dynamic_cast<const MeshType&>(*i_mesh1);
const MeshType& mesh2 = dynamic_cast<const MeshType&>(*i_mesh2);
-
+
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;
+ 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);
}
- return std::make_tuple(map1,map2);
+ return std::make_tuple(map1, map2);
}
public:
@@ -475,11 +480,11 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final : public
}
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>& aL1_v,
@@ -493,8 +498,8 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final : public
const std::shared_ptr<const DiscreteFunctionVariant>& u2_v,
const std::shared_ptr<const DiscreteFunctionVariant>& sigma2_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, aL1_v, aT1_v, u1_v, sigma1_v});
std::shared_ptr i_mesh2 = getCommonMesh({rho2_v, aL2_v, aT2_v, u2_v, sigma2_v});
@@ -505,11 +510,11 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final : public
if (not checkDiscretizationType({rho1_v, aL1_v, u1_v, sigma1_v}, DiscreteFunctionType::P0)) {
throw NormalError("hyperelastic solver expects P0 functions");
}
- if (not i_mesh2) {
+ if (not i_mesh2) {
throw NormalError("discrete functions are not defined on the same mesh");
}
- if (not checkDiscretizationType({rho2_v, aL2_v, u2_v, sigma2_v}, DiscreteFunctionType::P0)) {
+ if (not checkDiscretizationType({rho2_v, aL2_v, u2_v, sigma2_v}, DiscreteFunctionType::P0)) {
throw NormalError("acoustic solver expects P0 functions");
}
@@ -527,43 +532,39 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final : public
const DiscreteScalarFunction& aT2 = aT2_v->get<DiscreteScalarFunction>();
const DiscreteTensorFunction& sigma2 = sigma2_v->get<DiscreteTensorFunction>();
-
NodeValuePerCell<const Rdxd> Ajr1 = this->_computeAjr(solver_type, mesh1, rho1 * aL1, rho1 * aT1);
-
NodeValuePerCell<const Rdxd> Ajr2 = this->_computeAjr(solver_type, mesh2, rho2 * aL2, rho2 * aT2);
NodeValue<Rdxd> Ar1 = this->_computeAr(mesh1, Ajr1);
NodeValue<Rd> br1 = this->_computeBr(mesh1, Ajr1, u1, sigma1);
NodeValue<Rdxd> Ar2 = this->_computeAr(mesh2, Ajr2);
NodeValue<Rd> br2 = this->_computeBr(mesh2, Ajr2, u2, sigma2);
+ this->_applyCouplingBC(Ar1, Ar2, br1, br2, map1, map2);
const BoundaryConditionList bc_list1 = this->_getBCList(mesh1, bc_descriptor_list1);
const BoundaryConditionList bc_list2 = this->_getBCList(mesh2, bc_descriptor_list2);
this->_applyBoundaryConditions(bc_list1, mesh1, Ar1, br1);
this->_applyBoundaryConditions(bc_list2, mesh2, Ar2, br2);
- 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, sigma1);
- NodeValue<Rd> ur2 = this->_computeUr(mesh2, Ar2, br2);
- NodeValuePerCell<Rd> Fjr2 = this->_computeFjr(mesh2, Ajr2, ur2, u2, sigma2);
- NodeValue<Rd> CR_ur = this->_computeUr(mesh2, Ar2, br2);
- NodeValuePerCell<Rd> CR_Fjr = this->_computeFjr(mesh2, Ajr2, ur2, u2, sigma2);
+ NodeValue<Rd> ur1 = this->_computeUr(mesh1, Ar1, br1);
+ NodeValuePerCell<Rd> Fjr1 = this->_computeFjr(mesh1, Ajr1, ur1, u1, sigma1);
+ NodeValue<Rd> ur2 = this->_computeUr(mesh2, Ar2, br2);
+ NodeValuePerCell<Rd> Fjr2 = this->_computeFjr(mesh2, Ajr2, ur2, u2, sigma2);
+ NodeValue<Rd> CR_ur = this->_computeUr(mesh2, Ar2, br2);
+ NodeValuePerCell<Rd> CR_Fjr = this->_computeFjr(mesh2, Ajr2, ur2, u2, sigma2);
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>>
+ std::tuple<const std::shared_ptr<const ItemValueVariant>, const std::shared_ptr<const SubItemValuePerItemVariant>>
compute_fluxes(const SolverType& solver_type,
const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
const std::shared_ptr<const DiscreteFunctionVariant>& aL_v,
@@ -571,9 +572,9 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final : public
const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
const std::shared_ptr<const DiscreteFunctionVariant>& sigma_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, aL_v, aT_v, u_v, sigma_v});
if (not i_mesh) {
@@ -605,7 +606,7 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final : public
NodeValue<Rd> ur = this->_computeUr(mesh, Ar, br);
NodeValuePerCell<Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u, sigma);
- 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));
@@ -713,27 +714,27 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final : public
}
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 DiscreteFunctionVariant>,
+ 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,
- const DiscreteTensorFunction& CG,
- std::vector<NodeId>& map1,
- std::vector<NodeId>& map2) const
- {
+ const MeshType& mesh2,
+ const DiscreteScalarFunction& rho,
+ const DiscreteVectorFunction& u,
+ const DiscreteScalarFunction& E,
+ const DiscreteTensorFunction& CG,
+ 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];
}
}
@@ -744,35 +745,29 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final : public
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 {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::make_shared<DiscreteFunctionVariant>(DiscreteTensorFunction(new_mesh2, new_CG))};
}
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 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,
- const std::shared_ptr<const DiscreteFunctionVariant>& rho,
- const std::shared_ptr<const DiscreteFunctionVariant>& u,
- const std::shared_ptr<const DiscreteFunctionVariant>& E,
- const std::shared_ptr<const DiscreteFunctionVariant>& CG,
- std::vector<NodeId>& map1,
- std::vector<NodeId>& map2) const
+ 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,
+ const std::shared_ptr<const DiscreteFunctionVariant>& CG,
+ 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>(),
- CG->get<DiscreteTensorFunction>(),
- map1,
- map2);
+ 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>(), CG->get<DiscreteTensorFunction>(), map1, map2);
}
std::tuple<std::shared_ptr<const IMesh>,
@@ -780,32 +775,32 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final : public
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>>
+ 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>>
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
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2,
+ 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});
@@ -815,63 +810,73 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final : public
std::shared_ptr<const DiscreteFunctionVariant> new_E2 = E2;
std::shared_ptr<const DiscreteFunctionVariant> new_CG2 = CG2;
- 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, aL1, aT1, u1, sigma1, bc_descriptor_list1, rho2, aL2, aT2, u2, sigma2, bc_descriptor_list2, map1, map2);
- const auto [new_m1,new_rho1,new_u1,new_E1,new_CG1] = apply_fluxes(dt1, rho1, u1, E1, CG1, ur1, Fjr1);
+ auto [map1, map2] = _computeMapping(m1, new_m2, bc_descriptor_list1, bc_descriptor_list2);
- double dt2 = dt1/q;
- const double gamma = 1.4;
+ auto [ur1, Fjr1, ur2, Fjr2, CR_ur, CR_Fjr] =
+ compute_fluxes(solver_type, rho1, aL1, aT1, u1, sigma1, bc_descriptor_list1, rho2, aL2, aT2, u2, sigma2,
+ bc_descriptor_list2, map1, map2);
+ const auto [new_m1, new_rho1, new_u1, new_E1, new_CG1] = apply_fluxes(dt1, rho1, u1, E1, CG1, ur1, Fjr1);
+
+ double dt2 = dt1 / q;
+ const double gamma = 1.4;
const TinyMatrix<Dimension> I = identity;
- double sum_dt = 0;
-
+ double sum_dt = 0;
+
size_t fluid = 0;
- if((mu == 0) and (lambda == 0)){
+ if ((mu == 0) and (lambda == 0)) {
fluid = 1;
}
- 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 DiscreteTensorFunction& CG_d = new_CG2->get<DiscreteTensorFunction>();
- 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);
- const DiscreteScalarFunction& aT_d = sqrt(mu/rho_d);
+ 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 DiscreteTensorFunction& CG_d = new_CG2->get<DiscreteTensorFunction>();
+ 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);
- const std::shared_ptr<const DiscreteFunctionVariant>& new_sigma2 = std::make_shared<const DiscreteFunctionVariant>(sigma_d);
- const std::shared_ptr<const DiscreteFunctionVariant>& new_aL2 = std::make_shared<const DiscreteFunctionVariant>(aL_d);
- const std::shared_ptr<const DiscreteFunctionVariant>& new_aT2 = std::make_shared<const DiscreteFunctionVariant>(aT_d);
+ const DiscreteScalarFunction& aT_d = sqrt(mu / rho_d);
- auto [ur2,Fjr2] = compute_fluxes(solver_type,new_rho2,new_aL2,new_aT2,new_u2,new_sigma2,bc_descriptor_list2,CR_ur,CR_Fjr,map2);
-
- std::tie(new_m2,new_rho2,new_u2,new_E2,new_CG2) = apply_fluxes(dt2,new_rho2,new_u2,new_E2,new_CG2,ur2,Fjr2);
+ const std::shared_ptr<const DiscreteFunctionVariant>& new_sigma2 =
+ std::make_shared<const DiscreteFunctionVariant>(sigma_d);
+ const std::shared_ptr<const DiscreteFunctionVariant>& new_aL2 =
+ std::make_shared<const DiscreteFunctionVariant>(aL_d);
+ const std::shared_ptr<const DiscreteFunctionVariant>& new_aT2 =
+ std::make_shared<const DiscreteFunctionVariant>(aT_d);
+
+ auto [ur2, Fjr2] = compute_fluxes(solver_type, new_rho2, new_aL2, new_aT2, new_u2, new_sigma2,
+ bc_descriptor_list2, CR_ur, CR_Fjr, map2);
+
+ std::tie(new_m2, new_rho2, new_u2, new_E2, new_CG2) =
+ apply_fluxes(dt2, new_rho2, new_u2, new_E2, new_CG2, ur2, Fjr2);
sum_dt += dt2;
}
- std::tie(new_m2,new_rho2,new_u2,new_E2,new_CG2) = mesh_correction(new_m1,new_m2,new_rho2,new_u2,new_E2,new_CG2,map1,map2);
+ std::tie(new_m2, new_rho2, new_u2, new_E2, new_CG2) =
+ mesh_correction(new_m1, new_m2, new_rho2, new_u2, new_E2, new_CG2, map1, map2);
- return {new_m1,new_rho1,new_u1,new_E1,new_CG1,new_m2,new_rho2,new_u2,new_E2,new_CG2};
+ return {new_m1, new_rho1, new_u1, new_E1, new_CG1, new_m2, new_rho2, new_u2, new_E2, new_CG2};
}
- LocalDtHyperelasticSolver() = default;
+ LocalDtHyperelasticSolver() = default;
LocalDtHyperelasticSolver(LocalDtHyperelasticSolver&&) = default;
- ~LocalDtHyperelasticSolver() = default;
+ ~LocalDtHyperelasticSolver() = default;
};
template <size_t Dimension>
void
- LocalDtHyperelasticSolverHandler:: LocalDtHyperelasticSolver<Dimension>::_applyPressureBC(const BoundaryConditionList& bc_list,
- const MeshType& mesh,
- NodeValue<Rd>& br) const
+LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::_applyPressureBC(
+ const BoundaryConditionList& bc_list,
+ const MeshType& mesh,
+ NodeValue<Rd>& br) const
{
for (const auto& boundary_condition : bc_list) {
std::visit(
@@ -934,9 +939,10 @@ void
template <size_t Dimension>
void
- LocalDtHyperelasticSolverHandler:: LocalDtHyperelasticSolver<Dimension>::_applyNormalStressBC(const BoundaryConditionList& bc_list,
- const MeshType& mesh,
- NodeValue<Rd>& br) const
+LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::_applyNormalStressBC(
+ const BoundaryConditionList& bc_list,
+ const MeshType& mesh,
+ NodeValue<Rd>& br) const
{
for (const auto& boundary_condition : bc_list) {
std::visit(
@@ -999,9 +1005,10 @@ void
template <size_t Dimension>
void
- LocalDtHyperelasticSolverHandler:: LocalDtHyperelasticSolver<Dimension>::_applySymmetryBC(const BoundaryConditionList& bc_list,
- NodeValue<Rdxd>& Ar,
- NodeValue<Rd>& br) const
+LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::_applySymmetryBC(
+ const BoundaryConditionList& bc_list,
+ NodeValue<Rdxd>& Ar,
+ NodeValue<Rd>& br) const
{
for (const auto& boundary_condition : bc_list) {
std::visit(
@@ -1030,9 +1037,10 @@ void
template <size_t Dimension>
void
- LocalDtHyperelasticSolverHandler:: LocalDtHyperelasticSolver<Dimension>::_applyVelocityBC(const BoundaryConditionList& bc_list,
- NodeValue<Rdxd>& Ar,
- NodeValue<Rd>& br) const
+LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::_applyVelocityBC(
+ const BoundaryConditionList& bc_list,
+ NodeValue<Rdxd>& Ar,
+ NodeValue<Rd>& br) const
{
for (const auto& boundary_condition : bc_list) {
std::visit(
@@ -1066,17 +1074,18 @@ void
}
template <size_t Dimension>
-void LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<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
+void
+LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<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
{
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];
@@ -1085,33 +1094,33 @@ void LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::_ap
br1[node_id1] += br2[node_id2];
br2[node_id2] = br1[node_id1];
-
}
}
template <size_t Dimension>
-void LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<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
+void
+LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<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
{
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];
}
@@ -1138,7 +1147,7 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::Fi
};
template <size_t Dimension>
-class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::PressureBoundaryCondition
+class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::PressureBoundaryCondition
{
private:
const MeshFaceBoundary<Dimension> m_mesh_face_boundary;
@@ -1166,7 +1175,7 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::P
};
template <>
-class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<1>::PressureBoundaryCondition
+class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<1>::PressureBoundaryCondition
{
private:
const MeshNodeBoundary<1> m_mesh_node_boundary;
@@ -1193,7 +1202,7 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<1>::PressureB
};
template <size_t Dimension>
-class LocalDtHyperelasticSolverHandler:: LocalDtHyperelasticSolver<Dimension>::NormalStressBoundaryCondition
+class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::NormalStressBoundaryCondition
{
private:
const MeshFaceBoundary<Dimension> m_mesh_face_boundary;
@@ -1221,7 +1230,7 @@ class LocalDtHyperelasticSolverHandler:: LocalDtHyperelasticSolver<Dimension>::
};
template <>
-class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<1>::NormalStressBoundaryCondition
+class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<1>::NormalStressBoundaryCondition
{
private:
const MeshNodeBoundary<1> m_mesh_node_boundary;
@@ -1248,7 +1257,7 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<1>::NormalStr
};
template <size_t Dimension>
-class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::VelocityBoundaryCondition
+class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::VelocityBoundaryCondition
{
private:
const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
@@ -1277,7 +1286,7 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::V
};
template <size_t Dimension>
-class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::SymmetryBoundaryCondition
+class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::SymmetryBoundaryCondition
{
public:
using Rd = TinyVector<Dimension, double>;
@@ -1316,7 +1325,7 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::S
template <size_t Dimension>
class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::CouplingBoundaryCondition
{
- private:
+ private:
const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
const size_t m_label;
@@ -1328,21 +1337,22 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<Dimension>::Co
}
size_t
- label() const {
+ label() const
+ {
return m_label;
}
-
- CouplingBoundaryCondition(const MeshNodeBoundary<Dimension>& mesh_node_boundary,const size_t label)
+
+ CouplingBoundaryCondition(const MeshNodeBoundary<Dimension>& mesh_node_boundary, const size_t label)
: m_mesh_node_boundary{mesh_node_boundary}, m_label{label}
{
;
}
~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 IMesh>& i_mesh1,
+ const std::shared_ptr<const IMesh>& i_mesh2)
{
if (not i_mesh1) {
throw NormalError("mesh1 not defined");
@@ -1352,7 +1362,7 @@ LocalDtHyperelasticSolverHandler:: LocalDtHyperelasticSolverHandler(const std::s
throw NormalError("mesh2 not defined");
}
- if (not (i_mesh1->dimension()==i_mesh2->dimension())) {
+ if (not(i_mesh1->dimension() == i_mesh2->dimension())) {
throw NormalError("mesh1 and mesh2 must have the same dimension");
}