EshelbianTestingMonitor.cpp

Go to the documentation of this file.

/** @file
    @brief Contains definition of EshelbianTestingMonitor class.
    @ingroup EshelbianPlasticty
*/
 
struct EshelbianTestingMonitor : public FEMethod {
  using SetPtsData = FieldEvaluatorInterface::SetPtsData;
 
  EshelbianTestingMonitor(EshelbianCore &ep,
                          boost::shared_ptr<EshelbianMonitor> base_monitor_ptr)
      : baseMonitorPtr(base_monitor_ptr), eP(baseMonitorPtr->getEpCore()),
        ptsHashMap(baseMonitorPtr->getHashMap()),
        reactionForcesMap(baseMonitorPtr->getReactionMap()),
        gEnergy(baseMonitorPtr->getEnergy()),
        dataFieldEval(baseMonitorPtr->getDataField()) {
 
    PetscBool test_cook_flg = PETSC_FALSE;
    PetscInt atom_test = 0;
    CHK_THROW_MESSAGE(PetscOptionsGetBool(PETSC_NULLPTR, "", "-test_cook_pts",
                                          &test_cook_flg, PETSC_NULLPTR),
                      "get post proc points");
    CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-atom_test", &atom_test,
                              PETSC_NULLPTR);
 
    if (test_cook_flg) {
      ptsHashMap["Point A"] = {48., 60., 4.999};
      ptsHashMap["Point B"] = {48. / 2., 44. + (60. - 44.) / 2., 0.};
      ptsHashMap["Point C"] = {48. / 2., (44. - 0.) / 2., 0.};
    }
    if (atom_test == 14) {
      // Points for atom test 14: check external strain
      ptsHashMap["Point (2.5, 0., 0.)"] = {2.5, 0, 0.};
    }
  }
 
  MoFEMErrorCode preProcess() { return 0; }
 
  MoFEMErrorCode operator()() { return 0; }
 
  MoFEMErrorCode postProcess() {
    MoFEMFunctionBegin;
 
    MOFEM_LOG("EP", Sev::inform) << "Testing Monitor postProcess";
 
    PetscInt atom_test = 0;
    char reaction_block_name[255] = "FIX_ALL";
    CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-atom_test", &atom_test,
                              PETSC_NULLPTR);
    CHKERR PetscOptionsGetString(PETSC_NULLPTR, "", "-atom_test_reaction_block",
                                 reaction_block_name, 255, PETSC_NULLPTR);
 
    switch (atom_test) {
    case 14:
      // Points for external strain
      for (auto &pts : ptsHashMap) {
        CHKERR checkExternalStrain(pts.second, pts.first, atom_test);
      }
      break;
    case 15:
      if (ts_t > 0.0) {
        // L = 5, ty = tz = 1, A = 1, UDL Top face 0.1,
        // Moment at x = 0 should be 0.25.
        if (std::abs(reactionForcesMap[reaction_block_name][5] - 0.25) > 1e-3) {
          MOFEM_LOG_SYNCHRONISE(eP.mField.get_comm());
          SETERRQ(
              PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
              "Atom test 15 failed: reaction moment does not match expected "
              "value. Got [%3.6e], expected [0.25].",
              reactionForcesMap[reaction_block_name][5]);
        }
      }
      break;
    case 16:
      if (ts_t > 0.0) {
        // L = 5, ty = tz = 1, A = 1, E = 1000 * exp(0.1*x), u = 0.1, F = 1000 *
        // 0.1/(1-exp(-0.5)) = 25.4149, so reaction force should
        // be -25.4149
        if (std::abs(reactionForcesMap[reaction_block_name][0] - (-25.4149)) >
            1e-4) {
          MOFEM_LOG_SYNCHRONISE(eP.mField.get_comm());
          SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                  "Atom test 16 failed: reaction force does not match expected "
                  "value. Got [%3.6e, %3.6e, %3.6e], expected [-25.4149, 0.0, "
                  "0.0].",
                  reactionForcesMap[reaction_block_name][0],
                  reactionForcesMap[reaction_block_name][1],
                  reactionForcesMap[reaction_block_name][2]);
        }
      }
      break;
    case 17:
      if (ts_t > 0.0) {
        if (std::abs(*gEnergy - 1.27096) > 1e-5) {
          MOFEM_LOG_SYNCHRONISE(eP.mField.get_comm());
          SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                  "Atom test 17 failed: strain energy does not match expected "
                  "value. Got %3.6e, expected 1.27096.",
                  *gEnergy);
        }
      }
      break;
    case 18:
      if (ts_step == 7) {
        if (std::abs(eP.loadFactor - 16.0393) > 1e-5) {
          MOFEM_LOG_SYNCHRONISE(eP.mField.get_comm());
          SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                  "Atom test 18 failed: load factor does not match expected "
                  "value. Got %3.6e, expected 16.093.",
                  eP.loadFactor);
        }
      }
      break;
    case 19:
      if (ts_t > 0.0) {
        if (std::abs(reactionForcesMap["SPRING_BC"][2] - 2.016e1) > 1e-3) {
          MOFEM_LOG_SYNCHRONISE(eP.mField.get_comm());
          SETERRQ(
              PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
              "Atom test 19 failed: reaction force does not match expected "
              "value. Got [%3.6e], expected [2.016e1].",
              reactionForcesMap["SPRING_BC"][2]);
        }
      }
      break;
    case 20:
      if (ts_t > 0.0) {
        // cook.jou: end face height 16, thickness 5. The COOK traction
        // profile is 4*eta*(1-eta), with eta=(y-44)/16 and load scale t.
        // Its integral is 16*5*(2/3)*t, balanced by the clamp y reaction.
        const double expected_reaction = -(160.0 / 3.0) * ts_t;
        const double reaction = reactionForcesMap[reaction_block_name][1];
        if (!(std::abs(reaction - expected_reaction) <= 1e-7)) {
          MOFEM_LOG_SYNCHRONISE(eP.mField.get_comm());
          SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                  "Atom test 20 failed: Cook clamp reaction does not balance "
                  "the applied traction. Got %.12g, expected %.12g.",
                  reaction, expected_reaction);
        }
        MOFEM_LOG_C("EP", Sev::inform,
                    "Cook cantilever reaction balance passed at time %g", ts_t);
      }
      break;
    default:
      break;
    }
 
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode checkExternalStrain(std::array<double, 3> point,
                                     std::string str, PetscInt atom_test);
 
  protected:
  boost::shared_ptr<EshelbianMonitor> baseMonitorPtr;
  EshelbianCore &eP;
  std::map<std::string, std::array<double, 6>> &reactionForcesMap;
  boost::shared_ptr<SetPtsData> dataFieldEval;
  std::map<std::string, std::array<double, 3>> &ptsHashMap;
  boost::shared_ptr<double> gEnergy;
};
 
