diff --git a/src/scheme/HypoelasticSolver.cpp b/src/scheme/HypoelasticSolver.cpp new file mode 100644 index 0000000000000000000000000000000000000000..adf548b2448f3c7f478341df704dd84c4d95b0a3 --- /dev/null +++ b/src/scheme/HypoelasticSolver.cpp @@ -0,0 +1,982 @@ +#include <scheme/HyperelasticSolver.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/SubItemValuePerItemVariant.hpp> +#include <scheme/DirichletBoundaryConditionDescriptor.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/SymmetryBoundaryConditionDescriptor.hpp> +#include <utils/Socket.hpp> + +#include <variant> +#include <vector> + +template <size_t Dimension> +double +hypoelastic_dt(const DiscreteFunctionP0<Dimension, const double>& c) +{ + const Mesh<Connectivity<Dimension>>& mesh = dynamic_cast<const Mesh<Connectivity<Dimension>>&>(*c.mesh()); + + const auto Vj = MeshDataManager::instance().getMeshData(mesh).Vj(); + const auto Sj = MeshDataManager::instance().getMeshData(mesh).sumOverRLjr(); + + CellValue<double> local_dt{mesh.connectivity()}; + parallel_for( + mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { local_dt[j] = 2 * Vj[j] / (Sj[j] * c[j]); }); + + return min(local_dt); +} + +double +hypoelastic_dt(const std::shared_ptr<const DiscreteFunctionVariant>& c) +{ + std::shared_ptr mesh = getCommonMesh({c}); + + switch (mesh->dimension()) { + case 1: { + return hypoelastic_dt(c->get<DiscreteFunctionP0<1, const double>>()); + } + case 2: { + return hypoelastic_dt(c->get<DiscreteFunctionP0<2, const double>>()); + } + case 3: { + return hypoelastic_dt(c->get<DiscreteFunctionP0<3, const double>>()); + } + default: { + throw UnexpectedError("invalid mesh dimension"); + } + } +} + +template <size_t Dimension> +class HypoelasticSolverHandler::HypoelasticSolver final : public HypoelasticSolverHandler::IHypoelasticSolver +{ + private: + using Rdxd = TinyMatrix<Dimension>; + using Rd = TinyVector<Dimension>; + + using MeshType = Mesh<Connectivity<Dimension>>; + using MeshDataType = MeshData<Dimension>; + + using DiscreteScalarFunction = DiscreteFunctionP0<Dimension, const double>; + using DiscreteVectorFunction = DiscreteFunctionP0<Dimension, const Rd>; + using DiscreteTensorFunction = DiscreteFunctionP0<Dimension, 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<Connectivity<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; + } + + 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<Dimension> 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<Dimension> 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<Dimension> 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<Dimension> 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<Dimension> 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 hypoelastic 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, + 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; + } + + 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 i_mesh = getCommonMesh({rho_v, aL_v, aT_v, u_v, sigma_v}); + if (not i_mesh) { + 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("hypoelastic solver expects P0 functions"); + } + + const MeshType& mesh = dynamic_cast<const MeshType&>(*i_mesh); + 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<std::shared_ptr<const IMesh>, + std::shared_ptr<const DiscreteFunctionVariant>, + std::shared_ptr<const DiscreteFunctionVariant>, + std::shared_ptr<const DiscreteFunctionVariant>, + std::shared_ptr<const DiscreteFunctionVariant>> + apply_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 {new_mesh, std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_rho)), + std::make_shared<DiscreteFunctionVariant>(DiscreteVectorFunction(new_mesh, new_u)), + std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_E)), + std::make_shared<DiscreteFunctionVariant>(DiscreteTensorFunction(new_mesh, new_CG))}; + } + + std::tuple<std::shared_ptr<const IMesh>, + std::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 i_mesh = getCommonMesh({rho, u, E}); + if (not i_mesh) { + throw NormalError("discrete functions are not defined on the same mesh"); + } + + if (not checkDiscretizationType({rho, u, E}, DiscreteFunctionType::P0)) { + throw NormalError("hypoelastic solver expects P0 functions"); + } + + return this->apply_fluxes(dt, // + dynamic_cast<const MeshType&>(*i_mesh), // + rho->get<DiscreteScalarFunction>(), // + u->get<DiscreteVectorFunction>(), // + E->get<DiscreteScalarFunction>(), // + CG->get<DiscreteTensorFunction>(), // + ur->get<NodeValue<const Rd>>(), // + Fjr->get<NodeValuePerCell<const Rd>>()); + } + + std::tuple<std::shared_ptr<const IMesh>, + std::shared_ptr<const DiscreteFunctionVariant>, + std::shared_ptr<const DiscreteFunctionVariant>, + std::shared_ptr<const DiscreteFunctionVariant>, + std::shared_ptr<const DiscreteFunctionVariant>> + apply(const SolverType& solver_type, + 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 + { + auto [ur, Fjr] = compute_fluxes(solver_type, rho, aL, aT, u, sigma, bc_descriptor_list); + return apply_fluxes(dt, rho, u, E, CG, ur, Fjr); + } + + HypoelasticSolver() = default; + HypoelasticSolver(HypoelasticSolver&&) = default; + ~HypoelasticSolver() = default; +}; + +template <size_t Dimension> +void +HypoelasticSolverHandler::HypoelasticSolver<Dimension>::_applyPressureBC(const BoundaryConditionList& bc_list, + const MeshType& mesh, + NodeValue<Rd>& br) const +{ + for (const auto& boundary_condition : bc_list) { + std::visit( + [&](auto&& bc) { + using T = std::decay_t<decltype(bc)>; + if constexpr (std::is_same_v<PressureBoundaryCondition, T>) { + MeshData<Dimension>& mesh_data = MeshDataManager::instance().getMeshData(mesh); + 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 <size_t Dimension> +void +HypoelasticSolverHandler::HypoelasticSolver<Dimension>::_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>) { + MeshData<Dimension>& 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 <size_t Dimension> +void +HypoelasticSolverHandler::HypoelasticSolver<Dimension>::_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 <size_t Dimension> +void +HypoelasticSolverHandler::HypoelasticSolver<Dimension>::_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 <size_t Dimension> +class HypoelasticSolverHandler::HypoelasticSolver<Dimension>::FixedBoundaryCondition +{ + private: + const MeshNodeBoundary<Dimension> m_mesh_node_boundary; + + public: + const Array<const NodeId>& + nodeList() const + { + return m_mesh_node_boundary.nodeList(); + } + + FixedBoundaryCondition(const MeshNodeBoundary<Dimension> mesh_node_boundary) + : m_mesh_node_boundary{mesh_node_boundary} + {} + + ~FixedBoundaryCondition() = default; +}; + +template <size_t Dimension> +class HypoelasticSolverHandler::HypoelasticSolver<Dimension>::PressureBoundaryCondition +{ + private: + const MeshFaceBoundary<Dimension> 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<Dimension>& 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 HypoelasticSolverHandler::HypoelasticSolver<1>::PressureBoundaryCondition +{ + private: + const MeshNodeBoundary<1> 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<1>& mesh_node_boundary, const Array<const double>& value_list) + : m_mesh_node_boundary{mesh_node_boundary}, m_value_list{value_list} + {} + + ~PressureBoundaryCondition() = default; +}; + +template <size_t Dimension> +class HypoelasticSolverHandler::HypoelasticSolver<Dimension>::NormalStressBoundaryCondition +{ + private: + const MeshFaceBoundary<Dimension> 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<Dimension>& 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 HypoelasticSolverHandler::HypoelasticSolver<1>::NormalStressBoundaryCondition +{ + private: + const MeshNodeBoundary<1> 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<1>& mesh_node_boundary, const Array<const Rdxd>& value_list) + : m_mesh_node_boundary{mesh_node_boundary}, m_value_list{value_list} + {} + + ~NormalStressBoundaryCondition() = default; +}; + +template <size_t Dimension> +class HypoelasticSolverHandler::HypoelasticSolver<Dimension>::VelocityBoundaryCondition +{ + private: + const MeshNodeBoundary<Dimension> 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<Dimension>& 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 <size_t Dimension> +class HypoelasticSolverHandler::HypoelasticSolver<Dimension>::SymmetryBoundaryCondition +{ + public: + using Rd = TinyVector<Dimension, double>; + + private: + const MeshFlatNodeBoundary<Dimension> 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<Dimension>& mesh_flat_node_boundary) + : m_mesh_flat_node_boundary(mesh_flat_node_boundary) + { + ; + } + + ~SymmetryBoundaryCondition() = default; +}; + +HypoelasticSolverHandler::HypoelasticSolverHandler(const std::shared_ptr<const IMesh>& i_mesh) +{ + if (not i_mesh) { + throw NormalError("discrete functions are not defined on the same mesh"); + } + + switch (i_mesh->dimension()) { + case 1: { + m_hypoelastic_solver = std::make_unique<HypoelasticSolver<1>>(); + break; + } + case 2: { + m_hypoelastic_solver = std::make_unique<HypoelasticSolver<2>>(); + break; + } + case 3: { + m_hypoelastic_solver = std::make_unique<HypoelasticSolver<3>>(); + break; + } + default: { + throw UnexpectedError("invalid mesh dimension"); + } + } +} diff --git a/src/scheme/HypoelasticSolver.hpp b/src/scheme/HypoelasticSolver.hpp new file mode 100644 index 0000000000000000000000000000000000000000..8664f47ebe173ec04c40b05763ce5057e17d4605 --- /dev/null +++ b/src/scheme/HypoelasticSolver.hpp @@ -0,0 +1,90 @@ +#ifndef HYPOELASTIC_SOLVER_HPP +#define HYPOELASTIC_SOLVER_HPP + +#include <memory> +#include <tuple> +#include <vector> + +class IBoundaryConditionDescriptor; +class IMesh; +class ItemValueVariant; +class SubItemValuePerItemVariant; +class DiscreteFunctionVariant; + +double hypoelastic_dt(const std::shared_ptr<const DiscreteFunctionVariant>& c); + +class HypoelasticSolverHandler +{ + public: + enum class SolverType + { + Glace, + Eucclhyd + }; + + private: + struct IHypoelasticSolver + { + virtual 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, + const std::shared_ptr<const DiscreteFunctionVariant>& aL, + const std::shared_ptr<const DiscreteFunctionVariant>& aT, + const std::shared_ptr<const DiscreteFunctionVariant>& u, + const std::shared_ptr<const DiscreteFunctionVariant>& sigma, + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const = 0; + + virtual std::tuple<std::shared_ptr<const IMesh>, + 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 = 0; + + virtual std::tuple<std::shared_ptr<const IMesh>, + std::shared_ptr<const DiscreteFunctionVariant>, + std::shared_ptr<const DiscreteFunctionVariant>, + std::shared_ptr<const DiscreteFunctionVariant>, + std::shared_ptr<const DiscreteFunctionVariant>> + apply(const SolverType& solver_type, + 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>& p, + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const = 0; + + IHypoelasticSolver() = default; + IHypoelasticSolver(IHypoelasticSolver&&) = default; + IHypoelasticSolver& operator=(IHypoelasticSolver&&) = default; + + virtual ~IHypoelasticSolver() = default; + }; + + template <size_t Dimension> + class HypoelasticSolver; + + std::unique_ptr<IHypoelasticSolver> m_hypoelastic_solver; + + public: + const IHypoelasticSolver& + solver() const + { + return *m_hypoelastic_solver; + } + + HypoelasticSolverHandler(const std::shared_ptr<const IMesh>& mesh); +}; + +#endif // HYPOELASTIC_SOLVER_HPP