PostProcHookStresses.hpp

Go to the documentation of this file.

/**
 * \file PostProcHookStress.hpp
 * \brief Post-proc stresses for linear Hooke isotropic material
 *
 * \ingroup nonlinear_elastic_elem
 */
 
 
 
/**
 * \brief Operator post-procesing stresses for Hook isotropic material
 
 * Example how to use it
 
 \code
 PostProcBrokenMeshInMoab<VolumeElementForcesAndSourcesCore> post_proc(m_field);
 {
   auto disp_ptr = boost::make_shared<MatrixDouble>();
   auto disp_grad_ptr = boost::make_shared<MatrixDouble>();
   post_proc.getOpPtrVector().push_back(
       new OpCalculateVectorFieldValues<3>("DISPLACEMENT", disp_ptr));
   post_proc.getOpPtrVector().push_back(
       new OpCalculateVectorFieldGradient<3, 3>(
           "DISPLACEMENT", disp_grad_ptr));
   //add postprocessing for stresses
   post_proc.getOpPtrVector().push_back(
     new PostProcHookStress(
       m_field,
       post_proc.getPostProcMesh(),
       post_proc.getMapGaussPts(),
       post_proc.getPostProcElements(),
       "DISPLACEMENT",
       disp_grad_ptr,
       &elastic.setOfBlocks
     )
   );
   using OpPPMap = OpPostProcMapInMoab<3, 3>;
   post_proc.getOpPtrVector().push_back(
       new OpPPMap(post_proc.getPostProcMesh(), post_proc.getMapGaussPts(),
                   {}, {{"DISPLACEMENT", disp_ptr}},
                   {{"DISPLACEMENT_GRAD", disp_grad_ptr}}, {}));
   CHKERR DMoFEMLoopFiniteElements(dm,"ELASTIC",&post_proc);
   CHKERR post_proc.writeFile("out.h5m");
 }
 
 \endcode
 
 */
struct PostProcHookStress
    : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
 
  MoFEM::Interface &mField;
  moab::Interface &postProcMesh;
  std::vector<EntityHandle> &mapGaussPts;
  Range &postProcElements;
  bool isFieldDisp;
 
#ifdef __NONLINEAR_ELASTIC_HPP
  /// Material block data, ket is block id
  const std::map<int, NonlinearElasticElement::BlockData>
      *setOfBlocksMaterialDataPtr;
#endif //__NONLINEAR_ELASTIC_HPP
 
  boost::shared_ptr<MatrixDouble> fieldGradientPtr;
 
  /**
   * Constructor
   */
  PostProcHookStress(MoFEM::Interface &m_field, moab::Interface &post_proc_mesh,
                     std::vector<EntityHandle> &map_gauss_pts,
                     Range &post_proc_elements,
                     const std::string field_name,
                     boost::shared_ptr<MatrixDouble> field_gradient_ptr,
#ifdef __NONLINEAR_ELASTIC_HPP
                     const std::map<int, NonlinearElasticElement::BlockData>
                         *set_of_block_data_ptr = NULL,
#endif
                     const bool is_field_disp = true)
      : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
            field_name, ForcesAndSourcesCore::UserDataOperator::OPROW),
        mField(m_field), postProcMesh(post_proc_mesh),
        mapGaussPts(map_gauss_pts), postProcElements(post_proc_elements),
#ifdef __NONLINEAR_ELASTIC_HPP
        setOfBlocksMaterialDataPtr(set_of_block_data_ptr),