MoFEMErrorCode EshelbianTestingMonitor::checkExternalStrain(
    std::array<double, 3> point, std::string str, PetscInt atom_test) {
  MoFEMFunctionBegin;
 
  dataFieldEval->setEvalPoints(point.data(), point.size() / 3);
  auto vec = createVectorMPI(eP.mField.get_comm(), 1, 1);
  auto w_l2_at_pts = eP.dataAtPts->getSmallWL2AtPts();
  w_l2_at_pts->resize(0, 0, false);
 
  if (auto fe_ptr = dataFieldEval->feMethodPtr) {
    CHKERR eP.mField.getInterface<FieldEvaluatorInterface>()
        ->evalFEAtThePoint<SPACE_DIM>(
            point.data(), 1e-12, problemPtr->getName(), "EP", dataFieldEval,
            eP.mField.get_comm_rank(), eP.mField.get_comm_rank(), nullptr,
            MF_EXIST, QUIET);
  }
  if (w_l2_at_pts->size1() == 0 || w_l2_at_pts->size2() == 0) {
    CHKERR VecSetValue(vec, 0, 0.0, ADD_VALUES);
  } else if (w_l2_at_pts->size1() == 1 || w_l2_at_pts->size2() == 1) {
    auto add = [&]() {
      std::ostringstream s;
      s << str << " elem " << getFEEntityHandle() << " ";
      return s.str();
    };
    MOFEM_LOG("EPSYNC", Sev::inform)
        << add() << "comm rank " << eP.mField.get_comm_rank();
    MOFEM_LOG("EPSYNC", Sev::inform)
        << add() << "point " << getVectorAdaptor(point.data(), 3);
    MOFEM_LOG("EPSYNC", Sev::inform)
        << add() << "w " << *w_l2_at_pts;
    double disp_at_point = (*w_l2_at_pts)(0, 0);
    CHKERR VecSetValue(vec, 0, disp_at_point, ADD_VALUES);
  } else {
    SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
            "Unexpected displacement data shape (%ld, %ld).",
            static_cast<long>(w_l2_at_pts->size1()),
            static_cast<long>(w_l2_at_pts->size2()));
  }
 
  MOFEM_LOG_SEVERITY_SYNC(eP.mField.get_comm(), Sev::inform);
  if (ts_t > 0.0) {
    CHKERR VecAssemblyBegin(vec);
    CHKERR VecAssemblyEnd(vec);
    double error;
    PetscInt idx = 0;
    if (eP.mField.get_comm_rank() == 0) {
      CHKERR VecGetValues(vec, 1, &idx, &error);
    }
    MPI_Bcast(&error, 1, MPI_DOUBLE, 0, PETSC_COMM_WORLD);
    // Check if the displacement is correct for the applied external
    // strain For the bar problem, we expect a displacement of 0.25 at
    // point A
    if (std::abs(error - 0.25) > 1e-5) {
      SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
              "Atom test %d failed: wrong displacement %.12g.", atom_test,
              error);
    }
  }
  MoFEMFunctionReturn(0);
}