auxiliary_logarithmic_stress_patch_atom.cpp

Go to the documentation of this file.

/**
 * @file auxiliary_logarithmic_stress_patch_atom.cpp
 * @brief Homogeneous finite-deformation patch and boundary flux verification.
 */
 
#include <MoFEM.hpp>
using namespace MoFEM;
#include "AuxiliaryLogarithmicStressTestSetup.hpp"
#include <AuxiliaryLogarithmicStress.hpp>
#include <Lie.hpp>
using namespace EshelbianPlasticity;
 
namespace {
 
using VolumeElement = VolumeElementForcesAndSourcesCore;
using FaceElement = FaceElementForcesAndSourcesCore;
 
// Independent forward reference: diagonal logarithmic stretch followed by a
// finite rotation. The boundary callback constructs the same F with Rodrigues.
struct Reference {
  const double c10 = 1.7;
  const double bulk = 8.5;
  const double theta = .03;
  FTensor::Tensor1<double, 3> omega{.02, -.015, .01};
  FTensor::Tensor2<double, 3, 3> deformation, piola;
  FTensor::Tensor1<double, 5> deviator, stress;
 
  Reference() {
    FTENSOR_INDEXES(3, i, j, k);
    double total_angle = omega.l2();
    CHK_THROW_MESSAGE(PetscOptionsGetReal(nullptr, nullptr,
                                          "-test_patch_total_rotation",
                                          &total_angle, nullptr),
                      "get patch rotation");
    omega(i) *= total_angle / omega.l2();
    FTensor::Tensor2_symmetric<double, 3> t_h{.04, 0., 0., -.015, 0., .005};
    FTensor::Tensor2_symmetric<double, 3> t_u{
        std::exp(.04), 0., 0., std::exp(-.015), 0., std::exp(.005)};
    FTensor::Tensor2_symmetric<double, 3> t_inverse_u{
        std::exp(-.04), 0., 0., std::exp(.015), 0., std::exp(-.005)};
    const auto t_rotation = LieGroups::SO3::exp(omega, omega.l2());
    FTensor::Tensor2<double, 3, 3> t_inverse_transpose;
    deformation(i, j) = t_rotation(i, k) * t_u(k, j);
    t_inverse_transpose(i, j) = t_rotation(i, k) * t_inverse_u(k, j);
    const double jacobian = std::exp(theta);
    const double invariant = deformation(i, j) * deformation(i, j);
    piola(i, j) =
        (2. * c10 * std::exp(-2. * theta / 3.)) *
            (deformation(i, j) - (invariant / 3.) * t_inverse_transpose(i, j)) +
        (bulk * jacobian * (jacobian - 1.)) * t_inverse_transpose(i, j);
    deviator = Tensor2SymmetricDeviatorBasis::getCoordinates(t_h);
    FTensor::Tensor2_symmetric<double, 3> t_exp_2d{
        std::exp(.06), 0., 0., std::exp(-.05), 0., std::exp(-.01)};
    stress = Tensor2SymmetricDeviatorBasis::getCoordinates(t_exp_2d);
    FTensor::Index<'A', 5> A;
    stress(A) *= 2. * c10;
  }
};
 
struct Audit {
  enum Error { F, P, D, THETA, TD, OMEGA, W, FLUX, FACE_FLUX, COUNT };
  std::array<double, COUNT> error{};
  std::array<double, 3> reaction{};
  PetscInt volumePoints = 0;
  PetscInt faces = 0;
  double boundaryArea = 0.;
  MatrixDouble exactTraction;
};
 
struct OpCheckVolume : VolumeElement::UserDataOperator {
  OpCheckVolume(boost::shared_ptr<DataAtIntegrationPts> data,
                boost::shared_ptr<Audit> audit,
                boost::shared_ptr<Reference> reference)
      : VolumeElement::UserDataOperator(NOSPACE, OPSPACE), dataPtr(data),
        auditPtr(audit), referencePtr(reference) {}
 
