Select Git revision
Order2HyperelasticSolver.cpp
Order2HyperelasticSolver.cpp 59.80 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;
}
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_sigma) const
{
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
const NodeValue<const Rd>& xr = mesh.xr();
//const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix(); //const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
NodeValuePerCell<Rd> F{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);
//
// 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];
// F(J, R) = -Ajr(J, R) * (DPk_uh[J](xr[r]) - ur[r]) + DPk_sigma[J](xr[r]) * Cjr(J, R);
// }
// });
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) * (DPk_uh[j](xr[cell_nodes[r]]) - ur[cell_nodes[r]]) + DPk_sigma[j](xr[cell_nodes[r]])*Cjr(j,r);
}
});
return F;
}
void
_tensor_limiter(const MeshType& mesh,
const DiscreteFunctionP0<const Rdxd>& f,
DiscreteFunctionDPk<Dimension, Rdxd>& 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();
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
const Rdxd fj = f[cell_id];
Rdxd f_min = fj;
Rdxd f_max = fj;
const auto cell_stencil = stencil[cell_id];
for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
for (size_t row = 0; row < Dimension; ++row) {
for (size_t col = 0; col < Dimension; ++col) {
f_min(row, col) = std::min(f_min(row, col), f[cell_stencil[i_cell]](row, col));
f_max(row, col) = std::max(f_max(row, col), f[cell_stencil[i_cell]](row, col));
}
}
}
Rdxd f_bar_min = fj;
Rdxd 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 Rdxd f_xk = DPk_fh[cell_id](xj[cell_k_id]);
for (size_t row = 0; row < Dimension; ++row) {
for (size_t col = 0; col < Dimension; ++col) {
f_bar_min(row, col) = std::min(f_bar_min(row, col), f_xk(row, col));
f_bar_max(row, col) = std::max(f_bar_max(row, col), f_xk(row, col)); }
}
}
const double eps = 1E-14;
Rdxd coef1;
for (size_t row = 0; row < Dimension; ++row) {
for (size_t col = 0; col < Dimension; ++col) {
coef1(row, col) = 1;
if (std::abs(f_bar_max(row, col) - fj(row, col)) > eps) {
coef1(row, col) = (f_max(row, col) - fj(row, col)) / (f_bar_max(row, col) - fj(row, col));
}
}
}
Rdxd coef2;
for (size_t row = 0; row < Dimension; ++row) {
for (size_t col = 0; col < Dimension; ++col) {
coef2(row, col) = 1;
if (std::abs(f_bar_min(row, col) - fj(row, col)) > eps) {
coef2(row, col) = (f_min(row, col) - fj(row, col)) / ((f_bar_min(row, col) - fj(row, col)));
}
}
}
double min_coef1 = coef1(0, 0);
double min_coef2 = coef2(0, 0);
for (size_t row = 0; row < Dimension; ++row) {
for (size_t col = 0; col < Dimension; ++col) {
min_coef1 = std::min(min_coef1, coef1(row, col));
min_coef2 = std::min(min_coef2, coef2(row, col));
}
}
const double lambda = std::max(0., std::min(1., std::min(min_coef1, min_coef2)));
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
{
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();
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
const Rd fj = f[cell_id];
Rd f_min = fj;
Rd f_max = fj;
const auto cell_stencil = stencil[cell_id];
for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
for (size_t dim = 0; dim < Dimension; ++dim) {
f_min[dim] = std::min(f_min[dim], f[cell_stencil[i_cell]][dim]);
f_max[dim] = std::max(f_max[dim], f[cell_stencil[i_cell]][dim]); }
}
Rd f_bar_min = fj;
Rd 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 Rd f_xk = DPk_fh[cell_id](xj[cell_k_id]);
for (size_t dim = 0; dim < Dimension; ++dim) {
f_bar_min[dim] = std::min(f_bar_min[dim], f_xk[dim]);
f_bar_max[dim] = std::max(f_bar_max[dim], f_xk[dim]);
}
}
const double eps = 1E-14;
Rd coef1;
for (size_t dim = 0; dim < Dimension; ++dim) {
coef1[dim] = 1;
if (std::abs(f_bar_max[dim] - fj[dim]) > eps) {
coef1[dim] = (f_max[dim] - fj[dim]) / ((f_bar_max[dim] - fj[dim]));
}
}
Rd coef2;
for (size_t dim = 0; dim < Dimension; ++dim) {
coef2[dim] = 1;
if (std::abs(f_bar_min[dim] - fj[dim]) > eps) {
coef2[dim] = (f_min[dim] - fj[dim]) / ((f_bar_min[dim] - fj[dim]));
}
}
double min_coef1 = coef1[0];
double min_coef2 = coef2[0];
for (size_t dim = 1; dim < Dimension; ++dim) {
min_coef1 = std::min(min_coef1, coef1[dim]);
min_coef2 = std::min(min_coef2, coef2[dim]);
}
const double lambda = std::max(0., std::min(1., std::min(min_coef1, min_coef2)));
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;
}
});
}
CellValue<double>
_scalar_limiter_coefficient(const MeshType& mesh,
const DiscreteFunctionP0<const double>& f,
const DiscreteFunctionDPk<Dimension, const 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> lambda{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));
}
lambda[cell_id] = std::max(0., std::min(1., std::min(coef1, coef2)));
});
return lambda;
}
std::tuple<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>(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>>
_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) + 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>(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>>
_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");
} // - p* identity; // - chi_solid*P_zero_air_repo*identity;
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>();
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);
// CellValue<double> limiter_j{mesh.connectivity()};
// parallel_for(
// mesh.numberOfCells(), PUGS_LAMBDA(CellId j){
// limiter_j[j] = 1;
// }
//);
//
// const CellValue<const Rdxd>& sigma_j = copy(sigma.cellValues());
// std::vector<DiscreteFunctionDPk<Dimension, double>> sigma_coef;
// std::vector<size_t> row;
// std::vector<size_t> col;
// for(size_t i = 0; i < Dimension; ++i){
// for(size_t k = i; k < Dimension; ++k){
//
// CellValue<double> coef{mesh.connectivity()};
// for(CellId j = 0; j <mesh.numberOfCells(); ++j){
// coef[j] = sigma_j[j](i,k);
// }
//
// const auto& coef_scalar_function = DiscreteScalarFunction(mesh_v, coef);
// auto coef_v = std::make_shared<DiscreteFunctionVariant>(coef_scalar_function);
// auto reconstruction = PolynomialReconstruction{reconstruction_descriptor}.build(coef_v);
// auto DPk_coef = reconstruction[0]->get<DiscreteFunctionDPk<Dimension, const double>>();
//
// const CellValue<const double>& limiter_ik =
// _scalar_limiter_coefficient(mesh,coef_scalar_function,DPk_coef); parallel_for(
// mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
// limiter_j[j] = std::min(limiter_j[j],limiter_ik[j]);
// });
//
// sigma_coef.push_back(copy(DPk_coef));
// row.push_back(i);
// col.push_back(k);
// }
// }
//
// NodeValuePerCell<Rdxd> sigma_lim{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 i = 0; i < Dimension; ++i){
// for(size_t k = 0; k < Dimension; ++k){
// for(size_t R = 0; R < cell_nodes.size(); ++R){
// sigma_lim(j,R)(i,k) = 0;
// }
// }
// } // }
//);
//
// const NodeValue<const Rd>& xr = copy(mesh.xr());
// 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 NodeId r = cell_nodes[R];
//
// for(size_t l = 0; l < sigma_coef.size(); ++l){
// const size_t& i = row[l];
// const size_t& k = col[l];
//
// auto coefficients = sigma_coef[l].coefficients(j);
//
// coefficients[0] = (1-limiter_j[j])*sigma_j[j](i,k) + limiter_j[j] * coefficients[0];
// for (size_t indice = 1; indice < coefficients.size(); ++indice) {
// coefficients[indice] *= limiter_j[j];
// }
//
// sigma_lim(j,R)(i,k) = sigma_coef[l][j](xr[r]);
// if(i != k){
// sigma_lim(j,R)(i,k) *= 2;
// }
// }
//
// sigma_lim(j,R) += transpose(sigma_lim(j,R));
// sigma_lim(j,R) *= 0.5;
// }
//
// }
//);
auto reconstruction = PolynomialReconstruction{reconstruction_descriptor}.build(u_v, sigma_v);
auto DPk_uh = reconstruction[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
auto DPk_sigmah = reconstruction[1]->get<DiscreteFunctionDPk<Dimension, const Rdxd>>();
DiscreteFunctionDPk<Dimension, Rd> u_lim = copy(DPk_uh);
DiscreteFunctionDPk<Dimension, Rdxd> sigma_lim = copy(DPk_sigmah);
this->_vector_limiter(mesh, u, u_lim);
this->_tensor_limiter(mesh, sigma, sigma_lim);
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, sigma_lim);
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, sigma_lim);
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,
const NodeValuePerCell<const Rd>& Fjr) const
{
const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
if ((mesh.shared_connectivity() != ur.connectivity_ptr()) or
(mesh.shared_connectivity() != Fjr.connectivity_ptr())) {
throw NormalError("fluxes are not defined on the same connectivity than the mesh");
}
NodeValue<Rd> new_xr = copy(mesh.xr());
parallel_for(
mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });
std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr);
CellValue<const double> Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
const NodeValuePerCell<const Rd> Cjr = MeshDataManager::instance().getMeshData(mesh).Cjr();
CellValue<double> new_rho = copy(rho.cellValues());
CellValue<Rd> new_u = copy(u.cellValues());
CellValue<double> new_E = copy(E.cellValues());
CellValue<Rdxd> new_CG = copy(CG.cellValues());
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
const auto& cell_nodes = cell_to_node_matrix[j];
Rd momentum_fluxes = zero;
double energy_fluxes = 0;
Rdxd gradv = zero;
for (size_t R = 0; R < cell_nodes.size(); ++R) {
const NodeId r = cell_nodes[R];
gradv += tensorProduct(ur[r], Cjr(j, R));
momentum_fluxes += Fjr(j, R);
energy_fluxes += dot(Fjr(j, R), ur[r]);
}
const Rdxd cauchy_green_fluxes = gradv * CG[j] + CG[j] * transpose(gradv);
const double dt_over_Mj = dt / (rho[j] * Vj[j]);
const 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;
});
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 SubItemValuePerItemVariant>& Fjr) const
{
std::shared_ptr mesh_v = getCommonMesh({rho, u, E});
if (not mesh_v) {
throw NormalError("discrete functions are not defined on the same mesh");
}
if (not checkDiscretizationType({rho, u, E}, DiscreteFunctionType::P0)) {
throw NormalError("hyperelastic solver expects P0 functions");
}
return this->apply_fluxes(dt, //
*mesh_v->get<MeshType>(), //
rho->get<DiscreteScalarFunction>(), //
u->get<DiscreteVectorFunction>(), //
E->get<DiscreteScalarFunction>(), //
CG->get<DiscreteTensorFunction>(), //
ur->get<NodeValue<const Rd>>(), //
Fjr->get<NodeValuePerCell<const Rd>>());
}
std::tuple<std::shared_ptr<const MeshVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
order2_apply_fluxes(const double& dt,
const MeshType& mesh,
const DiscreteScalarFunction& rho,
const DiscreteVectorFunction& u,
const DiscreteScalarFunction& E,
const DiscreteTensorFunction& CG,
const DiscreteScalarFunction& rho_app,
const DiscreteTensorFunction& CG_app,
const NodeValue<const Rd>& ur,
const NodeValuePerCell<const Rd>& Fjr) const
{
const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
if ((mesh.shared_connectivity() != ur.connectivity_ptr()) or
(mesh.shared_connectivity() != Fjr.connectivity_ptr())) {
throw NormalError("fluxes are not defined on the same connectivity than the mesh");
}
NodeValue<Rd> new_xr = copy(mesh.xr());
parallel_for(
mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });
std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr);
CellValue<const double> Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
const NodeValuePerCell<const Rd> Cjr = MeshDataManager::instance().getMeshData(mesh).Cjr();
CellValue<double> new_rho = copy(rho.cellValues());
CellValue<Rd> new_u = copy(u.cellValues());
CellValue<double> new_E = copy(E.cellValues());
CellValue<Rdxd> new_CG = copy(CG.cellValues());
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
const auto& cell_nodes = cell_to_node_matrix[j];
Rd momentum_fluxes = zero;
double energy_fluxes = 0;
Rdxd gradv = zero;
for (size_t R = 0; R < cell_nodes.size(); ++R) {
const NodeId r = cell_nodes[R];
gradv += tensorProduct(ur[r], Cjr(j, R)); momentum_fluxes += Fjr(j, R);
energy_fluxes += dot(Fjr(j, R), ur[r]);
}
const Rdxd cauchy_green_fluxes = rho_app[j] * (gradv * CG_app[j] + CG_app[j] * transpose(gradv));
const double dt_over_Mj = dt / (rho[j] * Vj[j]);
new_u[j] += dt_over_Mj * momentum_fluxes;
new_E[j] += dt_over_Mj * energy_fluxes;
new_CG[j] += dt_over_Mj * cauchy_green_fluxes;
new_CG[j] += transpose(new_CG[j]);
new_CG[j] *= 0.5;
});
CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { new_rho[j] *= Vj[j] / new_Vj[j]; });
return {std::make_shared<MeshVariant>(new_mesh),
std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_rho)),
std::make_shared<DiscreteFunctionVariant>(DiscreteVectorFunction(new_mesh, new_u)),
std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_E)),
std::make_shared<DiscreteFunctionVariant>(DiscreteTensorFunction(new_mesh, new_CG))};
}
std::tuple<std::shared_ptr<const MeshVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
order2_apply_fluxes(const double& dt,
const std::shared_ptr<const DiscreteFunctionVariant>& rho,
const std::shared_ptr<const DiscreteFunctionVariant>& u,
const std::shared_ptr<const DiscreteFunctionVariant>& E,
const std::shared_ptr<const DiscreteFunctionVariant>& CG,
const std::shared_ptr<const DiscreteFunctionVariant>& rho_app,
const std::shared_ptr<const DiscreteFunctionVariant>& CG_app,
const std::shared_ptr<const ItemValueVariant>& ur,
const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr) const
{
std::shared_ptr mesh_v = getCommonMesh({rho, u, E});
if (not mesh_v) {
throw NormalError("discrete functions are not defined on the same mesh");
}
if (not checkDiscretizationType({rho, u, E}, DiscreteFunctionType::P0)) {
throw NormalError("hyperelastic solver expects P0 functions");
}
return this->order2_apply_fluxes(dt, //
*mesh_v->get<MeshType>(), //
rho->get<DiscreteScalarFunction>(), //
u->get<DiscreteVectorFunction>(), //
E->get<DiscreteScalarFunction>(), //
CG->get<DiscreteTensorFunction>(), //
rho_app->get<DiscreteScalarFunction>(), //
CG_app->get<DiscreteTensorFunction>(), //
ur->get<NodeValue<const Rd>>(), //
Fjr->get<NodeValuePerCell<const Rd>>());
}
std::tuple<std::shared_ptr<const MeshVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
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::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, bc_descriptor_list);
// auto [m_app, rho_app, u_app, E_app, CG_app] = apply_fluxes(dt/2., rho, u, E, CG, ur, Fjr);
// std::shared_ptr<const DiscreteFunctionVariant> chi_solid_app = shallowCopy(m_app, chi_solid);
// auto [sigma_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,bc_descriptor_list);
return apply_fluxes(dt, rho, u, E, CG, ur, Fjr);
}
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());
}