#endif //__NONLINEAR_ELASTIC_HPP
        fieldGradientPtr(field_gradient_ptr), isFieldDisp(is_field_disp) {
  }
 
  /**
   * \brief get material parameter
 
   * Material parameters are read form BlockSet, however if block data are
   present,
   * use data how are set for elastic element operators.
 
   * @param  _lambda   elastic material constant
   * @param  _mu       elastic material constant
   * @param  _block_id  block id
   * @return           error code
 
   */
  MoFEMErrorCode getMatParameters(double *_lambda, double *_mu,
                                  int *_block_id) {
    MoFEMFunctionBegin;
 
    *_lambda = 1;
    *_mu = 1;
 
    EntityHandle ent = getNumeredEntFiniteElementPtr()->getEnt();
    for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
             mField, BLOCKSET | MAT_ELASTICSET, it)) {
      Mat_Elastic mydata;
      CHKERR it->getAttributeDataStructure(mydata);
 
      Range meshsets;
      CHKERR mField.get_moab().get_entities_by_type(it->meshset, MBENTITYSET,
                                                    meshsets, false);
      meshsets.insert(it->meshset);
      for (Range::iterator mit = meshsets.begin(); mit != meshsets.end();
           mit++) {
        if (mField.get_moab().contains_entities(*mit, &ent, 1)) {
          *_lambda = LAMBDA(mydata.data.Young, mydata.data.Poisson);
          *_mu = MU(mydata.data.Young, mydata.data.Poisson);
          *_block_id = it->getMeshsetId();
#ifdef __NONLINEAR_ELASTIC_HPP
          if (setOfBlocksMaterialDataPtr) {
            *_lambda =
                LAMBDA(setOfBlocksMaterialDataPtr->at(*_block_id).E,
                       setOfBlocksMaterialDataPtr->at(*_block_id).PoissonRatio);
            *_mu = MU(setOfBlocksMaterialDataPtr->at(*_block_id).E,
                      setOfBlocksMaterialDataPtr->at(*_block_id).PoissonRatio);
          }
#endif //__NONLINEAR_ELASTIC_HPP
          MoFEMFunctionReturnHot(0);
        }
      }
    }
 
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Element is not in elastic block, however you run linear elastic "
            "analysis with that element\n"
            "top tip: check if you update block sets after mesh refinements or "
            "interface insertion");
 
    MoFEMFunctionReturn(0);
  }
 
  /**
   * \brief Here real work is done
   */
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data) {
    MoFEMFunctionBegin;
 
    if (type != MBVERTEX)
      MoFEMFunctionReturnHot(0);
    if (data.getFieldData().size() == 0)
      MoFEMFunctionReturnHot(0);
 
    int id;
    double lambda, mu;
    CHKERR getMatParameters(&lambda, &mu, &id);
 
    MatrixDouble D_lambda, D_mu, D;
    D_lambda.resize(6, 6);
    D_lambda.clear();
    for (int rr = 0; rr < 3; rr++) {
      for (int cc = 0; cc < 3; cc++) {
        D_lambda(rr, cc) = 1;
      }
    }
    D_mu.resize(6, 6);
    D_mu.clear();
    for (int rr = 0; rr < 6; rr++) {
      D_mu(rr, rr) = rr < 3 ? 2 : 1;
    }
    D = lambda * D_lambda + mu * D_mu;
 
    int tag_length = 9;
    double def_VAL[tag_length];
    bzero(def_VAL, tag_length * sizeof(double));
    Tag th_stress;
    CHKERR postProcMesh.tag_get_handle("STRESS", 9, MB_TYPE_DOUBLE, th_stress,
                                       MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
 
    Tag th_id;
    int def_block_id = -1;
    CHKERR postProcMesh.tag_get_handle("BLOCK_ID", 1, MB_TYPE_INTEGER, th_id,
                                       MB_TAG_CREAT | MB_TAG_SPARSE,
                                       &def_block_id);
    CHKERR postProcMesh.tag_clear_data(th_id, postProcElements, &id);
 
    VectorDouble strain;
    VectorDouble stress;
    MatrixDouble Stress;
 
    int nb_gauss_pts = data.getN().size1();
    if (mapGaussPts.size() != (unsigned int)nb_gauss_pts) 
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
 
    if (!fieldGradientPtr ||
        fieldGradientPtr->size1() != (unsigned int)nb_gauss_pts ||
        fieldGradientPtr->size2() != 9) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Gradient of field %s is not available", rowFieldName.c_str());
    }
    const auto &gradient = *fieldGradientPtr;
 
    for (int gg = 0; gg != nb_gauss_pts; ++gg) {
 
      strain.resize(6);
      strain[0] = gradient(gg, 0);
      strain[1] = gradient(gg, 4);
      strain[2] = gradient(gg, 8);
      strain[3] = gradient(gg, 1) + gradient(gg, 3);
      strain[4] = gradient(gg, 5) + gradient(gg, 7);
      strain[5] = gradient(gg, 2) + gradient(gg, 6);
 
      if (!isFieldDisp) {
        strain[0] -= 1.0;
        strain[1] -= 1.0;
        strain[2] -= 1.0;
      }
 
      stress.resize(6);
      noalias(stress) = prod(D, strain);
 
      Stress.resize(3, 3);
      Stress(0, 0) = stress[0];
      Stress(1, 1) = stress[1];
      Stress(2, 2) = stress[2];
      Stress(0, 1) = Stress(1, 0) = stress[3];
      Stress(1, 2) = Stress(2, 1) = stress[4];
      Stress(2, 0) = Stress(0, 2) = stress[5];
 
      CHKERR postProcMesh.tag_set_data(th_stress, &mapGaussPts[gg], 1,
                                       &Stress(0, 0));
    }
 
    MoFEMFunctionReturn(0);
  }
};
 
/// \deprecated Class name with spelling mistake
DEPRECATED typedef PostProcHookStress PostPorcHookStress;