  MoFEMErrorCode doWork(int, EntityType, EntData &) override {
    MoFEMFunctionBegin;
    FTENSOR_INDEXES(3, i, j);
    FTensor::Index<'A', 5> A;
    constexpr auto t_identity = FTensor::Kronecker_Delta<int>();
    const int points = getGaussPts().size2();
    if (points <= 1)
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "Patch requires multiple volume quadrature points");
    auto t_f = dataPtr->getFTensorSmallH(points);
    auto t_p = dataPtr->getFTensorApproxP(points);
    auto rotation = dataPtr->getFTensorRotMat(points);
    auto t_w = dataPtr->getFTensorSmallWL2(points);
    auto t_x = getFTensor1CoordsAtGaussPts();
    auto t_d = MatrixSizeHelper<GetFTensor1FromMatType<5, -1, DL>, DL>::get(
        *dataPtr->auxiliaryData->logDeviator, points)();
    auto t_theta = MatrixSizeHelper<GetFTensor1FromMatType<1, -1, DL>, DL>::get(
        *dataPtr->auxiliaryData->logJacobian, points)();
    auto t_stress =
        MatrixSizeHelper<GetFTensor1FromMatType<5, -1, DL>, DL>::get(
            *dataPtr->auxiliaryData->stress, points)();
    const auto &reference = *referencePtr;
    FTensor::Tensor2<double, 3, 3> t_matrix_error;
    FTensor::Tensor1<double, 3> t_vector_error;
    FTensor::Tensor1<double, 5> t_material_error;
    const auto update = [&](Audit::Error index, double error) {
      auditPtr->error[index] = std::max(
          auditPtr->error[index],
          std::isfinite(error) ? error
                               : std::numeric_limits<double>::infinity());
    };
    for (int gg = 0; gg != points; ++gg) {
      t_matrix_error(i, j) = t_f(i, j) - reference.deformation(i, j);
      update(Audit::F, std::sqrt(t_matrix_error(i, j) * t_matrix_error(i, j)));
      t_matrix_error(i, j) = t_p(i, j) - reference.piola(i, j);
      update(Audit::P, std::sqrt(t_matrix_error(i, j) * t_matrix_error(i, j)));
      t_material_error(A) = t_d(A) - reference.deviator(A);
      update(Audit::D, std::sqrt(t_material_error(A) * t_material_error(A)));
      update(Audit::THETA, std::abs(t_theta(0) - reference.theta));
      t_material_error(A) = t_stress(A) - reference.stress(A);
      update(Audit::TD, std::sqrt(t_material_error(A) * t_material_error(A)));
      // Compare physical rotations; branch-equivalent coefficients may differ.
      const auto expected_rotation =
          LieGroups::SO3::exp(reference.omega, reference.omega.l2());
      t_matrix_error(i, j) = rotation(i, j) - expected_rotation(i, j);
      update(Audit::OMEGA,
             std::sqrt(t_matrix_error(i, j) * t_matrix_error(i, j)));
      t_vector_error(i) =
          t_w(i) - (reference.deformation(i, j) - t_identity(i, j)) * t_x(j);
      update(Audit::W, t_vector_error.l2());
      ++t_f;
      ++t_p;
      ++rotation;
      ++t_w;
      ++t_x;
      ++t_d;
      ++t_theta;
      ++t_stress;
    }
    auditPtr->volumePoints += points;
    MoFEMFunctionReturn(0);
  }
 
  boost::shared_ptr<DataAtIntegrationPts> dataPtr;
  boost::shared_ptr<Audit> auditPtr;
  boost::shared_ptr<Reference> referencePtr;
};
 
struct OpPatchTraction : OpCalculateTractionFromSideEle {
  OpPatchTraction(boost::shared_ptr<DataAtIntegrationPts> data,
                  boost::shared_ptr<Audit> audit,
                  boost::shared_ptr<Reference> reference)
      : OpCalculateTractionFromSideEle(data), auditPtr(audit),
        referencePtr(reference) {}
 
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) override {
    MoFEMFunctionBegin;
    CHKERR OpCalculateTractionFromSideEle::doWork(side, type, data);
    if (getLoopSize() != 1)
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "Patch flux check requires an exterior face");
    FTENSOR_INDEXES(3, i, j);
    const int points = getGaussPts().size2();
    auto t_normal = getFTensor1NormalsAtGaussPts();
    auto t_exact =
        MatrixSizeHelper<GetFTensor1FromMatType<3, -1, DL>, DL>::size(
            auditPtr->exactTraction, points)();
    for (int gg = 0; gg != points; ++gg) {
      t_exact(i) = referencePtr->piola(i, j) *
                   (getSkeletonSense() * t_normal(j) / t_normal.l2());
      ++t_normal;
      ++t_exact;
    }
    MoFEMFunctionReturn(0);
  }
 
  boost::shared_ptr<Audit> auditPtr;
  boost::shared_ptr<Reference> referencePtr;
};
 
