9#ifndef __HOOKE_OPS_HPP__ 
   10#define __HOOKE_OPS_HPP__ 
   14struct CommonData : 
public boost::enable_shared_from_this<CommonData> {
 
   25    return boost::shared_ptr<MatrixDouble>(shared_from_this(),
 
 
   33    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &
matStrainPtr);
 
 
 
   54               boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
 
   55               std::string block_name,
 
   56               boost::shared_ptr<MatrixDouble> mat_D_Ptr, Sev sev,
 
   64                std::vector<const CubitMeshSets *> meshset_vec_ptr,
 
   69          scaleYoungModulus(
scale) {
 
   71                        "Can not get data from block");
 
   74    MoFEMErrorCode doWork(
int side, EntityType type,
 
   75                          EntitiesFieldData::EntData &data) {
 
   78      for (
auto &b : blockData) {
 
   80        if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
 
   81          CHKERR getMatDPtr(matDPtr, b.bulkModulusK * scaleYoungModulus,
 
   82                            b.shearModulusG * scaleYoungModulus, isPlaneStrain);
 
   87      CHKERR getMatDPtr(matDPtr, bulkModulusKDefault * scaleYoungModulus,
 
   88                        shearModulusGDefault * scaleYoungModulus,
 
   94    boost::shared_ptr<MatrixDouble> matDPtr;
 
   95    const double scaleYoungModulus;
 
  103    double bulkModulusKDefault;
 
  104    double shearModulusGDefault;
 
  105    std::vector<BlockData> blockData;
 
  110                     std::vector<const CubitMeshSets *> meshset_vec_ptr,
 
  114      for (
auto m : meshset_vec_ptr) {
 
  116        std::vector<double> block_data;
 
  117        CHKERR m->getAttributes(block_data);
 
  118        if (block_data.size() < 2) {
 
  120                  "Expected that block has two attribute");
 
  122        auto get_block_ents = [&]() {
 
  125          m_field.
get_moab().get_entities_by_handle(
m->meshset, ents, 
true);
 
  145    MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
 
  150      auto set_material_stiffness = [&]() {
 
  169        auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*mat_D_ptr);
 
  176      constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
 
  178      set_material_stiffness();
 
  187  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, 
"", 
"-young_modulus", &
E,
 
  189  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, 
"", 
"-poisson_ratio", &nu,
 
  193  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, 
"", 
"-plane_strain",
 
  198  pip.push_back(
new OpMatBlocks(
 
  202      m_field.
getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
 
  204          (boost::format(
"%s(.*)") % block_name).str()
 
 
  214template <
int DIM, IntegrationType I, 
typename DomainEleOp>
 
  217    boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
 
  218    std::string 
field_name, std::string block_name, Sev sev, 
double scale = 1) {
 
  220  auto common_ptr = boost::make_shared<HookeOps::CommonData>();
 
  221  common_ptr->matDPtr = boost::make_shared<MatrixDouble>();
 
  222  common_ptr->matGradPtr = boost::make_shared<MatrixDouble>();
 
  225                                        common_ptr->matDPtr, sev, 
scale),
 
  228  pip.push_back(
new OpCalculateVectorFieldGradient<DIM, DIM>(
 
  230  pip.push_back(
new OpSymmetrizeTensor<SPACE_DIM>(common_ptr->matGradPtr,
 
  231                                                  common_ptr->getMatStrain()));
 
  232  pip.push_back(
new OpTensorTimesSymmetricTensor<SPACE_DIM, SPACE_DIM>(
 
  233      common_ptr->getMatStrain(), common_ptr->getMatCauchyStress(),
 
  234      common_ptr->matDPtr));
 
 
  248template <
int DIM, AssemblyType A, IntegrationType I, 
typename DomainEleOp>
 
  251    boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
 
  252    std::string 
field_name, boost::shared_ptr<HookeOps::CommonData> common_ptr,
 
  253    Sev sev, 
bool is_non_linear = 
false) {
 
  256  using B = 
typename FormsIntegrators<DomainEleOp>::template Assembly<
 
  258  using OpInternalForce =
 
  259      typename B::template OpGradTimesSymTensor<1, DIM, DIM>;
 
  262      new OpInternalForce(
field_name, common_ptr->getMatCauchyStress(),
 
  263                          [is_non_linear](
double, 
double, 
double) 
constexpr {
 
  264                            return (is_non_linear ? 1 : -1);
 
 
  287template <
int DIM, AssemblyType A, IntegrationType I, 
typename DomainEleOp>
 
  290    boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
 
  291    std::string 
field_name, std::string block_name, Sev sev,
 
  292    bool is_non_linear = 
false, 
double scale = 1) {
 
  295  auto common_ptr = commonDataFactory<DIM, I, DomainEleOp>(
 
  297  CHKERR opFactoryDomainRhs<DIM, A, I, DomainEleOp>(
 
  298      m_field, pip, 
field_name, common_ptr, sev, is_non_linear);
 
 
  308template <
int DIM, AssemblyType A, IntegrationType I, 
typename DomainEleOp>
 
  311    boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
 
  312    std::string 
field_name, boost::shared_ptr<HookeOps::CommonData> common_ptr,
 
  316  using B = 
typename FormsIntegrators<DomainEleOp>::template Assembly<
 
  319  using OpKCauchy = 
typename B::template OpGradSymTensorGrad<1, DIM, DIM, 0>;
 
 
  332template <
int DIM, AssemblyType A, IntegrationType I, 
typename DomainEleOp>
 
  335    boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
 
  336    std::string 
field_name, std::string block_name, Sev sev, 
double scale = 1) {
 
  339  auto common_ptr = commonDataFactory<DIM, I, DomainEleOp>(
 
 
 
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
PetscBool is_plane_strain
DomainEle::UserDataOperator DomainEleOp
Finire element operator type.
Kronecker Delta class symmetric.
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
MoFEMErrorCode addMatBlockOps(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string block_name, boost::shared_ptr< MatrixDouble > mat_D_Ptr, Sev sev, double scale=1)
// Add operator to compute material stiffness tensor D based on block attributes.
MoFEMErrorCode opFactoryDomainLhs(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, boost::shared_ptr< HookeOps::CommonData > common_ptr, Sev sev)
Assemble domain LHS K factory (LHS first overload) Initializes and pushes operators to assemble the L...
auto commonDataFactory(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, std::string block_name, Sev sev, double scale=1)
MoFEMErrorCode opFactoryDomainRhs(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, boost::shared_ptr< HookeOps::CommonData > common_ptr, Sev sev, bool is_non_linear=false)
Factory function to create and push internal force operators for the domain RHS. (RHS first overload)...
double poisson_ratio
Poisson ratio.
double young_modulus
Young modulus.
constexpr auto field_name
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Only used with Hooke equation (linear material model)]
FTensor::Index< 'm', 3 > m
boost::shared_ptr< MatrixDouble > matDPtr
auto getMatStrain()
Get shared strain.
auto getMatCauchyStress()
Get shared Cauchy stress.
MatrixDouble matStrainPtr
boost::shared_ptr< MatrixDouble > matGradPtr
MatrixDouble matCauchyStressPtr
virtual moab::Interface & get_moab()=0
Deprecated interface functions.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.