Order2HyperelasticSolver.cpp

#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;