Select Git revision
Order2HyperelasticSolver.cpp
Order2HyperelasticSolver.cpp 74.00 KiB
#include <scheme/Order2HyperelasticSolver.hpp>
#include <language/utils/InterpolateItemValue.hpp>
#include <mesh/ItemValueUtils.hpp>
#include <mesh/ItemValueVariant.hpp>
#include <mesh/MeshFaceBoundary.hpp>
#include <mesh/MeshFlatNodeBoundary.hpp>
#include <mesh/MeshNodeBoundary.hpp>
#include <mesh/MeshTraits.hpp>
#include <mesh/StencilArray.hpp>
#include <mesh/StencilManager.hpp>
#include <mesh/SubItemValuePerItemVariant.hpp>
#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
#include <scheme/DiscreteFunctionDPk.hpp>
#include <scheme/DiscreteFunctionDPkVariant.hpp>
#include <scheme/DiscreteFunctionP0.hpp>
#include <scheme/DiscreteFunctionUtils.hpp>
#include <scheme/DiscreteFunctionVariant.hpp>
#include <scheme/ExternalBoundaryConditionDescriptor.hpp>
#include <scheme/FixedBoundaryConditionDescriptor.hpp>
#include <scheme/IBoundaryConditionDescriptor.hpp>
#include <scheme/IDiscreteFunctionDescriptor.hpp>
#include <scheme/PolynomialReconstruction.hpp>
#include <scheme/PolynomialReconstructionDescriptor.hpp>
#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
#include <utils/Socket.hpp>
#include <variant>
#include <vector>
template <MeshConcept MeshType>
class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final
: public Order2HyperelasticSolverHandler::IOrder2HyperelasticSolver
{
private:
static constexpr size_t Dimension = MeshType::Dimension;
using Rdxd = TinyMatrix<Dimension>;
using Rd = TinyVector<Dimension>;
using MeshDataType = MeshData<MeshType>;
using DiscreteScalarFunction = DiscreteFunctionP0<const double>;
using DiscreteVectorFunction = DiscreteFunctionP0<const Rd>;
using DiscreteTensorFunction = DiscreteFunctionP0<const Rdxd>;
class FixedBoundaryCondition;
class PressureBoundaryCondition;
class NormalStressBoundaryCondition;
class SymmetryBoundaryCondition;
class VelocityBoundaryCondition;
using BoundaryCondition = std::variant<FixedBoundaryCondition,
PressureBoundaryCondition,
NormalStressBoundaryCondition,
SymmetryBoundaryCondition,
VelocityBoundaryCondition>;
using BoundaryConditionList = std::vector<BoundaryCondition>;
NodeValuePerCell<const Rdxd>
_computeGlaceAjr(const MeshType& mesh, const DiscreteScalarFunction& rhoaL, const DiscreteScalarFunction& rhoaT) const
{
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
const NodeValuePerCell<const Rd> njr = mesh_data.njr();
NodeValuePerCell<Rdxd> Ajr{mesh.connectivity()};
const Rdxd I = identity; parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
const size_t& nb_nodes = Ajr.numberOfSubValues(j);
const double& rhoaL_j = rhoaL[j];
const double& rhoaT_j = rhoaT[j];
for (size_t r = 0; r < nb_nodes; ++r) {
const Rdxd& M = tensorProduct(Cjr(j, r), njr(j, r));
Ajr(j, r) = rhoaL_j * M + rhoaT_j * (l2Norm(Cjr(j, r)) * I - M);
}
});
return Ajr;
}
NodeValuePerCell<const Rdxd>
_computeEucclhydAjr(const MeshType& mesh,
const DiscreteScalarFunction& rhoaL,
const DiscreteScalarFunction& rhoaT) const
{
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();
const NodeValuePerFace<const Rd> nlr = mesh_data.nlr();
const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
NodeValuePerCell<Rdxd> Ajr{mesh.connectivity()};
parallel_for(
Ajr.numberOfValues(), PUGS_LAMBDA(size_t jr) { Ajr[jr] = zero; });
const Rdxd I = identity;
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
const auto& cell_nodes = cell_to_node_matrix[j];
const auto& cell_faces = cell_to_face_matrix[j];
const double& rho_aL = rhoaL[j];
const double& rho_aT = rhoaT[j];
for (size_t L = 0; L < cell_faces.size(); ++L) {
const FaceId& l = cell_faces[L];
const auto& face_nodes = face_to_node_matrix[l];
auto local_node_number_in_cell = [&](NodeId node_number) {
for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
if (node_number == cell_nodes[i_node]) {
return i_node;
}
}
return std::numeric_limits<size_t>::max();
};
for (size_t rl = 0; rl < face_nodes.size(); ++rl) {
const size_t R = local_node_number_in_cell(face_nodes[rl]);
const Rdxd& M = tensorProduct(Nlr(l, rl), nlr(l, rl));
Ajr(j, R) += rho_aL * M + rho_aT * (l2Norm(Nlr(l, rl)) * I - M);
}
}
});
return Ajr;
}
NodeValuePerCell<const Rdxd>
_computeAjr(const SolverType& solver_type,
const MeshType& mesh,
const DiscreteScalarFunction& rhoaL, const DiscreteScalarFunction& rhoaT) const
{
if constexpr (Dimension == 1) {
return _computeGlaceAjr(mesh, rhoaL, rhoaT);
} else {
switch (solver_type) {
case SolverType::Glace: {
return _computeGlaceAjr(mesh, rhoaL, rhoaT);
}
case SolverType::Eucclhyd: {
return _computeEucclhydAjr(mesh, rhoaL, rhoaT);
}
default: {
throw UnexpectedError("invalid solver type");
}
}
}
}
NodeValue<Rdxd>
_computeAr(const MeshType& mesh, const NodeValuePerCell<const Rdxd>& Ajr) const
{
const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
NodeValue<Rdxd> Ar{mesh.connectivity()};
parallel_for(
mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
Rdxd sum = zero;
const auto& node_to_cell = node_to_cell_matrix[r];
const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(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];
sum += Ajr(J, R);
}
Ar[r] = sum;
});
return Ar;
}
NodeValue<Rd>
_computeBr(const Mesh<Dimension>& mesh,
const NodeValuePerCell<const Rdxd>& Ajr,
DiscreteFunctionDPk<Dimension, const Rd> DPk_u,
DiscreteFunctionDPk<Dimension, const Rdxd> DPk_sigma) const
{
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
const NodeValuePerCell<const Rd>& Cjr = mesh_data.Cjr();
const auto& xr = mesh.xr();
const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
NodeValue<Rd> b{mesh.connectivity()};
parallel_for(
mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
const auto& node_to_cell = node_to_cell_matrix[r];
const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(r);
Rd br = zero;
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];
// br += Ajr(J, R) * DPk_u[J](xr[r]) - DPk_sigma[j](xr[r]) * Cjr(J, R); br += Ajr(J, R) * DPk_u[J](xr[r]) - DPk_sigma[J](xr[r]) * Cjr(J, R);
}
b[r] = br;
});
return b;
}
BoundaryConditionList
_getBCList(const MeshType& mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
{
BoundaryConditionList bc_list;
for (const auto& bc_descriptor : bc_descriptor_list) {
bool is_valid_boundary_condition = true;
switch (bc_descriptor->type()) {
case IBoundaryConditionDescriptor::Type::symmetry: {
bc_list.emplace_back(
SymmetryBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())));
break;
}
case IBoundaryConditionDescriptor::Type::fixed: {
bc_list.emplace_back(FixedBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())));
break;
}
case IBoundaryConditionDescriptor::Type::dirichlet: {
const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
if (dirichlet_bc_descriptor.name() == "velocity") {
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(),
mesh.xr(),
mesh_node_boundary.nodeList());
bc_list.emplace_back(VelocityBoundaryCondition{mesh_node_boundary, value_list});
} else if (dirichlet_bc_descriptor.name() == "pressure") {
const FunctionSymbolId pressure_id = dirichlet_bc_descriptor.rhsSymbolId();
if constexpr (Dimension == 1) {
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});
} else {
MeshFaceBoundary mesh_face_boundary = getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor());
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
Array<const double> face_values =
InterpolateItemValue<double(Rd)>::template interpolate<ItemType::face>(pressure_id, mesh_data.xl(),
mesh_face_boundary.faceList());
bc_list.emplace_back(PressureBoundaryCondition{mesh_face_boundary, face_values});
}
} else if (dirichlet_bc_descriptor.name() == "normal-stress") {
const FunctionSymbolId normal_stress_id = dirichlet_bc_descriptor.rhsSymbolId();
if constexpr (Dimension == 1) {
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(),
mesh_node_boundary.nodeList());
bc_list.emplace_back(NormalStressBoundaryCondition{mesh_node_boundary, node_values});
} else {
MeshFaceBoundary mesh_face_boundary = getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor());
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
Array<const Rdxd> face_values =
InterpolateItemValue<Rdxd(Rd)>::template interpolate<ItemType::face>(normal_stress_id, mesh_data.xl(),
mesh_face_boundary.faceList());
bc_list.emplace_back(NormalStressBoundaryCondition{mesh_face_boundary, face_values});
}
} else {
is_valid_boundary_condition = false;
}
break;
}
default: {
is_valid_boundary_condition = false;
}
}
if (not is_valid_boundary_condition) {
std::ostringstream error_msg;
error_msg << *bc_descriptor << " is an invalid boundary condition for hyperelastic solver";
throw NormalError(error_msg.str());
}
}
return bc_list;
}
void _applyPressureBC(const BoundaryConditionList& bc_list, const MeshType& mesh, NodeValue<Rd>& br) const;
void _applyNormalStressBC(const BoundaryConditionList& bc_list, const MeshType& mesh, NodeValue<Rd>& br) const;
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
_applyBoundaryConditions(const BoundaryConditionList& bc_list,
const MeshType& mesh,
NodeValue<Rdxd>& Ar,
NodeValue<Rd>& br) const
{
this->_applyPressureBC(bc_list, mesh, br);
this->_applyNormalStressBC(bc_list, mesh, br);
this->_applySymmetryBC(bc_list, Ar, br);
this->_applyVelocityBC(bc_list, Ar, br);
}
NodeValue<const Rd>
_computeUr(const MeshType& mesh, const NodeValue<Rdxd>& Ar, const NodeValue<Rd>& br) const
{
NodeValue<Rd> u{mesh.connectivity()};
parallel_for(
mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { u[r] = inverse(Ar[r]) * br[r]; });
return u;
}
NodeValue<Rd>
_computeBr(const Mesh<Dimension>& mesh,
const NodeValuePerCell<const Rdxd>& Ajr,
const DiscreteVectorFunction& u,
const DiscreteTensorFunction& sigma) const
{
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
const NodeValuePerCell<const Rd>& Cjr = mesh_data.Cjr();
const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
NodeValue<Rd> b{mesh.connectivity()};
parallel_for(
mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
const auto& node_to_cell = node_to_cell_matrix[r];
const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(r);
Rd br = zero;
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];
br += Ajr(J, R) * u[J] - sigma[J] * Cjr(J, R);
}
b[r] = br;
});
return b;
}
NodeValuePerCell<Rd>
_computeFjr(const MeshType& mesh,
const NodeValuePerCell<const Rdxd>& Ajr,
const NodeValue<const Rd>& ur,
const DiscreteVectorFunction& u,
const DiscreteTensorFunction& sigma) const
{
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
NodeValuePerCell<Rd> F{mesh.connectivity()};
parallel_for(
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) {
F(j, r) = -Ajr(j, r) * (u[j] - ur[cell_nodes[r]]) + sigma[j] * Cjr(j, r);
}
});
return F;
}
NodeValue<Rd>
_computeBr(const Mesh<Dimension>& mesh,
const NodeValuePerCell<const Rdxd>& Ajr,
DiscreteFunctionDPk<Dimension, const Rd> DPk_u,
DiscreteFunctionDPk<Dimension, const Rdxd> DPk_S,
DiscreteFunctionDPk<Dimension, const double> DPk_p,
CellValue<double> lambda, NodeValuePerCell<double> coef) const
{
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
const NodeValuePerCell<const Rd>& Cjr = mesh_data.Cjr();
const auto& xr = mesh.xr();
const TinyMatrix<Dimension> I = identity;
const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
NodeValue<Rd> b{mesh.connectivity()};
parallel_for(
mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
const auto& node_to_cell = node_to_cell_matrix[r];
const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(r);
Rd br = zero;
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];
const Rdxd sigma_r = sqrt(lambda[J]+(1-lambda[J])*coef(J,R))*DPk_S[J](xr[r]) - DPk_p[J](xr[r])*I;
br += Ajr(J, R) * DPk_u[J](xr[r]) - sigma_r * Cjr(J, R);
}
b[r] = br;
});
return b;
}
NodeValuePerCell<Rd>
_computeFjr(const MeshType& mesh,
const NodeValuePerCell<const Rdxd>& Ajr,
const NodeValue<const Rd>& ur,
DiscreteFunctionDPk<Dimension, const Rd> DPk_uh,
DiscreteFunctionDPk<Dimension, const Rdxd> DPk_S,
DiscreteFunctionDPk<Dimension, const double> DPk_p,
CellValue<double> lambda, NodeValuePerCell<double> coef) const
{
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
const NodeValue<const Rd>& xr = mesh.xr();
const TinyMatrix<Dimension> I = identity;
const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
NodeValuePerCell<Rd> F{mesh.connectivity()};
parallel_for(
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) {
Rdxd sigma_r = sqrt(lambda[j]+(1-lambda[j])*coef(j,r))*DPk_S[j](xr[cell_nodes[r]]) - DPk_p[j](xr[cell_nodes[r]])*I;
F(j,r) = -Ajr(j,r) * (DPk_uh[j](xr[cell_nodes[r]]) - ur[cell_nodes[r]]) + sigma_r*Cjr(j,r);
}
});
return F;
}
void
_scalar_limiter(const MeshType& mesh,
const DiscreteFunctionP0<const double>& f,
DiscreteFunctionDPk<Dimension, double>& DPk_fh) const
{
MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
StencilManager::BoundaryDescriptorList symmetry_boundary_descriptor_list;
StencilDescriptor stencil_descriptor{1, StencilDescriptor::ConnectionType::by_nodes};
auto stencil = StencilManager::instance().getCellToCellStencilArray(mesh.connectivity(), stencil_descriptor, symmetry_boundary_descriptor_list);
const auto xj = mesh_data.xj();
CellValue<double> alph_J2{mesh.connectivity()};
parallel_for(mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){
const double fj = f[cell_id];
double f_min = fj;
double f_max = fj;
const auto cell_stencil = stencil[cell_id];
for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
f_min = std::min(f_min, f[cell_stencil[i_cell]]); f_max = std::max(f_max, f[cell_stencil[i_cell]]);
}
double f_bar_min = fj;
double f_bar_max = fj;
for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
const CellId cell_k_id = cell_stencil[i_cell];
const double f_xk = DPk_fh[cell_id](xj[cell_k_id]);
f_bar_min = std::min(f_bar_min, f_xk);
f_bar_max = std::max(f_bar_max, f_xk);
}
const double eps = 1E-14;
double coef1 = 1;
if (std::abs(f_bar_max - fj) > eps) {
coef1 = (f_max - fj) / ((f_bar_max - fj));
}
double coef2 = 1.;
if (std::abs(f_bar_min - fj) > eps) {
coef2 = (f_min - fj) / ((f_bar_min - fj));
}
const double lambda = std::max(0., std::min(1., std::min(coef1, coef2)));
alph_J2[cell_id] = lambda;
auto coefficients = DPk_fh.coefficients(cell_id);
coefficients[0] = (1 - lambda) * f[cell_id] + lambda * coefficients[0];
for (size_t i = 1; i < coefficients.size(); ++i) {
coefficients[i] *= lambda;
}
});
}
void
_vector_limiter(const MeshType& mesh,
const DiscreteFunctionP0<const Rd>& f,
DiscreteFunctionDPk<Dimension, Rd>& DPK_fh,
const DiscreteFunctionDPk<Dimension, const double>& DPk_ph) const
{
MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
StencilManager::BoundaryDescriptorList symmetry_boundary_descriptor_list;
StencilDescriptor stencil_descriptor{1, StencilDescriptor::ConnectionType::by_nodes};
auto stencil = StencilManager::instance().getCellToCellStencilArray(mesh.connectivity(), stencil_descriptor,
symmetry_boundary_descriptor_list);
const auto xj = mesh_data.xj();
CellValue<Rd> n{mesh.connectivity()};
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
auto coefficients_p = DPk_ph.coefficients(cell_id);
Rd grad_p = zero;
for(size_t i = 1; i < coefficients_p.size(); ++i){
grad_p[i-1] = coefficients_p[i];
}
const double norm_grad_p = l2Norm(grad_p);
if(norm_grad_p == 0){
n[cell_id] = zero;
}
else {
n[cell_id] = (1/norm_grad_p) * grad_p;
}
});
const CellValue<Rd> t{mesh.connectivity()};
const CellValue<Rd> l{mesh.connectivity()};
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
const Rd nj = n[cell_id];
Rd a = zero;
if(l2Norm(nj) != 0){
if((nj[0] / l2Norm(nj) != 1) and (nj[0] / l2Norm(nj) != -1)){
a[0] = 1.;
}
else {
a[1] = 1.;
}
}
Rd tj = a - dot(a,nj) * nj;
const double& norm_tj = l2Norm(tj);
if(norm_tj == 0){
tj = zero;
}
else {
tj = (1/norm_tj) * tj;
}
t[cell_id] = tj;
Rd lj = zero;
if(Dimension == 3){
lj[0] = nj[1]*tj[2] - nj[2]*tj[1];
lj[1] = nj[2]*tj[0] - nj[0]*tj[2];
lj[2] = nj[0]*tj[1] - nj[1]*tj[0];
const double& norm_lj = l2Norm(lj);
if(norm_lj != 0){
lj = (1/norm_lj) * lj;
}
}
l[cell_id] = lj;
});
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){
const double fn = dot(f[cell_id],n[cell_id]);
const double ft = dot(f[cell_id],t[cell_id]);
const double fl = dot(f[cell_id],l[cell_id]);
double fn_min = fn;
double fn_max = fn;
double ft_min = ft;
double ft_max = ft;
double fl_min = fl;
double fl_max = fl;
const auto cell_stencil = stencil[cell_id];
for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
const double fn_k = dot(f[cell_stencil[i_cell]],n[cell_stencil[i_cell]]);
fn_min = std::min(fn_min,fn_k);
fn_max = std::max(fn_max,fn_k);
const double ft_k = dot(f[cell_stencil[i_cell]],t[cell_stencil[i_cell]]);
ft_min = std::min(ft_min,ft_k);
ft_max = std::max(ft_max,ft_k);
const double fl_k = dot(f[cell_stencil[i_cell]],l[cell_stencil[i_cell]]);
fl_min = std::min(fl_min,fl_k);
fl_max = std::max(fl_max,fl_k);
}
double fn_bar_min = fn;
double fn_bar_max = fn;
double ft_bar_min = ft; double ft_bar_max = ft;
double fl_bar_min = fl;
double fl_bar_max = fl;
for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
const CellId cell_k_id = cell_stencil[i_cell];
const Rd f_xk = DPK_fh[cell_id](xj[cell_k_id]);
const double fn_xk = dot(f_xk,n[cell_k_id]);
fn_bar_min = std::min(fn_bar_min,fn_xk);
fn_bar_max = std::max(fn_bar_max,fn_xk);
const double ft_xk = dot(f_xk,t[cell_k_id]);
ft_bar_min = std::min(ft_bar_min,ft_xk);
ft_bar_max = std::max(ft_bar_max,ft_xk);
const double fl_xk = dot(f_xk,l[cell_k_id]);
fl_bar_min = std::min(fl_bar_min,fl_xk);
fl_bar_max = std::max(fl_bar_max,fl_xk);
}
const double eps = 1E-14;
double coef1_n = 1;
if(std::abs(fn_bar_max - fn) > eps){
coef1_n = (fn_max - fn) / (fn_bar_max - fn);
}
double coef2_n = 1;
if(std::abs(fn_bar_min - fn) > eps){
coef2_n = (fn_min - fn) / (fn_bar_min - fn);
}
double coef1_t = 1;
if(std::abs(ft_bar_max - ft) > eps){
coef1_t = (ft_max - ft) / (ft_bar_max - ft);
}
double coef2_t = 1;
if(std::abs(ft_bar_min - ft) > eps){
coef2_t = (ft_min - ft) / (ft_bar_min - ft);
}
double coef1_l = 1;
if(std::abs(fl_bar_max - fl) > eps){
coef1_l = (fl_max - fl) / (fl_bar_max - fl);
}
double coef2_l = 1;
if(std::abs(fl_bar_min - fl) > eps){
coef2_l = (fl_min - fl) / (fl_bar_min - fl);
}
const double lambda_n = std::max(0., std::min(1., std::min(coef1_n,coef2_n)));
const double lambda_t = std::max(0., std::min(1., std::min(coef1_t,coef2_t)));
const double lambda_l = std::max(0., std::min(1., std::min(coef1_l,coef2_l)));
auto coefficients = DPK_fh.coefficients(cell_id);
coefficients[0] = (1 - lambda_n)*fn*n[cell_id] + lambda_n*dot(coefficients[0],n[cell_id])*n[cell_id]
+ (1 - lambda_t)*ft*t[cell_id] + lambda_t*dot(coefficients[0],t[cell_id])*t[cell_id]
+ (1 - lambda_l)*fl*l[cell_id] + lambda_l*dot(coefficients[0],l[cell_id])*l[cell_id];
for(size_t i = 1; i < coefficients.size(); ++i){
coefficients[i] = lambda_n*dot(coefficients[i],n[cell_id])*n[cell_id]
+ lambda_t*dot(coefficients[i],t[cell_id])*t[cell_id]
+ lambda_l*dot(coefficients[i],l[cell_id])*l[cell_id];
}
});
}
CellValue<double>
_matrix_limiter(const MeshType& mesh,
const DiscreteFunctionP0<const Rdxd>& S,
DiscreteFunctionDPk<Dimension, Rdxd>& DPk_Sh) const
{ MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
StencilManager::BoundaryDescriptorList symmetry_boundary_descriptor_list;
StencilDescriptor stencil_descriptor{1, StencilDescriptor::ConnectionType::by_nodes};
auto stencil = StencilManager::instance().getCellToCellStencilArray(mesh.connectivity(), stencil_descriptor,
symmetry_boundary_descriptor_list);
const auto xj = mesh_data.xj();
//A retirer + tard
const auto xr = mesh.xr();
const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
//
CellValue<double> lambda{mesh.connectivity()};
const DiscreteScalarFunction& inv2 = 0.5*trace(transpose(S)*S);
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){
const double inv2j = inv2[cell_id];
double inv2_min = inv2j;
double inv2_max = inv2j;
const auto cell_stencil = stencil[cell_id];
for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
inv2_min = std::min(inv2_min, inv2[cell_stencil[i_cell]]);
inv2_max = std::max(inv2_max, inv2[cell_stencil[i_cell]]);
}
double inv2_bar_min = inv2j;
double inv2_bar_max = inv2j;
const auto& cell_nodes = cell_to_node_matrix[cell_id];
for(size_t r = 0; r < cell_nodes.size(); ++r){
const double inv2_xr = 0.5*trace(transpose(DPk_Sh[cell_id](xr[cell_nodes[r]]))*DPk_Sh[cell_id](xr[cell_nodes[r]]));
inv2_bar_min = std::min(inv2_bar_min, inv2_xr);
inv2_bar_max = std::max(inv2_bar_max, inv2_xr);
}
//for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
// const CellId cell_k_id = cell_stencil[i_cell];
// const double inv2_xk = 0.5*trace(transpose(DPk_Sh[cell_id](xj[cell_k_id]))*DPk_Sh[cell_id](xj[cell_k_id]));
//
// inv2_bar_min = std::min(inv2_bar_min, inv2_xk);
// inv2_bar_max = std::max(inv2_bar_max, inv2_xk);
//}
const double eps = 1E-14;
double coef1 = 1.;
double coef2 = 1.;
if(std::abs(inv2_bar_max - inv2j) > eps){
coef1 = (inv2_max - inv2j) / (inv2_bar_max - inv2j);
}
if(std::abs(inv2_bar_min - inv2j) > eps){
coef2 = (inv2_min - inv2j) / (inv2_bar_min - inv2j);
}
//const double lambda_inv2 = sqrt(std::max(0., std::min(1., std::min(coef1, coef2))));
const double lambda_inv2 = std::max(0., std::min(1., std::min(coef1, coef2)));
lambda[cell_id] = lambda_inv2;
//auto coefficients = DPk_Sh.coefficients(cell_id);
//coefficients[0] = (1-lambda_inv2)*S[cell_id] + lambda_inv2 * coefficients[0];
//for (size_t i = 1; i < coefficients.size(); ++i) {
// coefficients[i] *= lambda_inv2;
//}
});
return lambda;
}
std::tuple<const std::shared_ptr<const DiscreteFunctionVariant>,
const std::shared_ptr<const DiscreteFunctionVariant>,
const std::shared_ptr<const DiscreteFunctionVariant>,
const std::shared_ptr<const DiscreteFunctionVariant>>
_apply_NeoHook(const DiscreteScalarFunction rho,
const DiscreteVectorFunction u,
const DiscreteScalarFunction E,
const DiscreteTensorFunction CG,
const DiscreteScalarFunction chi_solid,
const double& lambda,
const double& mu,
const double& gamma,
const double& p_inf) const
{
const TinyMatrix<Dimension> I = identity;
const DiscreteScalarFunction& eps = E - 0.5 * dot(u, u);
const DiscreteScalarFunction& p_d = (1 - chi_solid) * (gamma - 1) * rho * eps - p_inf * gamma;
const DiscreteTensorFunction& sigma_d =
chi_solid * (mu / sqrt(det(CG)) * (CG - I) + lambda / sqrt(det(CG)) * log(sqrt(det(CG))) * I) - p_d * I;
const DiscreteScalarFunction& aL_d = sqrt(chi_solid * (lambda + 2 * mu) / rho + gamma * (p_d + p_inf) / rho);
const DiscreteScalarFunction& aT_d = sqrt(chi_solid * mu / rho);
return {std::make_shared<DiscreteFunctionVariant>(sigma_d),
std::make_shared<DiscreteFunctionVariant>(p_d),
std::make_shared<DiscreteFunctionVariant>(aL_d),
std::make_shared<DiscreteFunctionVariant>(aT_d)};
}
std::tuple<const std::shared_ptr<const DiscreteFunctionVariant>,
const std::shared_ptr<const DiscreteFunctionVariant>,
const std::shared_ptr<const DiscreteFunctionVariant>,
const std::shared_ptr<const DiscreteFunctionVariant>>
_apply_NeoHook2(const DiscreteScalarFunction rho,
const DiscreteVectorFunction u,
const DiscreteScalarFunction E,
const DiscreteTensorFunction CG,
const DiscreteScalarFunction chi_solid,
const double& mu,
const double& gamma_f,
const double& p_inf_f,
const double& gamma_s,
const double& p_inf_s) const
{
const TinyMatrix<Dimension> I = identity;
const DiscreteScalarFunction& eps = E - 0.5 * dot(u, u);
const DiscreteScalarFunction& p_d =
(1 - chi_solid) * ((gamma_f - 1) * rho * eps - gamma_f * p_inf_f) + chi_solid * ((gamma_s - 1.) * rho * eps - gamma_s * p_inf_s);
const DiscreteTensorFunction& sigma_d =
chi_solid * 2 * mu / sqrt(det(CG)) * (CG - 1. / 3. * trace(CG) * I) - p_d * I;
const DiscreteScalarFunction& aL_d = sqrt(
chi_solid * (2 * mu) / rho + ((1 - chi_solid) * gamma_f * p_d + chi_solid * (gamma_s * (p_d + p_inf_s))) / rho);
const DiscreteScalarFunction& aT_d = sqrt(chi_solid * mu / rho);
return {std::make_shared<DiscreteFunctionVariant>(sigma_d),
std::make_shared<DiscreteFunctionVariant>(p_d),
std::make_shared<DiscreteFunctionVariant>(aL_d),
std::make_shared<DiscreteFunctionVariant>(aT_d)};
}
std::tuple<const std::shared_ptr<const DiscreteFunctionVariant>,
const std::shared_ptr<const DiscreteFunctionVariant>,
const std::shared_ptr<const DiscreteFunctionVariant>,
const std::shared_ptr<const DiscreteFunctionVariant>>
_apply_state_law(const size_t law,
const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
const std::shared_ptr<const DiscreteFunctionVariant>& E_v, const std::shared_ptr<const DiscreteFunctionVariant>& CG_v,
const std::shared_ptr<const DiscreteFunctionVariant>& chi_solid_v,
const double& lambda,
const double& mu,
const double& gamma_f,
const double& p_inf_f,
const double& gamma_s,
const double& p_inf_s) const
{
const DiscreteScalarFunction& chi_solid_d = chi_solid_v->get<DiscreteScalarFunction>();
const DiscreteScalarFunction& rho_d = rho_v->get<DiscreteScalarFunction>();
const DiscreteVectorFunction& u_d = u_v->get<DiscreteVectorFunction>();
const DiscreteScalarFunction& E_d = E_v->get<DiscreteScalarFunction>();
const DiscreteTensorFunction& CG_d = CG_v->get<DiscreteTensorFunction>();
// const std::shared_ptr<const DiscreteFunctionVariant> sigma;
// const std::shared_ptr<const DiscreteFunctionVariant> aL;
// const std::shared_ptr<const DiscreteFunctionVariant> aT;
if (law == 1) {
return _apply_NeoHook(rho_d, u_d, E_d, CG_d, chi_solid_d, lambda, mu, gamma_f, p_inf_f);
} else {
return _apply_NeoHook2(rho_d, u_d, E_d, CG_d, chi_solid_d, mu, gamma_f, p_inf_f, gamma_s, p_inf_s);
}
}
public:
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,
const std::shared_ptr<const DiscreteFunctionVariant>& aT_v,
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) const
{
std::shared_ptr mesh_v = getCommonMesh({rho_v, aL_v, aT_v, u_v, sigma_v});
if (not mesh_v) {
throw NormalError("discrete functions are not defined on the same mesh");
}
if (not checkDiscretizationType({rho_v, aL_v, u_v, sigma_v}, DiscreteFunctionType::P0)) {
throw NormalError("hyperelastic solver expects P0 functions");
}
const MeshType& mesh = *mesh_v->get<MeshType>();
const DiscreteScalarFunction& rho = rho_v->get<DiscreteScalarFunction>();
const DiscreteVectorFunction& u = u_v->get<DiscreteVectorFunction>();
const DiscreteScalarFunction& aL = aL_v->get<DiscreteScalarFunction>();
const DiscreteScalarFunction& aT = aT_v->get<DiscreteScalarFunction>();
const DiscreteTensorFunction& sigma = sigma_v->get<DiscreteTensorFunction>();
NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * aL, rho * aT);
NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
NodeValue<Rd> br = this->_computeBr(mesh, Ajr, u, sigma);
const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
this->_applyBoundaryConditions(bc_list, mesh, Ar, br);
synchronize(Ar);
synchronize(br);
NodeValue<const Rd> ur = this->_computeUr(mesh, Ar, br);
NodeValuePerCell<const Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u, sigma);
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>>
compute_fluxes(const SolverType& solver_type,
const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
const std::shared_ptr<const DiscreteFunctionVariant>& aL_v,
const std::shared_ptr<const DiscreteFunctionVariant>& aT_v,
const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
const std::shared_ptr<const DiscreteFunctionVariant>& sigma_v,
const std::shared_ptr<const DiscreteFunctionVariant>& p_v,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
{
std::shared_ptr mesh_v = getCommonMesh({rho_v, aL_v, aT_v, u_v, sigma_v, p_v});
if (not mesh_v) {
throw NormalError("discrete functions are not defined on the same mesh");
} // - p* identity; // - chi_solid*P_zero_air_repo*identity;
if (not checkDiscretizationType({rho_v, aL_v, u_v, sigma_v, p_v}, DiscreteFunctionType::P0)) {
throw NormalError("hyperelastic solver expects P0 functions");
}
const MeshType& mesh = *mesh_v->get<MeshType>();
const DiscreteScalarFunction& rho = rho_v->get<DiscreteScalarFunction>();
const DiscreteVectorFunction& u = u_v->get<DiscreteVectorFunction>();
const DiscreteScalarFunction& aL = aL_v->get<DiscreteScalarFunction>();
const DiscreteScalarFunction& aT = aT_v->get<DiscreteScalarFunction>();
const DiscreteTensorFunction& sigma = sigma_v->get<DiscreteTensorFunction>();
const DiscreteScalarFunction& p = p_v->get<DiscreteScalarFunction>();
const TinyMatrix<Dimension> I = identity;
const DiscreteTensorFunction& S = sigma + p*I;
const std::shared_ptr<const DiscreteFunctionVariant>& S_v = std::make_shared<DiscreteFunctionVariant>(S);
std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list;
for (auto&& bc_descriptor : bc_descriptor_list) {
if (bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry) {
symmetry_boundary_descriptor_list.push_back(bc_descriptor->boundaryDescriptor_shared());
}
}
PolynomialReconstructionDescriptor reconstruction_descriptor(IntegrationMethodType::cell_center, 1,
symmetry_boundary_descriptor_list);
auto reconstruction = PolynomialReconstruction{reconstruction_descriptor}.build(u_v, S_v, p_v);
auto DPk_uh = reconstruction[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
auto DPk_Sh = reconstruction[1]->get<DiscreteFunctionDPk<Dimension, const Rdxd>>();
auto DPk_ph = reconstruction[2]->get<DiscreteFunctionDPk<Dimension, const double>>();
DiscreteFunctionDPk<Dimension, Rd> u_lim = copy(DPk_uh);
DiscreteFunctionDPk<Dimension, Rdxd> S_lim = copy(DPk_Sh);
DiscreteFunctionDPk<Dimension, double> p_lim = copy(DPk_ph);
this->_vector_limiter(mesh, u, u_lim, DPk_ph);
CellValue<double> lambda = this->_matrix_limiter(mesh, S, S_lim);
this->_scalar_limiter(mesh, p, p_lim);
const NodeValue<const Rd>& xr = mesh.xr();
const DiscreteScalarFunction& J2 = 0.5*trace(transpose(S)*S);
NodeValuePerCell<double> coef{mesh.connectivity()};
const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
parallel_for(
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){
const double J2_p = 0.5*trace(transpose(DPk_Sh[j](xr[cell_nodes[r]]))*DPk_Sh[j](xr[cell_nodes[r]]));
if(J2_p != 0.){
coef(j,r) = J2[j]/J2_p;
}else{
coef(j,r) = 0;
}
} });
NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * aL, rho * aT);
NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
NodeValue<Rd> br = this->_computeBr(mesh, Ajr, u_lim, S_lim, p_lim, lambda, coef);
const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
this->_applyBoundaryConditions(bc_list, mesh, Ar, br);
synchronize(Ar);
synchronize(br);
NodeValue<const Rd> ur = this->_computeUr(mesh, Ar, br);
NodeValuePerCell<const Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u_lim, S_lim, p_lim, lambda, coef);
return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
std::make_shared<const SubItemValuePerItemVariant>(Fjr));
}
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,
const DiscreteVectorFunction& u,
const DiscreteScalarFunction& E,
const DiscreteTensorFunction& CG,
const NodeValue<const Rd>& ur_o2,
const NodeValue<const Rd>& ur_o1,
const NodeValuePerCell<const Rd>& Fjr_o2,
const NodeValuePerCell<const Rd>& Fjr_o1) const
{
const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
NodeValue<Rd> ur = copy(ur_o2);
NodeValuePerCell<Rd> Fjr = copy(Fjr_o2);
if ((mesh.shared_connectivity() != ur.connectivity_ptr()) or
(mesh.shared_connectivity() != Fjr.connectivity_ptr())) {
throw NormalError("fluxes are not defined on the same connectivity than the mesh");
}
StencilManager::BoundaryDescriptorList symmetry_boundary_descriptor_list;
StencilDescriptor stencil_descriptor{1, StencilDescriptor::ConnectionType::by_nodes};
auto stencil = StencilManager::instance().getCellToCellStencilArray(mesh.connectivity(), stencil_descriptor, symmetry_boundary_descriptor_list);
CellValue<const double> Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
const NodeValuePerCell<const Rd> Cjr = MeshDataManager::instance().getMeshData(mesh).Cjr();
CellValue<double> new_rho = copy(rho.cellValues());
CellValue<Rd> new_u = copy(u.cellValues());
CellValue<double> new_E = copy(E.cellValues());
CellValue<Rdxd> new_CG = copy(CG.cellValues());
CellValue<double> new_rho_stencil = copy(rho.cellValues());
CellValue<Rd> new_u_stencil = copy(u.cellValues());
CellValue<double> new_E_stencil = copy(E.cellValues());
CellValue<Rdxd> new_CG_stencil = copy(CG.cellValues());
CellValue<int> cell_flag{mesh.connectivity()};
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
auto cell_nodes = cell_to_node_matrix[j];
Rd momentum_fluxes = zero; double energy_fluxes = 0;
Rdxd gradv = zero;
for (size_t R = 0; R < cell_nodes.size(); ++R) {
NodeId r = cell_nodes[R];
gradv += tensorProduct(ur[r], Cjr(j, R));
momentum_fluxes += Fjr(j, R);
energy_fluxes += dot(Fjr(j, R), ur[r]);
}
Rdxd cauchy_green_fluxes = gradv * CG[j] + CG[j] * transpose(gradv);
double dt_over_Mj = dt / (rho[j] * Vj[j]);
double dt_over_Vj = dt / Vj[j];
new_u[j] += dt_over_Mj * momentum_fluxes;
new_E[j] += dt_over_Mj * energy_fluxes;
new_CG[j] += dt_over_Vj * cauchy_green_fluxes;
new_CG[j] += transpose(new_CG[j]);
new_CG[j] *= 0.5;
double new_eps = new_E[j] - 0.5 * dot(new_u[j],new_u[j]);
if(new_eps < 0.){
cell_flag[j] = 1;
}else{
cell_flag[j] = 0;
}
});
for(CellId j = 0; j < mesh.numberOfCells(); ++j){
if(cell_flag[j] == 1){
const auto& cell_nodes = cell_to_node_matrix[j];
Rd momentum_fluxes = zero;
double energy_fluxes = 0;
Rdxd gradv = zero;
std::vector<NodeId> NodeList;
for (size_t R = 0; R < cell_nodes.size(); ++R) {
NodeId r = cell_nodes[R];
NodeList.push_back(r);
Fjr(j,R) = Fjr_o1(j,R);
ur[r] = ur_o1[r];
gradv += tensorProduct(ur[r], Cjr(j, R));
momentum_fluxes += Fjr(j, R);
energy_fluxes += dot(Fjr(j, R), ur[r]);
}
Rdxd cauchy_green_fluxes = gradv * CG[j] + CG[j] * transpose(gradv);
double dt_over_Mj = dt / (rho[j] * Vj[j]);
double dt_over_Vj = dt / Vj[j];
new_u[j] = new_u_stencil[j] + dt_over_Mj * momentum_fluxes;
new_E[j] = new_E_stencil[j] + dt_over_Mj * energy_fluxes;
new_CG[j] = new_CG_stencil[j] + dt_over_Vj * cauchy_green_fluxes;
new_CG[j] += transpose(new_CG[j]);
new_CG[j] *= 0.5;
auto cell_stencil = stencil[j];
for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
CellId J = cell_stencil[i_cell];
auto i_cell_nodes = cell_to_node_matrix[J];
momentum_fluxes = zero;
energy_fluxes = 0;
gradv = zero;
for (size_t R = 0; R < i_cell_nodes.size(); ++R) {
NodeId i_node = i_cell_nodes[R];
for(size_t node_test = 0; node_test < NodeList.size(); ++node_test){
if(NodeList[node_test] == i_node){
Fjr(J,R) = Fjr_o1(J,R);
}
}
gradv += tensorProduct(ur[i_node], Cjr(J, R));
momentum_fluxes += Fjr(J, R);
energy_fluxes += dot(Fjr(J, R), ur[i_node]);
}
cauchy_green_fluxes = gradv * CG[J] + CG[J] * transpose(gradv); dt_over_Mj = dt / (rho[J] * Vj[J]);
dt_over_Vj = dt / Vj[J];
new_u[J] = new_u_stencil[J] + dt_over_Mj * momentum_fluxes;
new_E[J] = new_E_stencil[J] + dt_over_Mj * energy_fluxes;
new_CG[J] = new_CG_stencil[J] + dt_over_Vj * cauchy_green_fluxes;
new_CG[J] += transpose(new_CG[J]);
new_CG[J] *= 0.5;
}
}
}
NodeValue<Rd> new_xr = copy(mesh.xr());
parallel_for(
mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });
std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr);
CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { new_rho[j] *= Vj[j] / new_Vj[j]; });
return {std::make_shared<MeshVariant>(new_mesh),
std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_rho)),
std::make_shared<DiscreteFunctionVariant>(DiscreteVectorFunction(new_mesh, new_u)),
std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_E)),
std::make_shared<DiscreteFunctionVariant>(DiscreteTensorFunction(new_mesh, new_CG))};
}
std::tuple<std::shared_ptr<const MeshVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
apply_fluxes(const double& dt,
const std::shared_ptr<const DiscreteFunctionVariant>& rho,
const std::shared_ptr<const DiscreteFunctionVariant>& u,
const std::shared_ptr<const DiscreteFunctionVariant>& E,
const std::shared_ptr<const DiscreteFunctionVariant>& CG,
const std::shared_ptr<const ItemValueVariant>& ur,
const std::shared_ptr<const ItemValueVariant>& ur_o1,
const std::shared_ptr<const SubItemValuePerItemVariant> Fjr,
const std::shared_ptr<const SubItemValuePerItemVariant> Fjr_o1) const
{
std::shared_ptr mesh_v = getCommonMesh({rho, u, E});
if (not mesh_v) {
throw NormalError("discrete functions are not defined on the same mesh");
}
if (not checkDiscretizationType({rho, u, E}, DiscreteFunctionType::P0)) {
throw NormalError("hyperelastic solver expects P0 functions");
}
return this->apply_fluxes(dt, //
*mesh_v->get<MeshType>(), //
rho->get<DiscreteScalarFunction>(), //
u->get<DiscreteVectorFunction>(), //
E->get<DiscreteScalarFunction>(), //
CG->get<DiscreteTensorFunction>(), //
ur->get<NodeValue<const Rd>>(),
ur_o1->get<NodeValue<const Rd>>(), //
Fjr->get<NodeValuePerCell<const Rd>>(),
Fjr_o1->get<NodeValuePerCell<const Rd>>()
);
}
std::tuple<std::shared_ptr<const MeshVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>>
order2_apply_fluxes(const double& dt,
const MeshType& mesh,
const DiscreteScalarFunction& rho,
const DiscreteVectorFunction& u,
const DiscreteScalarFunction& E,
const DiscreteTensorFunction& CG,
const DiscreteScalarFunction& rho_app,
const DiscreteTensorFunction& CG_app,
const NodeValue<const Rd>& ur_o2,
const NodeValue<const Rd>& ur_o1,
const NodeValuePerCell<const Rd>& Fjr_o2,
const NodeValuePerCell<const Rd>& Fjr_o1) const
{
const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
NodeValue<Rd> ur = copy(ur_o2);
NodeValuePerCell<Rd> Fjr = copy(Fjr_o2);
if ((mesh.shared_connectivity() != ur.connectivity_ptr()) or
(mesh.shared_connectivity() != Fjr.connectivity_ptr())) {
throw NormalError("fluxes are not defined on the same connectivity than the mesh");
}
StencilManager::BoundaryDescriptorList symmetry_boundary_descriptor_list;
StencilDescriptor stencil_descriptor{1, StencilDescriptor::ConnectionType::by_nodes};
auto stencil = StencilManager::instance().getCellToCellStencilArray(mesh.connectivity(), stencil_descriptor, symmetry_boundary_descriptor_list);
CellValue<const double> Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
const NodeValuePerCell<const Rd> Cjr = MeshDataManager::instance().getMeshData(mesh).Cjr();
CellValue<double> new_rho = copy(rho.cellValues());
CellValue<Rd> new_u = copy(u.cellValues());
CellValue<double> new_E = copy(E.cellValues());
CellValue<Rdxd> new_CG = copy(CG.cellValues());
CellValue<double> new_rho_stencil = copy(rho.cellValues());
CellValue<Rd> new_u_stencil = copy(u.cellValues());
CellValue<double> new_E_stencil = copy(E.cellValues());
CellValue<Rdxd> new_CG_stencil = copy(CG.cellValues());
CellValue<int> cell_flag{mesh.connectivity()};
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
const auto& cell_nodes = cell_to_node_matrix[j];
Rd momentum_fluxes = zero;
double energy_fluxes = 0;
Rdxd gradv = zero;
for (size_t R = 0; R < cell_nodes.size(); ++R) {
const NodeId r = cell_nodes[R];
gradv += tensorProduct(ur[r], Cjr(j, R));
momentum_fluxes += Fjr(j, R);
energy_fluxes += dot(Fjr(j, R), ur[r]);
}
const Rdxd cauchy_green_fluxes = rho_app[j] * (gradv * CG_app[j] + CG_app[j] * transpose(gradv));
const double dt_over_Mj = dt / (rho[j] * Vj[j]);
new_u[j] += dt_over_Mj * momentum_fluxes;
new_E[j] += dt_over_Mj * energy_fluxes;
new_CG[j] += dt_over_Mj * cauchy_green_fluxes;
new_CG[j] += transpose(new_CG[j]);
new_CG[j] *= 0.5;
const double& new_eps = new_E[j] - 0.5 * dot(new_u[j],new_u[j]);
if(new_eps < 0.){
cell_flag[j] = 1;
}else{
cell_flag[j] = 0;
} });
for(CellId j = 0; j < mesh.numberOfCells(); ++j){
if(cell_flag[j] == 1){
const auto& cell_nodes = cell_to_node_matrix[j];
Rd momentum_fluxes = zero;
double energy_fluxes = 0;
Rdxd gradv = zero;
std::vector<NodeId> NodeList;
for (size_t R = 0; R < cell_nodes.size(); ++R) {
NodeId r = cell_nodes[R];
NodeList.push_back(r);
Fjr(j,R) = Fjr_o1(j,R);
ur[r] = ur_o1[r];
gradv += tensorProduct(ur[r], Cjr(j, R));
momentum_fluxes += Fjr(j, R);
energy_fluxes += dot(Fjr(j, R), ur[r]);
}
Rdxd cauchy_green_fluxes = rho_app[j] * (gradv * CG_app[j] + CG_app[j] * transpose(gradv));
double dt_over_Mj = dt / (rho[j] * Vj[j]);
new_u[j] = new_u_stencil[j] + dt_over_Mj * momentum_fluxes;
new_E[j] = new_E_stencil[j] + dt_over_Mj * energy_fluxes;
new_CG[j] = new_CG_stencil[j] + dt_over_Mj * cauchy_green_fluxes;
new_CG[j] += transpose(new_CG[j]);
new_CG[j] *= 0.5;
auto cell_stencil = stencil[j];
for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
CellId J = cell_stencil[i_cell];
auto i_cell_nodes = cell_to_node_matrix[J];
momentum_fluxes = zero;
energy_fluxes = 0;
gradv = zero;
for (size_t R = 0; R < i_cell_nodes.size(); ++R) {
NodeId i_node = i_cell_nodes[R];
for(size_t node_test = 0; node_test < NodeList.size(); ++node_test){
if(NodeList[node_test] == i_node){
Fjr(J,R) = Fjr_o1(J,R);
}
}
gradv += tensorProduct(ur[i_node], Cjr(J, R));
momentum_fluxes += Fjr(J, R);
energy_fluxes += dot(Fjr(J, R), ur[i_node]);
}
cauchy_green_fluxes = rho_app[J] * (gradv * CG_app[J] + CG_app[J] * transpose(gradv));
dt_over_Mj = dt / (rho[J] * Vj[J]);
new_u[J] = new_u_stencil[J] + dt_over_Mj * momentum_fluxes;
new_E[J] = new_E_stencil[J] + dt_over_Mj * energy_fluxes;
new_CG[J] = new_CG_stencil[J] + dt_over_Mj * cauchy_green_fluxes;
new_CG[J] += transpose(new_CG[J]);
new_CG[J] *= 0.5;
}
}
}
NodeValue<Rd> new_xr = copy(mesh.xr());
parallel_for(
mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });
std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr);
CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { new_rho[j] *= Vj[j] / new_Vj[j]; });
return {std::make_shared<MeshVariant>(new_mesh),
std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_rho)),
std::make_shared<DiscreteFunctionVariant>(DiscreteVectorFunction(new_mesh, new_u)), std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_E)),
std::make_shared<DiscreteFunctionVariant>(DiscreteTensorFunction(new_mesh, new_CG))};
}
std::tuple<std::shared_ptr<const MeshVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
order2_apply_fluxes(const double& dt,
const std::shared_ptr<const DiscreteFunctionVariant>& rho,
const std::shared_ptr<const DiscreteFunctionVariant>& u,
const std::shared_ptr<const DiscreteFunctionVariant>& E,
const std::shared_ptr<const DiscreteFunctionVariant>& CG,
const std::shared_ptr<const DiscreteFunctionVariant>& rho_app,
const std::shared_ptr<const DiscreteFunctionVariant>& CG_app,
const std::shared_ptr<const ItemValueVariant>& ur,
const std::shared_ptr<const ItemValueVariant>& ur_o1,
const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr,
const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr_o1) const
{
std::shared_ptr mesh_v = getCommonMesh({rho, u, E});
if (not mesh_v) {
throw NormalError("discrete functions are not defined on the same mesh");
}
if (not checkDiscretizationType({rho, u, E}, DiscreteFunctionType::P0)) {
throw NormalError("hyperelastic solver expects P0 functions");
}
return this->order2_apply_fluxes(dt, //
*mesh_v->get<MeshType>(), //
rho->get<DiscreteScalarFunction>(), //
u->get<DiscreteVectorFunction>(), //
E->get<DiscreteScalarFunction>(), //
CG->get<DiscreteTensorFunction>(), //
rho_app->get<DiscreteScalarFunction>(), //
CG_app->get<DiscreteTensorFunction>(), //
ur->get<NodeValue<const Rd>>(), //
ur_o1->get<NodeValue<const Rd>>(),
Fjr->get<NodeValuePerCell<const Rd>>(),
Fjr_o1->get<NodeValuePerCell<const Rd>>());
}
std::tuple<std::shared_ptr<const MeshVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
apply(const SolverType& solver_type,
const size_t& law,
const double& dt,
const std::shared_ptr<const DiscreteFunctionVariant>& rho,
const std::shared_ptr<const DiscreteFunctionVariant>& u,
const std::shared_ptr<const DiscreteFunctionVariant>& E,
const std::shared_ptr<const DiscreteFunctionVariant>& CG,
const std::shared_ptr<const DiscreteFunctionVariant>& aL,
const std::shared_ptr<const DiscreteFunctionVariant>& aT,
const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
const std::shared_ptr<const DiscreteFunctionVariant>& p,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const std::shared_ptr<const DiscreteFunctionVariant>& chi_solid,
const double& lambda,
const double& mu,
const double& gamma_f,
const double& p_inf_f,
const double& gamma_s,
const double& p_inf_s) const
{
//std::shared_ptr m_app = getCommonMesh({rho, aL, aT, u, sigma}); //std::shared_ptr<const DiscreteFunctionVariant> rho_app = rho;
//std::shared_ptr<const DiscreteFunctionVariant> u_app = u;
//std::shared_ptr<const DiscreteFunctionVariant> E_app = E;
//std::shared_ptr<const DiscreteFunctionVariant> CG_app = CG;
auto [ur, Fjr] = compute_fluxes(solver_type, rho, aL, aT, u, sigma, p, bc_descriptor_list);
auto [ur_o1, Fjr_o1] = compute_fluxes(solver_type, rho, aL, aT, u, sigma, bc_descriptor_list);
auto [m_app, rho_app, u_app, E_app, CG_app] = apply_fluxes(dt/2., rho, u, E, CG, ur, ur_o1, Fjr, Fjr_o1);
std::shared_ptr<const DiscreteFunctionVariant> chi_solid_app = shallowCopy(m_app, chi_solid);
auto [sigma_app, p_app, aL_app, aT_app] =
_apply_state_law(law,rho_app,u_app,E_app,CG_app,chi_solid_app,lambda,mu,gamma_f,p_inf_f,gamma_s,p_inf_s);
auto [ur_app, Fjr_app] = compute_fluxes(solver_type, rho_app, aL_app,aT_app,u_app,sigma_app, p_app, bc_descriptor_list);
auto [ur_o1_app, Fjr_o1_app] = compute_fluxes(solver_type, rho_app, aL_app,aT_app,u_app,sigma_app,bc_descriptor_list);
return order2_apply_fluxes(dt, rho, u, E, CG, rho_app, CG_app, ur_app, ur_o1_app, Fjr_app, Fjr_o1_app);
//return apply_fluxes(dt, rho, u, E, CG, ur, ur_o1, Fjr, Fjr_o1);
}
Order2HyperelasticSolver() = default;
Order2HyperelasticSolver(Order2HyperelasticSolver&&) = default;
~Order2HyperelasticSolver() = default;
};
template <MeshConcept MeshType>
void
Order2HyperelasticSolverHandler::Order2HyperelasticSolver<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>) {
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
if constexpr (Dimension == 1) {
const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
const auto& node_list = bc.nodeList();
const auto& value_list = bc.valueList();
parallel_for(
node_list.size(), PUGS_LAMBDA(size_t i_node) {
const NodeId node_id = node_list[i_node];
const auto& node_cell_list = node_to_cell_matrix[node_id];
Assert(node_cell_list.size() == 1);
CellId node_cell_id = node_cell_list[0];
size_t node_local_number_in_cell = node_local_numbers_in_their_cells(node_id, 0);
br[node_id] -= value_list[i_node] * Cjr(node_cell_id, node_local_number_in_cell);
});
} else {
const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();
const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells();
const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed();
const auto& face_list = bc.faceList();
const auto& value_list = bc.valueList();
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face]; const auto& face_cell_list = face_to_cell_matrix[face_id];
Assert(face_cell_list.size() == 1);
CellId face_cell_id = face_cell_list[0];
size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);
const double sign = face_cell_is_reversed(face_cell_id, face_local_number_in_cell) ? -1 : 1;
const auto& face_nodes = face_to_node_matrix[face_id];
for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
NodeId node_id = face_nodes[i_node];
br[node_id] -= sign * value_list[i_face] * Nlr(face_id, i_node);
}
}
}
}
},
boundary_condition);
}
}
template <MeshConcept MeshType>
void
Order2HyperelasticSolverHandler::Order2HyperelasticSolver<MeshType>::_applyNormalStressBC(
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<NormalStressBoundaryCondition, T>) {
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
if constexpr (Dimension == 1) {
const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
const auto& node_list = bc.nodeList();
const auto& value_list = bc.valueList();
parallel_for(
node_list.size(), PUGS_LAMBDA(size_t i_node) {
const NodeId node_id = node_list[i_node];
const auto& node_cell_list = node_to_cell_matrix[node_id];
Assert(node_cell_list.size() == 1);
CellId node_cell_id = node_cell_list[0];
size_t node_local_number_in_cell = node_local_numbers_in_their_cells(node_id, 0);
br[node_id] += value_list[i_node] * Cjr(node_cell_id, node_local_number_in_cell);
});
} else {
const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();
const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells();
const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed();
const auto& face_list = bc.faceList();
const auto& value_list = bc.valueList();
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
const auto& face_cell_list = face_to_cell_matrix[face_id];
Assert(face_cell_list.size() == 1);
CellId face_cell_id = face_cell_list[0]; size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);
const double sign = face_cell_is_reversed(face_cell_id, face_local_number_in_cell) ? -1 : 1;
const auto& face_nodes = face_to_node_matrix[face_id];
for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
NodeId node_id = face_nodes[i_node];
br[node_id] += sign * value_list[i_face] * Nlr(face_id, i_node);
}
}
}
}
},
boundary_condition);
}
}
template <MeshConcept MeshType>
void
Order2HyperelasticSolverHandler::Order2HyperelasticSolver<MeshType>::_applySymmetryBC(
const BoundaryConditionList& bc_list,
NodeValue<Rdxd>& Ar,
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<SymmetryBoundaryCondition, T>) {
const Rd& n = bc.outgoingNormal();
const Rdxd I = identity;
const Rdxd nxn = tensorProduct(n, n);
const Rdxd P = I - nxn;
const Array<const NodeId>& node_list = bc.nodeList();
parallel_for(
bc.numberOfNodes(), PUGS_LAMBDA(int r_number) {
const NodeId r = node_list[r_number];
Ar[r] = P * Ar[r] * P + nxn;
br[r] = P * br[r];
});
}
},
boundary_condition);
}
}
template <MeshConcept MeshType>
void
Order2HyperelasticSolverHandler::Order2HyperelasticSolver<MeshType>::_applyVelocityBC(
const BoundaryConditionList& bc_list,
NodeValue<Rdxd>& Ar,
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<VelocityBoundaryCondition, T>) {
const auto& node_list = bc.nodeList();
const auto& value_list = bc.valueList();
parallel_for(
node_list.size(), PUGS_LAMBDA(size_t i_node) {
NodeId node_id = node_list[i_node];
const auto& value = value_list[i_node];
Ar[node_id] = identity;
br[node_id] = value;
});
} else if constexpr (std::is_same_v<FixedBoundaryCondition, T>) {
const auto& node_list = bc.nodeList();
parallel_for(
node_list.size(), PUGS_LAMBDA(size_t i_node) {
NodeId node_id = node_list[i_node];
Ar[node_id] = identity;
br[node_id] = zero;
});
}
},
boundary_condition);
}
}
template <MeshConcept MeshType>
class Order2HyperelasticSolverHandler::Order2HyperelasticSolver<MeshType>::FixedBoundaryCondition
{
private:
const MeshNodeBoundary m_mesh_node_boundary;
public:
const Array<const NodeId>&
nodeList() const
{
return m_mesh_node_boundary.nodeList();
}
FixedBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary) : m_mesh_node_boundary{mesh_node_boundary} {}
~FixedBoundaryCondition() = default;
};
template <MeshConcept MeshType>
class Order2HyperelasticSolverHandler::Order2HyperelasticSolver<MeshType>::PressureBoundaryCondition
{
private:
const MeshFaceBoundary m_mesh_face_boundary;
const Array<const double> m_value_list;
public:
const Array<const FaceId>&
faceList() const
{
return m_mesh_face_boundary.faceList();
}
const Array<const double>&
valueList() const
{
return m_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}
{}
~PressureBoundaryCondition() = default;
};
template <>
class Order2HyperelasticSolverHandler::Order2HyperelasticSolver<Mesh<1>>::PressureBoundaryCondition
{
private:
const MeshNodeBoundary m_mesh_node_boundary;
const Array<const double> m_value_list;
public:
const Array<const NodeId>&
nodeList() const
{
return m_mesh_node_boundary.nodeList();
}
const Array<const double>&
valueList() const
{
return m_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 <MeshConcept MeshType>
class Order2HyperelasticSolverHandler::Order2HyperelasticSolver<MeshType>::NormalStressBoundaryCondition
{
private:
const MeshFaceBoundary m_mesh_face_boundary;
const Array<const Rdxd> m_value_list;
public:
const Array<const FaceId>&
faceList() const
{
return m_mesh_face_boundary.faceList();
}
const Array<const Rdxd>&
valueList() const
{
return m_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}
{}
~NormalStressBoundaryCondition() = default;
};
template <>
class Order2HyperelasticSolverHandler::Order2HyperelasticSolver<Mesh<1>>::NormalStressBoundaryCondition
{
private:
const MeshNodeBoundary m_mesh_node_boundary;
const Array<const Rdxd> m_value_list;
public:
const Array<const NodeId>&
nodeList() const
{
return m_mesh_node_boundary.nodeList();
}
const Array<const Rdxd>&
valueList() const
{
return m_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 <MeshConcept MeshType>
class Order2HyperelasticSolverHandler::Order2HyperelasticSolver<MeshType>::VelocityBoundaryCondition
{
private:
const MeshNodeBoundary m_mesh_node_boundary;
const Array<const TinyVector<Dimension>> m_value_list;
public:
const Array<const NodeId>&
nodeList() const
{
return m_mesh_node_boundary.nodeList();
}
const Array<const TinyVector<Dimension>>&
valueList() const
{
return m_value_list;
}
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}
{}
~VelocityBoundaryCondition() = default;
};
template <MeshConcept MeshType>
class Order2HyperelasticSolverHandler::Order2HyperelasticSolver<MeshType>::SymmetryBoundaryCondition
{
public:
using Rd = TinyVector<Dimension, double>;
private:
const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_boundary;
public:
const Rd&
outgoingNormal() const
{
return m_mesh_flat_node_boundary.outgoingNormal();
}
size_t
numberOfNodes() const
{
return m_mesh_flat_node_boundary.nodeList().size();
}
const Array<const NodeId>&
nodeList() const
{
return m_mesh_flat_node_boundary.nodeList();
}
SymmetryBoundaryCondition(const MeshFlatNodeBoundary<MeshType>& mesh_flat_node_boundary)
: m_mesh_flat_node_boundary(mesh_flat_node_boundary)
{
;
}
~SymmetryBoundaryCondition() = default;
};
Order2HyperelasticSolverHandler::Order2HyperelasticSolverHandler(const std::shared_ptr<const MeshVariant>& mesh_v)
{
if (not mesh_v) {
throw NormalError("discrete functions are not defined on the same mesh");
}
std::visit(
[&](auto&& mesh) {
using MeshType = mesh_type_t<decltype(mesh)>;
if constexpr (is_polygonal_mesh_v<MeshType>) {
m_hyperelastic_solver = std::make_unique<Order2HyperelasticSolver<MeshType>>();
} else {
throw NormalError("unexpected mesh type");
}
},
mesh_v->variant());
}