EshelbianTestingMonitor.cpp

Go to the documentation of this file.

/** @file
    @brief TestingInterface checks for the Eshelbian solvers.
    @ingroup EshelbianPlasticty
 
    Two kinds of checks are provided. Monitor checks are registered with
    TestingInterface::addTSMonitor and fire after every accepted TS step, once
    EshelbianMonitor has refreshed reactions and energy; they suit quantities
    that must hold along the whole load path. Final checks are registered with
    TestingInterface::addTestCommandLine and run by a single
    TestingInterface::runTests() call issued by the solver after it finishes.
    Both kinds are selected on the command line by "-test_name" or "-test_nb".
*/
 
struct EshelbianTesting {
 
  using MonitorPtr = boost::shared_ptr<EshelbianMonitor>;
 
  /// Register per-step checks. Call after TSMonitorSet(ts, TsMonitorSet, ...)
  /// so PETSc invokes them after EshelbianMonitor::postProcess.
  static MoFEMErrorCode addMonitorTests(EshelbianCore &ep, TS ts,
                                        MonitorPtr monitor);
 
  /// Register final-state checks; the solver runs them with runTests().
  static MoFEMErrorCode addFinalTests(EshelbianCore &ep, MonitorPtr monitor);
 
private:
  static MoFEMErrorCode getReactionBlockName(std::string &name);
  static MoFEMErrorCode getReaction(MonitorPtr monitor,
                                    const std::string &block, int component,
                                    double &value);
};
 