struct OpCheckFace : FaceElement::UserDataOperator {
  OpCheckFace(boost::shared_ptr<DataAtIntegrationPts> data,
              boost::shared_ptr<Audit> audit)
      : FaceElement::UserDataOperator(NOSPACE, OPSPACE), dataPtr(data),
        auditPtr(audit) {}
 
  MoFEMErrorCode doWork(int, EntityType, EntData &) override {
    MoFEMFunctionBegin;
    FTENSOR_INDEX(3, i);
    const int points = getGaussPts().size2();
    auto t_weight = getFTensor0IntegrationWeight();
    auto t_traction = dataPtr->getFTensorTraction(points);
    auto t_exact = MatrixSizeHelper<GetFTensor1FromMatType<3, -1, DL>, DL>::get(
        auditPtr->exactTraction, points)();
    FTensor::Tensor1<double, 3> t_error, t_face_error{0., 0., 0.};
    FTensor::Tensor1<double, 3> t_reaction{0., 0., 0.};
    double area = 0.;
    for (int gg = 0; gg != points; ++gg) {
      const double weight = t_weight * getMeasure();
      t_error(i) = t_traction(i) - t_exact(i);
      if (!std::isfinite(t_error.l2()))
        SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                "Nonfinite patch traction error");
      auditPtr->error[Audit::FLUX] =
          std::max(auditPtr->error[Audit::FLUX], t_error.l2());
      t_face_error(i) += weight * t_error(i);
      t_reaction(i) += weight * t_traction(i);
      area += weight;
      ++t_weight;
      ++t_traction;
      ++t_exact;
    }
    if (!(area > 0.))
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "Nonpositive patch face area");
    auditPtr->error[Audit::FACE_FLUX] =
        std::max(auditPtr->error[Audit::FACE_FLUX], t_face_error.l2() / area);
    auditPtr->reaction[0] += t_reaction(0);
    auditPtr->reaction[1] += t_reaction(1);
    auditPtr->reaction[2] += t_reaction(2);
    auditPtr->boundaryArea += area;
    ++auditPtr->faces;
    MoFEMFunctionReturn(0);
  }
 
  boost::shared_ptr<DataAtIntegrationPts> dataPtr;
  boost::shared_ptr<Audit> auditPtr;
};
 
