Adds operators for evaluating stress and strain based on the material model: i.e. Visit the variant to handle either HookeOps::CommonData or HenckyOps::CommonData. i.e.
   36 
   39 
   40    using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
   41 
   42    struct OpScale : public ForcesAndSourcesCore::UserDataOperator {
   43      OpScale(boost::shared_ptr<MatrixDouble> m_ptr, double s)
   45            mPtr(m_ptr), 
scale(s) {}
 
   46      MoFEMErrorCode doWork(
int, EntityType, EntitiesFieldData::EntData &) {
 
   48        return 0;
   49      }
   50 
   51    private:
   52      boost::shared_ptr<MatrixDouble> mPtr;
   54    };
   55 
   56    auto push_domain_ops = [&](auto &pip) {
   58                            pip, {
H1, 
HDIV}, 
"GEOMETRY")),
 
   59                        "Apply base transform");
   60      
   61      
   62      using CommonDataVariant =
   63          std::variant<boost::shared_ptr<HookeOps::CommonData>,
   64                       boost::shared_ptr<HenckyOps::CommonData>>;
   65      auto contact_stress_ptr = boost::make_shared<MatrixDouble>();
   66      pip.push_back(new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
   67          "SIGMA", contact_stress_ptr));
   68 
   70        auto hooke_common_data_ptr =
   71            HookeOps::commonDataFactory<SPACE_DIM, GAUSS, DomainEleOp>(
   72                *m_field_ptr, pip, 
"U", 
"MAT_ELASTIC", Sev::inform, 
scale);
 
   73 
   74        pip.push_back(
new OpScale(contact_stress_ptr, 
scale));
 
   75        return std::make_tuple(CommonDataVariant(hooke_common_data_ptr),
   76                               contact_stress_ptr);
   77      } else { 
   78        auto hencky_common_data_ptr =
   79            HenckyOps::commonDataFactory<SPACE_DIM, GAUSS, DomainEleOp>(
   80                *m_field_ptr, pip, 
"U", 
"MAT_ELASTIC", Sev::inform, 
scale);
 
   81        pip.push_back(
new OpScale(contact_stress_ptr, 
scale));
 
   82        return std::make_tuple(CommonDataVariant(hencky_common_data_ptr),
   83                               contact_stress_ptr);
   84      }
   85    };
   86    auto push_bdy_ops = [&](auto &pip) {
   87      
   88      auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
   89      pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
   90          "U", common_data_ptr->contactDispPtr()));
   91      pip.push_back(new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
   92          "SIGMA", common_data_ptr->contactTractionPtr()));
   93      pip.push_back(
new OpScale(common_data_ptr->contactTractionPtr(), 
scale));
 
   94      using C = ContactIntegrators<BoundaryEleOp>;
   95      pip.push_back(new typename C::template OpEvaluateSDF<SPACE_DIM, GAUSS>(
   96          common_data_ptr));
   97      return common_data_ptr;
   98    };
   99 
  100    auto get_domain_pip = [&](auto &pip)
  101        -> boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> & {
  103        auto op_loop_side = new OpLoopSide<SideEle>(
  104            *m_field_ptr, 
"dFE", 
SPACE_DIM, Sev::noisy,
 
  105            boost::make_shared<
  106                ForcesAndSourcesCore::UserDataOperator::AdjCache>());
  107        pip.push_back(op_loop_side);
  108        return op_loop_side->getOpPtrVector();
  109      } else {
  110        return pip;
  111      }
  112    };
  113 
  114    auto get_post_proc_domain_fe = [&]() {
  115      auto post_proc_fe =
  116          boost::make_shared<PostProcEleDomain>(*m_field_ptr, 
postProcMesh);
 
  117      auto &pip = post_proc_fe->getOpPtrVector();
  118      
  119      auto [hencky_or_hooke_common_data_ptr, contact_stress_ptr] =
  120          push_domain_ops(get_domain_pip(pip));
  121 
  122      auto u_ptr = boost::make_shared<MatrixDouble>();
  123      pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
  124      auto X_ptr = boost::make_shared<MatrixDouble>();
  125      pip.push_back(
  126          new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", X_ptr));
  127
  128
  129
  130
  131
  132
  133
  134
  135
  136
  137      std::visit(
  138          [&](auto &common_data_ptr) {
  139            if constexpr (std::is_same_v<
  140                              std::decay_t<decltype(common_data_ptr)>,
  141                              boost::shared_ptr<HookeOps::CommonData>>) {
  143                  post_proc_fe->getPostProcMesh(),
  144                  post_proc_fe->getMapGaussPts(), {},
  145                  {{"U", u_ptr}, {"GEOMETRY", X_ptr}},
  146                  {{"SIGMA", contact_stress_ptr},
  147                   {"G", common_data_ptr->matGradPtr}}, 
  148                  {{"STRESS", common_data_ptr->getMatCauchyStress()},
  149                   {"STRAIN", common_data_ptr->getMatStrain()}}));
  150            }
  151            
  152            else if constexpr (std::is_same_v<
  153                                   std::decay_t<decltype(common_data_ptr)>,
  154                                   boost::shared_ptr<HenckyOps::CommonData>>) {
  155 
  157                  post_proc_fe->getPostProcMesh(),
  158                  post_proc_fe->getMapGaussPts(), {},
  159                  {{"U", u_ptr}, {"GEOMETRY", X_ptr}},
  160                  {{"SIGMA", contact_stress_ptr},
  161                   {"G", common_data_ptr->matGradPtr},
  162                   {"PK1", common_data_ptr->getMatFirstPiolaStress()}
  163 
  164                  },
  165                  {{"HENCKY_STRAIN", common_data_ptr->getMatLogC()}}));
  166            }
  167          },
  168          hencky_or_hooke_common_data_ptr);
  169 
  171 
  173                              pip, {
HDIV}, 
"GEOMETRY")),
 
  174                          "Apply transform");
  175        auto common_data_ptr = push_bdy_ops(pip);
  176 
  177        pip.push_back(
  178 
  180 
  181                post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
  182 
  183                {{"SDF", common_data_ptr->sdfPtr()},
  184                 {"CONSTRAINT_CONTACT", common_data_ptr->constraintPtr()}},
  185 
  186                {
  187 
  188                    {"TRACTION_CONTACT", common_data_ptr->contactTractionPtr()},
  189                    {"GRAD_SDF", common_data_ptr->gradSdfPtr()}
  190 
  191                },
  192 
  193                {},
  194 
  195                {{"HESS_SDF", common_data_ptr->hessSdfPtr()}}
  196 
  197                )
  198 
  199        );
  200      }
  201 
  202      return post_proc_fe;
  203    };
  204 
  205    auto get_post_proc_bdy_fe = [&]() {
  206      auto post_proc_fe =
  207          boost::make_shared<PostProcEleBdy>(*m_field_ptr, 
postProcMesh);
 
  208      auto &pip = post_proc_fe->getOpPtrVector();
  209 
  211                            pip, {
HDIV}, 
"GEOMETRY")),
 
  212                        "Apply transform");
  213      auto common_data_ptr = push_bdy_ops(pip);
  214 
  215      
  216      auto op_loop_side = new OpLoopSide<SideEle>(
  217          *m_field_ptr, 
"dFE", 
SPACE_DIM, Sev::noisy,
 
  218          boost::make_shared<
  219              ForcesAndSourcesCore::UserDataOperator::AdjCache>());
  220      pip.push_back(op_loop_side);
  221 
  222      auto [hencky_or_hooke_common_data_ptr, contact_stress_ptr] =
  223          push_domain_ops(op_loop_side->getOpPtrVector());
  224 
  225      auto X_ptr = boost::make_shared<MatrixDouble>();
  226      pip.push_back(
  227          new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", X_ptr));
  228 
  229      pip.push_back(
  230 
  232 
  233              post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
  234 
  235              {{"SDF", common_data_ptr->sdfPtr()},
  236               {"CONSTRAINT_CONTACT", common_data_ptr->constraintPtr()}},
  237 
  238              {{"U", common_data_ptr->contactDispPtr()},
  239               {"GEOMETRY", X_ptr},
  240               {"TRACTION_CONTACT", common_data_ptr->contactTractionPtr()},
  241               {"GRAD_SDF", common_data_ptr->gradSdfPtr()}
  242 
  243              },
  244 
  245              {},
  246 
  247              {{"HESS_SDF", common_data_ptr->hessSdfPtr()}}
  248 
  249              )
  250 
  251      );
  252 
  253      return post_proc_fe;
  254    };
  255 
  256    auto get_integrate_traction = [&]() {
  257      auto integrate_traction = boost::make_shared<BoundaryEle>(*m_field_ptr);
  260              integrate_traction->getOpPtrVector(), {HDIV}, "GEOMETRY")),
  261          "Apply transform");
  262      
  263      
  264      integrate_traction->getOpPtrVector().push_back(
  265          new OpSetHOWeightsOnSubDim<SPACE_DIM>());
  268      };
  269 
  271          (opFactoryCalculateTraction<SPACE_DIM, GAUSS, BoundaryEleOp>(
  273          "push operators to calculate traction");
  274 
  275      return integrate_traction;
  276    };
  277 
  278    auto get_integrate_area = [&]() {
  279      auto integrate_area = boost::make_shared<BoundaryEle>(*m_field_ptr);
  280 
  283              integrate_area->getOpPtrVector(), {HDIV}, "GEOMETRY")),
  284          "Apply transform");
  285      
  286      
  287      integrate_area->getOpPtrVector().push_back(
  288          new OpSetHOWeightsOnSubDim<SPACE_DIM>());
  291      };
  294           m_field_ptr->getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
  295               std::regex((boost::format(
"%s(.*)") % 
"CONTACT").str()))) {
 
  296        auto meshset = 
m->getMeshset();
 
  297        Range contact_meshset_range;
 
  299            meshset, 
SPACE_DIM - 1, contact_meshset_range, 
true);
 
  300 
  302            contact_meshset_range);
  303        contact_range.merge(contact_meshset_range);
  304      }
  305 
  306      auto contact_range_ptr = boost::make_shared<Range>(contact_range);
  307 
  308      auto op_loop_side = new OpLoopSide<SideEle>(
  309          *m_field_ptr, m_field_ptr->
getInterface<Simple>()->getDomainFEName(),
 
  312          op_loop_side->getOpPtrVector(), {H1}, "GEOMETRY");
  313 
  315          (opFactoryCalculateArea<SPACE_DIM, GAUSS, BoundaryEleOp>(
  316              integrate_area->getOpPtrVector(), op_loop_side, "SIGMA", "U",
  318          "push operators to calculate area");
  319 
  320      return integrate_area;
  321    };
  322 
  326 
  329 
  333  }
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
@ HDIV
field with continuous normal traction
#define CHKERR
Inline error check.
PetscErrorCode DMoFEMGetInterfacePtr(DM dm, MoFEM::Interface **m_field_ptr)
Get pointer to MoFEM::Interface.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
int geom_order
Order if fixed.
static constexpr int approx_order
FTensor::Index< 'm', 3 > m
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
Deprecated interface functions.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.