PostProcStresses.hpp

Go to the documentation of this file.

/** \file PostProcStresses.hpp
 * \brief Post-processing stresses for non-linear analysis
 * \ingroup nonlinear_elastic_elem
 *
 * Implementation of method for post-processing stresses.
 */
 
 
 
#ifndef __POSTPROCSTRESSES_HPP__
#define __POSTPROCSTRESSES_HPP__
 
#ifndef WITH_ADOL_C
#error "MoFEM need to be compiled with ADOL-C"
#endif
 
struct PostProcStress
    : public MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator {
 
  moab::Interface &postProcMesh;
  std::vector<EntityHandle> &mapGaussPts;
  Range &postProcElements;
 
  NonlinearElasticElement::BlockData &dAta;
  boost::shared_ptr<MatrixDouble> fieldValuesPtr;
  boost::shared_ptr<MatrixDouble> meshPositionsPtr;
  boost::shared_ptr<MatrixDouble> fieldGradientPtr;
  const std::string meshPositionsFieldName;
  const bool fieldDisp;
  const bool replaceNonANumberByMaxValue;
  const double maxVal;
  const bool printCauchy;
 
  PostProcStress(moab::Interface &post_proc_mesh,
                 std::vector<EntityHandle> &map_gauss_pts,
                 Range &post_proc_elements,
                 const std::string field_name,
                 NonlinearElasticElement::BlockData &data,
                 boost::shared_ptr<MatrixDouble> field_values_ptr,
                 boost::shared_ptr<MatrixDouble> mesh_positions_ptr,
                 boost::shared_ptr<MatrixDouble> field_gradient_ptr,
                 const std::string mesh_positions_field_name =
                     "MESH_NODE_POSITIONS",
                 const bool field_disp = false,
                 const bool replace_nonanumber_by_max_value = false,
                 const double max_val = 1e16,
                 const bool print_cauchy_stress = false)
      : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
            field_name, ForcesAndSourcesCore::UserDataOperator::OPROW),
        postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
        postProcElements(post_proc_elements), dAta(data),
        fieldValuesPtr(field_values_ptr), meshPositionsPtr(mesh_positions_ptr),
        fieldGradientPtr(field_gradient_ptr),
        meshPositionsFieldName(mesh_positions_field_name),
        fieldDisp(field_disp),
        replaceNonANumberByMaxValue(replace_nonanumber_by_max_value),
        maxVal(max_val), printCauchy(print_cauchy_stress) {}
 
  NonlinearElasticElement::CommonData nonLinearElementCommonData;
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data) {
    MoFEMFunctionBegin;
 
    if (type != MBVERTEX)
      MoFEMFunctionReturnHot(0);
    if (data.getIndices().size() == 0)
      MoFEMFunctionReturnHot(0);
    if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
        dAta.tEts.end()) {
      MoFEMFunctionReturnHot(0);
    }
 
    const auto &dof_ptr = data.getFieldDofs()[0];
 
    int id = dAta.iD;
 
    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);
 
    string tag_name_piola1 = dof_ptr->getName() + "_PIOLA1_STRESS";
    string tag_name_energy = dof_ptr->getName() + "_ENERGY_DENSITY";
 
    int tag_length = 9;
    double def_VAL[tag_length];
    bzero(def_VAL, tag_length * sizeof(double));
    Tag th_piola1, th_energy, th_cauchy;
    CHKERR postProcMesh.tag_get_handle(tag_name_piola1.c_str(), tag_length,
                                       MB_TYPE_DOUBLE, th_piola1,
                                       MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
    CHKERR postProcMesh.tag_get_handle(tag_name_energy.c_str(), 1,
                                       MB_TYPE_DOUBLE, th_energy,
                                       MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
 
    if (printCauchy) {
      string tag_name_cauchy = "MED_" + dof_ptr->getName() + "_CAUCHY_STRESS";
      CHKERR postProcMesh.tag_get_handle(tag_name_cauchy.c_str(), tag_length,
                                         MB_TYPE_DOUBLE, th_cauchy,
                                         MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
    }
 
    int nb_gauss_pts = data.getN().size1();
    if (mapGaussPts.size() != (unsigned int)nb_gauss_pts) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Nb. of integration points is not equal to number points on "
              "post-processing mesh");
    }
    if (!fieldGradientPtr ||
        fieldGradientPtr->size1() != (unsigned int)nb_gauss_pts ||
        fieldGradientPtr->size2() != 9) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Gradient of field not found, filed <%s> not found",
              rowFieldName.c_str());
    }
 
    auto make_field_values = [](auto values_ptr) {
      std::vector<VectorDouble> values;
      if (!values_ptr)
        return values;
      values.resize(values_ptr->size1());
      for (size_t gg = 0; gg != values_ptr->size1(); ++gg) {
        values[gg].resize(values_ptr->size2(), false);
        for (size_t rr = 0; rr != values_ptr->size2(); ++rr) {
          values[gg][rr] = (*values_ptr)(gg, rr);
        }
      }
      return values;
    };
 
    auto make_gradients = [](auto gradient_ptr) {
      std::vector<MatrixDouble> gradients;
      if (!gradient_ptr)
        return gradients;
      const size_t field_rank = gradient_ptr->size2() / 3;
      gradients.resize(gradient_ptr->size1());
      for (size_t gg = 0; gg != gradient_ptr->size1(); ++gg) {
        gradients[gg].resize(field_rank, 3, false);
        for (size_t rr = 0; rr != field_rank; ++rr) {
          for (size_t cc = 0; cc != 3; ++cc) {
            gradients[gg](rr, cc) = (*gradient_ptr)(gg, 3 * rr + cc);
          }
        }
      }
      return gradients;
    };
 
    std::map<std::string, std::vector<VectorDouble>> field_map{
        {rowFieldName, make_field_values(fieldValuesPtr)},
        {meshPositionsFieldName, make_field_values(meshPositionsPtr)}};
    std::map<std::string, std::vector<MatrixDouble>> gradient_map{
        {rowFieldName, make_gradients(fieldGradientPtr)}};
 
    MatrixDouble3by3 H, invH;
    double detH;
 
    dAta.materialDoublePtr->commonDataPtr = &nonLinearElementCommonData;
    dAta.materialDoublePtr->opPtr = this;
    CHKERR dAta.materialDoublePtr->getDataOnPostProcessor(field_map,
                                                          gradient_map);
 
    nonLinearElementCommonData.dataAtGaussPts = field_map;
    nonLinearElementCommonData.gradAtGaussPts = gradient_map;
 
    MatrixDouble3by3 maxP(3, 3);
    maxP.clear();
 
    for (int gg = 0; gg != nb_gauss_pts; ++gg) {
 
      dAta.materialDoublePtr->gG = gg;
      dAta.materialDoublePtr->F.resize(3, 3);
      noalias(dAta.materialDoublePtr->F) =
          gradient_map[rowFieldName][gg];
      if (fieldDisp) {
        for (int dd = 0; dd != 3; dd++) {
          dAta.materialDoublePtr->F(dd, dd) += 1;
        }
      }
      if (gradient_map[meshPositionsFieldName].size() ==
          (unsigned int)nb_gauss_pts) {
        H.resize(3, 3);
        invH.resize(3, 3);
        noalias(H) = gradient_map[meshPositionsFieldName][gg];
        detH = determinantTensor3by3(H);
        CHKERR invertTensor3by3(H, detH, invH);
        noalias(dAta.materialDoublePtr->F) =
            prod(dAta.materialDoublePtr->F, invH);
      }
 
      int nb_active_variables = 9;
      CHKERR dAta.materialDoublePtr->setUserActiveVariables(
          nb_active_variables);
      CHKERR dAta.materialDoublePtr->calculateP_PiolaKirchhoffI(
          dAta, getNumeredEntFiniteElementPtr());
      CHKERR dAta.materialDoublePtr->calculateElasticEnergy(
          dAta, getNumeredEntFiniteElementPtr());
      CHKERR postProcMesh.tag_set_data(th_piola1, &mapGaussPts[gg], 1,
                                       &dAta.materialDoublePtr->P(0, 0));
      CHKERR postProcMesh.tag_set_data(th_energy, &mapGaussPts[gg], 1,
                                       &dAta.materialDoublePtr->eNergy);
      if (printCauchy) {
        dAta.materialDoublePtr->sigmaCauchy.resize(3, 3);
        CHKERR dAta.materialDoublePtr->calculateCauchyStress(
            dAta, getNumeredEntFiniteElementPtr());
        CHKERR postProcMesh.tag_set_data(
            th_cauchy, &mapGaussPts[gg], 1,
            &dAta.materialDoublePtr->sigmaCauchy(0, 0));
      }
    }
 
    if (replaceNonANumberByMaxValue) {
      MatrixDouble3by3 P(3, 3);
      for (int gg = 0; gg != nb_gauss_pts; ++gg) {
        double val_energy;
        CHKERR postProcMesh.tag_get_data(th_energy, &mapGaussPts[gg], 1,
                                         &val_energy);
        if (!std::isnormal(val_energy)) {
          CHKERR postProcMesh.tag_set_data(th_energy, &mapGaussPts[gg], 1,
                                           &maxVal);
          CHKERR postProcMesh.tag_get_data(th_piola1, &mapGaussPts[gg], 1,
                                           &P(0, 0));
          for (unsigned int r = 0; r != P.size1(); ++r) {
            for (unsigned int c = 0; c != P.size2(); ++c) {
              if (!std::isnormal(P(r, c)))
                P(r, c) = copysign(maxVal, P(r, c));
            }
          }
          CHKERR postProcMesh.tag_set_data(th_piola1, &mapGaussPts[gg], 1,
                                           &P(0, 0));
        }
      }
    }
 
    MoFEMFunctionReturn(0);
  }
};
 
/// \deprecated Use PostProcStress
DEPRECATED typedef PostProcStress PostPorcStress;
 
#endif //__POSTPROCSTRESSES_HPP__