MoFEMErrorCode prepareMesh() {
  MoFEMFunctionBegin;
  int ranks;
  CHKERR MPI_Comm_size(PETSC_COMM_WORLD, &ranks);
  if (ranks != 1)
    SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
            "Prepare the patch mesh on one rank before partitioning");
  char source[PETSC_MAX_PATH_LEN] = "";
  char output[PETSC_MAX_PATH_LEN] = "";
  CHKERR PetscOptionsGetString(nullptr, nullptr, "-file_name", source,
                               sizeof(source), nullptr);
  CHKERR PetscOptionsGetString(nullptr, nullptr, "-output_file", output,
                               sizeof(output), nullptr);
  if (!source[0] || !output[0] || std::string(source) == output)
    SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
            "Patch mesh preparation requires distinct source and output files");
  moab::Core moab;
  new ParallelComm(&moab, PETSC_COMM_WORLD);
  CHKERR moab.load_file(source);
  MoFEM::Core core(moab);
  MoFEM::Interface &m_field = core;
  auto *meshsets = m_field.getInterface<MeshsetsManager>();
  CHKERR meshsets->setMeshsetFromFile();
  std::vector<std::pair<CubitBCType, int>> old_sets;
  for (auto it = meshsets->getBegin(); it != meshsets->getEnd(); ++it)
    old_sets.emplace_back(CubitBCType(it->getMaskedBcTypeULong()),
                          it->getMeshsetId());
  for (const auto &entry : old_sets)
    CHKERR meshsets->deleteMeshset(entry.first, entry.second);
  Range tets, skin;
  CHKERR moab.get_entities_by_type(0, MBTET, tets);
  if (tets.empty())
    SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
            "Patch fixture requires tetrahedra");
  Skinner skinner(&moab);
  CHKERR skinner.find_skin(0, tets, false, skin);
  const std::string name = "ANALYTICAL_DISPLACEMENT_1001";
  CHKERR meshsets->addMeshset(BLOCKSET, 1001, name);
  CHKERR meshsets->addEntitiesToMeshset(BLOCKSET, 1001, skin);
  CHKERR meshsets->setAttributes(BLOCKSET, 1001, {1., 1., 1.}, name);
  CHKERR moab.write_file(output);
  CHKERR PetscPrintf(
      PETSC_COMM_WORLD, "Prepared affine patch: %d tetrahedra, %d skin faces\n",
      static_cast<int>(tets.size()), static_cast<int>(skin.size()));
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode runPatch(EshelbianCore &ep) {
  MoFEMFunctionBegin;
  auto reference = boost::make_shared<Reference>();
  CHKERR ep.setElasticElementOps(1);
  CHKERR ep.setElasticElementToTs(ep.dmElastic);
  auto state = createDMVector(ep.dmElastic);
  CHKERR VecZeroEntries(state);
  auto ts = createTS(PETSC_COMM_WORLD);
  CHKERR TSSetType(ts, TSBEULER);
  CHKERR TSAdaptRegister(TSADAPTMOFEM, TSAdaptCreateMoFEM);
  TSAdapt adapt;
  CHKERR TSGetAdapt(ts, &adapt);
  CHKERR TSAdaptSetType(adapt, TSADAPTNONE);
  CHKERR ep.solveElastic(ts, state);
  TSConvergedReason reason;
  PetscReal time;
  CHKERR TSGetConvergedReason(ts, &reason);
  CHKERR TSGetTime(ts, &time);
  if (reason <= 0 || std::abs(time - 1.) > 1.e-12)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Patch solve did not reach t=1: reason %d, time %g", reason, time);
  CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, state, INSERT_VALUES,
                                 SCATTER_REVERSE);
 
  auto audit = boost::make_shared<Audit>();
  auto volume = boost::make_shared<VolumeElement>(ep.mField);
  CHKERR ep.setBaseVolumeElementOps(1, false, false, false, volume);
  volume->getOpPtrVector().push_back(
      new OpCheckVolume(ep.dataAtPts, audit, reference));
  CHKERR DMoFEMLoopFiniteElements(ep.dmElastic, ep.elementVolumeName, volume);
 
  auto face = boost::make_shared<FaceElement>(ep.mField);
  face->getRuleHook = [](int, int, int order) { return 2 * order + 1; };
  CHKERR EshelbianPlasticity::AddHOOps<2, 3, 3>::add(
      face->getOpPtrVector(), {L2}, ep.materialH1Positions, ep.frontAdjEdges);
  auto *side =
      new OpLoopSide<EleOnSide>(ep.mField, ep.elementVolumeName, 3, Sev::noisy);
  face->getOpPtrVector().push_back(side);
  auto side_fe = side->getSideFEPtr();
  side_fe->getUserPolynomialBase() =
      boost::make_shared<CGGUserPolynomialBase>();
  CHKERR EshelbianPlasticity::AddHOOps<3, 3, 3>::add(
      side_fe->getOpPtrVector(), {H1, HDIV, L2}, ep.materialH1Positions,
      ep.frontAdjEdges);
  side_fe->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
      ep.piolaStress, ep.dataAtPts->getApproxPAtPts(),
      boost::make_shared<double>(1.)));
  side_fe->getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
      ep.bubbleField, ep.dataAtPts->getApproxPAtPts(), MBMAXTYPE));
  side_fe->getOpPtrVector().push_back(
      new OpPatchTraction(ep.dataAtPts, audit, reference));
  face->getOpPtrVector().push_back(new OpCheckFace(ep.dataAtPts, audit));
  CHKERR DMoFEMLoopFiniteElements(ep.dmElastic, ep.skinElement, face);
 
  std::array<double, Audit::COUNT> error;
  std::array<double, 3> reaction;
  PetscInt points, faces;
  double area;
  CHKERR MPI_Allreduce(audit->error.data(), error.data(), Audit::COUNT,
                       MPI_DOUBLE, MPI_MAX, ep.mField.get_comm());
  CHKERR MPI_Allreduce(audit->reaction.data(), reaction.data(), 3, MPI_DOUBLE,
                       MPI_SUM, ep.mField.get_comm());
  CHKERR MPI_Allreduce(&audit->volumePoints, &points, 1, MPIU_INT, MPI_SUM,
                       ep.mField.get_comm());
  CHKERR MPI_Allreduce(&audit->faces, &faces, 1, MPIU_INT, MPI_SUM,
                       ep.mField.get_comm());
  CHKERR MPI_Allreduce(&audit->boundaryArea, &area, 1, MPI_DOUBLE, MPI_SUM,
                       ep.mField.get_comm());
  double tolerance = 1.e-7;
  CHKERR PetscOptionsGetReal(nullptr, nullptr, "-patch_tolerance", &tolerance,
                             nullptr);
  if (!(tolerance > 0.) || points == 0 || faces == 0 || !(area > 0.))
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Invalid patch tolerance or incomplete quadrature coverage");
  const std::array<const char *, Audit::COUNT> names{
      "F",     "P",           "D",
      "theta", "Td",          "R",
      "w",     "normal flux", "integrated face flux/area"};
  for (int index = 0; index != Audit::COUNT; ++index) {
    CHKERR PetscPrintf(ep.mField.get_comm(), "Patch %s error %.3e\n",
                       names[index], error[index]);
    if (!std::isfinite(error[index]) || error[index] > tolerance)
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Patch %s error %g exceeds %g", names[index], error[index],
              tolerance);
  }
  const FTensor::Tensor1<double, 3> t_reaction{reaction[0], reaction[1],
                                               reaction[2]};
  if (!std::isfinite(t_reaction.l2()) || t_reaction.l2() / area > tolerance)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Patch total flux/area %g exceeds %g", t_reaction.l2() / area,
            tolerance);
  CHKERR PetscPrintf(ep.mField.get_comm(),
                     "Auxiliary homogeneous patch passed: %d volume points, "
                     "%d faces, total flux/area %.3e\n",
                     static_cast<int>(points), static_cast<int>(faces),
                     t_reaction.l2() / area);
  MoFEMFunctionReturn(0);
}
 
} // namespace
 
