TestingContact.hpp

Go to the documentation of this file.

/** Final-state pressure checks for the Hertz-contact tests. */
 
namespace ContactOps {
 
/**
 * Compare the magnitude of the contact traction against an analytical
 * pressure distribution, integrated as an L2 error over the contact
 * boundary, using the final accepted field state.
 */
struct ContactPressureNormCheck {
 
  ContactPressureNormCheck(MoFEM::Interface &m_field, double tol_norm,
                           std::string check_name)
      : mField(m_field), tolNorm(tol_norm), checkName(std::move(check_name)) {
  }
 
  virtual ~ContactPressureNormCheck() = default;
 
  MoFEMErrorCode operator()(DM dm) {
    MoFEMFunctionBegin;
 
    Range contact_range;
    for (auto m : mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
             std::regex((boost::format("%s(.*)") % "CONTACT").str()))) {
      auto meshset = m->getMeshset();
      Range contact_meshset_range;
      CHKERR mField.get_moab().get_entities_by_dimension(
          meshset, SPACE_DIM - 1, contact_meshset_range, true);
      CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
          contact_meshset_range);
      contact_range.merge(contact_meshset_range);
    }
    auto contact_range_ptr = boost::make_shared<Range>(contact_range);
 
    struct OpCalcTractionMagnitude : public BoundaryEleOp {
      OpCalcTractionMagnitude(boost::shared_ptr<MatrixDouble> traction_ptr,
                              boost::shared_ptr<VectorDouble> mag_ptr)
          : BoundaryEleOp(NOSPACE, OPSPACE), tractionPtr(traction_ptr),
            magPtr(mag_ptr) {}
      MoFEMErrorCode doWork(int, EntityType, EntitiesFieldData::EntData &) {
        MoFEMFunctionBegin;
        const auto nb_gauss_pts = getGaussPts().size2();
        magPtr->resize(nb_gauss_pts);
        using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
        auto t_traction = MatrixSizeHelper<
            GetFTensor1FromMatType<SPACE_DIM, -1, DL>, DL>::get(
                *tractionPtr, nb_gauss_pts)();
        auto t_mag = getFTensor0FromVec(*magPtr);
        for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
          t_mag = t_traction.l2();
          ++t_traction;
          ++t_mag;
        }
        MoFEMFunctionReturn(0);
      }
 
    private:
      boost::shared_ptr<MatrixDouble> tractionPtr;
      boost::shared_ptr<VectorDouble> magPtr;
    };
 
    auto fe = boost::make_shared<BoundaryEle>(mField);
    CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
        fe->getOpPtrVector(), {HDIV}, "GEOMETRY");
    // We have to integrate on curved face geometry, thus integration weight
    // have to adjusted.
    fe->getOpPtrVector().push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
    fe->getRuleHook = [](int, int, int approx_order) {
      return 2 * approx_order + geom_order - 1;
    };
 
    struct OpResetTraction : public BoundaryEleOp {
      boost::shared_ptr<MatrixDouble> tractionPtr;
      OpResetTraction(boost::shared_ptr<MatrixDouble> ptr)
          : BoundaryEleOp(NOSPACE, OPSPACE), tractionPtr(ptr) {}
      MoFEMErrorCode doWork(int, EntityType, EntitiesFieldData::EntData &) override {
        MoFEMFunctionBegin;
        using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
        MatrixSizeHelper<GetFTensor1FromMatType<SPACE_DIM, -1, DL>, DL>::size(
            *tractionPtr, getGaussPts().size2());
        tractionPtr->clear();
        MoFEMFunctionReturn(0);
      }
    };
    auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
    fe->getOpPtrVector().push_back(
        new OpResetTraction(common_data_ptr->contactTractionPtr()));
    fe->getOpPtrVector().push_back(
        new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
            "SIGMA", common_data_ptr->contactTractionPtr()));
 
    auto mag_traction_ptr = boost::make_shared<VectorDouble>();
    fe->getOpPtrVector().push_back(new OpCalcTractionMagnitude(
        common_data_ptr->contactTractionPtr(), mag_traction_ptr));
 
    auto analytical_pressure_ptr = boost::make_shared<VectorDouble>();
    fe->getOpPtrVector().push_back(new OpGetTensor0fromFunc(
        analytical_pressure_ptr,
        [this](double x, double y, double z) {
          return analyticalPressure(x, y, z);
        }));
 
    enum { NORM, LAST_NORM };
    auto norm_vec = createVectorMPI(
        mField.get_comm(), (mField.get_comm_rank() == 0) ? LAST_NORM : 0,
        LAST_NORM);
    fe->getOpPtrVector().push_back(new OpCalcNormL2Tensor0(
        mag_traction_ptr, norm_vec, NORM, analytical_pressure_ptr,
        contact_range_ptr));
 
    CHKERR VecZeroEntries(norm_vec);
    CHKERR DMoFEMLoopFiniteElements(dm, "bFE", fe);
    CHKERR VecAssemblyBegin(norm_vec);
    CHKERR VecAssemblyEnd(norm_vec);
 
    double norm = 0.;
    if (!mField.get_comm_rank()) {
      const double *norm_ptr;
      CHKERR VecGetArrayRead(norm_vec, &norm_ptr);
      norm = std::sqrt(norm_ptr[NORM]);
      CHKERR VecRestoreArrayRead(norm_vec, &norm_ptr);
    }
    CHKERR MPI_Bcast(&norm, 1, MPI_DOUBLE, 0, mField.get_comm());
 
    CHKERR mField.getInterface<TestingInterface>()
        ->checkError<TestType::SCALAR>(checkName, norm, 0., {tolNorm, 0.});
 
    MoFEMFunctionReturn(0);
  }
 
protected:
  /** Analytical pressure at a boundary gauss point's physical coordinates. */
  virtual double analyticalPressure(double x, double y, double z) const = 0;
 
  MoFEM::Interface &mField;
  double tolNorm;
  std::string checkName;
};
 
/** Test 2: Hertz cylinder-on-plane contact, plane strain. */
struct HertzPlaneStrainPressureNormCheck : public ContactPressureNormCheck {
  using ContactPressureNormCheck::ContactPressureNormCheck;
 
protected:
  double analyticalPressure(double x, double y, double) const override {
    double E_star = young_modulus / (1 - poisson_ratio * poisson_ratio);
    double R = 100.;
    double a = 1;
    double r = std::sqrt((x * x) + (y * y));
    if (r > a)
      return 0.;
    return E_star / (2. * R) * std::sqrt(a * a - r * r);
  }
};
 
/** Test 5: Hertz sphere-on-plane contact, axisymmetric. */
struct HertzAxisymmetricPressureNormCheck : public ContactPressureNormCheck {
  using ContactPressureNormCheck::ContactPressureNormCheck;
 
protected:
  double analyticalPressure(double x, double y, double) const override {
    double E_star = young_modulus / (1 - poisson_ratio * poisson_ratio);
    double R = 100.;
    double d = 0.01;
    double F = (4. / 3.) * E_star * std::sqrt(R) * std::pow(d, 1.5);
    double a = std::pow((3. * F * R) / (4. * E_star), 1. / 3.);
    double p_max = (3. * F) / (2. * M_PI * a * a);
    double r = std::sqrt((x * x) + (y * y));
    if (r > a)
      return 0.;
    return p_max * std::sqrt(1 - ((r * r) / (a * a)));
  }
};
 
} // namespace ContactOps