MoFEMErrorCode EshelbianTesting::getReactionBlockName(std::string &name) {
  MoFEMFunctionBegin;
  char block_name[255] = "FIX_ALL";
  CHKERR PetscOptionsGetString(PETSC_NULLPTR, "", "-test_reaction_block",
                               block_name, sizeof(block_name), PETSC_NULLPTR);
  name = block_name;
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianTesting::getReaction(MonitorPtr monitor,
                                             const std::string &block,
                                             int component, double &value) {
  MoFEMFunctionBegin;
  auto &reactions = monitor->getReactionMap();
  auto it = reactions.find(block);
  if (it == reactions.end())
    SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
            "No reaction computed for block %s", block.c_str());
  value = it->second[component];
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianTesting::addMonitorTests(EshelbianCore &ep, TS ts,
                                                 MonitorPtr monitor) {
  MoFEMFunctionBegin;
  auto testing = ep.mField.getInterface<TestingInterface>();
 
  // Bar with imposed external strain: the L2 displacement at x = 2.5 is 0.25.
  CHKERR testing->addTSMonitor(
      "", 14, "external_strain", ts,
      [&ep, testing](TS, PetscInt, PetscReal, Vec) {
        MoFEMFunctionBegin;
        const std::array<double, 3> point{2.5, 0., 0.};
        auto w_ptr = boost::make_shared<MatrixDouble>();
        auto evaluator = ep.mField.getInterface<FieldEvaluatorInterface>();
        auto data = evaluator->getData<VolumeElementForcesAndSourcesCore>(
            point.data(), 1);
        auto fe = data->feMethodPtr;
        fe->getRuleHook = [](int, int, int) { return -1; };
        // EP prepares its USER_BASE bubble field even when only the L2
        // displacement is evaluated; no bubble values are consumed here.
        fe->getUserPolynomialBase() =
            boost::make_shared<CGGUserPolynomialBase>(nullptr, true);
        CHKERR EshelbianPlasticity::AddHOOps<
            SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
            fe->getOpPtrVector(), {H1, L2}, ep.materialH1Positions,
            ep.frontAdjEdges);
        fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
            ep.spatialL2Disp, w_ptr, MBTET));
        CHKERR evaluator->buildTree<SPACE_DIM>(data, "EP");
        CHKERR evaluator->evalFEAtThePoint<SPACE_DIM>(
            point.data(), 1e-12, getProblemPtr(ep.dmElastic)->getName(), "EP",
            data, ep.mField.get_comm_rank(), ep.mField.get_comm_rank(), nullptr,
            MF_EXIST, QUIET);
        const bool found = fe->gaussPts.size2() == 1;
        // The point is owned by one rank; share its value so the check is
        // collective with identical data everywhere.
        double owner_data[2] = {0., 0.};
        if (found && w_ptr->size1() && w_ptr->size2()) {
          owner_data[0] = (*w_ptr)(0, 0);
          owner_data[1] = 1.;
        }
        CHKERR MPI_Allreduce(MPI_IN_PLACE, owner_data, 2, MPI_DOUBLE, MPI_SUM,
                             ep.mField.get_comm());
        if (owner_data[1] == 0.)
          SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
                  "Evaluation point (2.5, 0, 0) not found on any rank");
        CHKERR testing->checkValue<TestType::POINTWISE>(
            "external strain displacement", owner_data[0] / owner_data[1], 0.25,
            {1e-5, 0.});
        MoFEMFunctionReturn(0);
      });
 
  // Cook cantilever (cook.jou): end face height 16, thickness 5, traction
  // profile 4*eta*(1-eta) with eta = (y-44)/16 and load scale t. Its integral
  // 16*5*(2/3)*t is balanced by the clamp y reaction at every step.
  CHKERR testing->addTSMonitor(
      "", 20, "cook_reaction", ts,
      [monitor, testing](TS, PetscInt, PetscReal t, Vec) {
        MoFEMFunctionBegin;
        std::string block;
        CHKERR getReactionBlockName(block);
        double reaction;
        CHKERR getReaction(monitor, block, 1, reaction);
        CHKERR testing->checkValue<TestType::SCALAR>(
            "cook clamp reaction", reaction, -(160.0 / 3.0) * t, {1e-7, 0.});
        MoFEMFunctionReturn(0);
      });
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianTesting::addFinalTests(EshelbianCore &ep,
                                               MonitorPtr monitor) {
  MoFEMFunctionBegin;
  auto testing = ep.mField.getInterface<TestingInterface>();
 
  auto add_reaction_test = [&](int number, const std::string &name,
                               const char *block_or_null, int component,
                               double expected, double tolerance) {
    return testing->addTestCommandLine(
        "", number, name,
        [monitor, testing, name, block_or_null, component, expected,
         tolerance]() {
          MoFEMFunctionBegin;
          std::string block;
          if (block_or_null)
            block = block_or_null;
          else
            CHKERR getReactionBlockName(block);
          double reaction;
          CHKERR getReaction(monitor, block, component, reaction);
          CHKERR testing->checkValue<TestType::SCALAR>(name, reaction, expected,
                                                       {tolerance, 0.});
          MoFEMFunctionReturn(0);
        });
  };
 
  // L = 5, ty = tz = 1, A = 1, UDL 0.1 on the top face: moment at x = 0.
  CHKERR add_reaction_test(15, "reaction_moment", nullptr, 5, 0.25, 1e-3);
  // L = 5, A = 1, E = 1000 exp(0.1 x), u = 0.1: F = 1000*0.1/(1-exp(-0.5)).
  CHKERR add_reaction_test(16, "reaction_force", nullptr, 0, -25.4149, 1e-4);
  CHKERR add_reaction_test(19, "spring_reaction", "SPRING_BC", 2, 2.016e1,
                           1e-3);
 
  CHKERR testing->addTestCommandLine(
      "", 17, "strain_energy", [monitor, testing]() {
        MoFEMFunctionBegin;
        CHKERR testing->checkValue<TestType::SCALAR>(
            "strain energy", *monitor->getEnergy(), 1.27096, {1e-5, 0.});
        MoFEMFunctionReturn(0);
      });
 
  // Load-factor solver: the factor reached at physical step 7 is fixed by
  // the Griffith energy; other steps are not checked.
  CHKERR testing->addTestCommandLine("", 18, "load_factor", [&ep, testing]() {
    MoFEMFunctionBegin;
    if (EshelbianCore::physicalStepNumber == 7)
      CHKERR testing->checkValue<TestType::SCALAR>(
          "load factor at step 7", ep.loadFactor, 16.0393, {1e-5, 0.});
    MoFEMFunctionReturn(0);
  });
 
  MoFEMFunctionReturn(0);
}