static char help[] =
    "Prepare and solve an affine auxiliary logarithmic patch.\n";
 
int main(int argc, char *argv[]) {
  MoFEM::Core::Initialize(&argc, &argv, nullptr, help);
#ifdef ENABLE_PYTHON_BINDING
  Py_Initialize();
  np::initialize();
#endif
  auto core_log = logging::core::get();
  core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(), "EP"));
  LogManager::setLog("EP");
  core_log->add_sink(
      LogManager::createSink(LogManager::getStrmSelf(), "EPSELF"));
  LogManager::setLog("EPSELF");
  core_log->add_sink(
      LogManager::createSink(LogManager::getStrmSync(), "EPSYNC"));
  LogManager::setLog("EPSYNC");
  try {
    CHKERR JsonConfigManager::setMeshFileFromJson();
    // Read the mode before loading its mesh. The normal test setup creates its
    // own MoFEM database and applies the same JSON configuration there.
    {
      moab::Core options_moab;
      MoFEM::Core options_core(options_moab);
      options_core.getInterface<JsonConfigManager>();
    }
    PetscBool prepare = PETSC_FALSE;
    CHKERR PetscOptionsGetBool(nullptr, nullptr, "-patch_prepare", &prepare,
                               nullptr);
    if (prepare)
      CHKERR prepareMesh();
    else
      CHKERR runAuxiliaryLogarithmicStressTest(runPatch);
  }
  CATCH_ERRORS;
  MoFEM::Core::Finalize();
  return 0;
}