auxiliary_logarithmic_stress_layout_atom.cpp

Go to the documentation of this file.

/**
 * @file auxiliary_logarithmic_stress_layout_atom.cpp
 * @brief Verify auxiliary material fields, distributed layouts and evaluation.
 */
 
#include <MoFEM.hpp>
using namespace MoFEM;
#include "AuxiliaryLogarithmicStressTestSetup.hpp"
#include <AuxiliaryLogarithmicStress.hpp>
using namespace EshelbianPlasticity;
 
namespace {
 
using Coordinates = FTensor::Tensor1<double, 5>;
using VolumeElement = VolumeElementForcesAndSourcesCore;
using MatrixPtr = boost::shared_ptr<MatrixDouble>;
 
MoFEMErrorCode checkNear(const double value, const double expected,
                         const char *description) {
  MoFEMFunctionBegin;
  if (!std::isfinite(value) ||
      std::abs(value - expected) > 1.e-11 * std::max(1., std::abs(expected)))
    SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
            "%s: expected %.16e, obtained %.16e", description, expected, value);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode getCellScale(moab::Interface &moab, const EntityHandle entity,
                            double &scale) {
  MoFEMFunctionBegin;
  const EntityHandle *connectivity;
  int count;
  CHKERR moab.get_connectivity(entity, connectivity, count, true);
  std::vector<double> coordinates(3 * count);
  CHKERR moab.get_coords(connectivity, count, coordinates.data());
  scale = 1.;
  for (const double coordinate : coordinates)
    scale += 0.001 * coordinate / count;
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode checkFieldLayout(EshelbianCore &ep) {
  MoFEMFunctionBegin;
  const auto expected_features = PhysicalEquations::Features{}.set(
      PhysicalEquations::AUXILIARY_LOGARITHMIC_STRESS);
  if (ep.physicalEquations->getFeatures() != expected_features ||
      ep.mField.check_field(ep.stretchTensor))
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Auxiliary mode must register D/theta/Td without the old u field");
  const std::vector<std::string> names{ep.logDeviator, ep.logJacobian,
                                       ep.auxiliaryLogStress};
  if (ep.physicalEquations->getMaterialFields(ep) != names)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Unexpected auxiliary material-field group");
  const auto *finite_element =
      ep.mField.get_finite_element_structure(ep.elementVolumeName);
  const int order =
      ep.spaceOrder + (ep.brokenHdivBase == AINSWORTH_LEGENDRE_BASE ? 1 : 0);
  const PetscInt basis_count = (order + 1) * (order + 2) * (order + 3) / 6;
  const Problem *problem;
  CHKERR DMMoFEMGetProblemPtr(ep.dmElastic, &problem);
  int rank;
  CHKERR MPI_Comm_rank(ep.mField.get_comm(), &rank);
  // Count active owned volume elements, excluding broken-space mesh copies.
  PetscInt local_tets = 0, global_tets;
  for (const auto &element : *problem->getNumeredFiniteElementsPtr())
    if (element->getName() == ep.elementVolumeName &&
        element->getPart() == rank && element->getEntType() == MBTET)
      ++local_tets;
  CHKERR MPI_Allreduce(&local_tets, &global_tets, 1, MPIU_INT, MPI_SUM,
                       ep.mField.get_comm());
  if (!global_tets)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "Layout atom requires tetrahedra");
  for (const auto &name : names) {
    const int components = name == ep.logJacobian ? 1 : 5;
    const auto *field = ep.mField.get_field_structure(name);
    if (field->getNbOfCoeffs() != components || field->getSpace() != L2 ||
        field->getApproxBase() != USER_BASE)
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Field %s must have %d USER_BASE L2 components", name.c_str(),
              components);
    for (const auto mask : {finite_element->getBitFieldIdRow(),
                            finite_element->getBitFieldIdCol(),
                            finite_element->getBitFieldIdData()})
      if ((mask & field->getId()) != field->getId())
        SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
                "Field %s missing from volume FE row, column or data",
                name.c_str());
    for (const auto entry :
         {std::make_pair(ep.dmElastic.get(), RowColData::ROW),
          std::make_pair(ep.dmElastic.get(), RowColData::COL),
          std::make_pair(ep.dmMaterial.get(), RowColData::COL)}) {
      IS raw_is;
      CHKERR DMMoFEMGetFieldIS(entry.first, entry.second, name.c_str(),
                               &raw_is);
      SmartPetscObj<IS> field_is(raw_is);
      PetscInt size;
      CHKERR ISGetSize(field_is, &size);
      if (size != global_tets * components * basis_count)
        SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
                "Field %s in %s/%s has %d DOFs, expected %d", name.c_str(),
                entry.first == ep.dmElastic.get() ? "dmElastic" : "dmMaterial",
                entry.second == RowColData::ROW ? "ROW" : "COL",
                static_cast<int>(size),
                static_cast<int>(global_tets * components * basis_count));
    }
  }
  for (const auto &dof : *problem->getNumeredRowDofsPtr()) {
    if (std::find(names.begin(), names.end(), dof->getName()) == names.end())
      continue;
    if (dof->getEntType() != MBTET ||
        dof->getFieldEntityPtr()->getMaxOrder() != order)
      SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
              "Material DOF has unexpected entity or approximation order");
  }
  CHKERR PetscPrintf(ep.mField.get_comm(),
                     "Auxiliary layout: %d tetrahedra, order %d, "
                     "%d scalar basis functions per field\n",
                     static_cast<int>(global_tets), order,
                     static_cast<int>(basis_count));
  MoFEMFunctionReturn(0);
}
 
double getCoefficient(const EshelbianCore &ep, const std::string &name,
                      const int component) {
  if (name == ep.logDeviator)
    return 0.02 * (component + 1);
  if (name == ep.logJacobian)
    return 0.18;
  return -0.07 * (component + 1);
}
 
MoFEMErrorCode setMaterialVector(EshelbianCore &ep, Vec vector,
                                 const double state_scale) {
  MoFEMFunctionBegin;
  CHKERR VecZeroEntries(vector);
  PetscInt first, last;
  CHKERR VecGetOwnershipRange(vector, &first, &last);
  const Problem *problem;
  CHKERR DMMoFEMGetProblemPtr(ep.dmElastic, &problem);
  const auto names = ep.physicalEquations->getMaterialFields(ep);
  for (const auto &dof : *problem->getNumeredRowDofsPtr()) {
    const auto index = dof->getPetscGlobalDofIdx();
    if (index < first || index >= last || dof->getDofOrder() != 0 ||
        std::find(names.begin(), names.end(), dof->getName()) == names.end())
      continue;
    double cell_scale;
    CHKERR getCellScale(ep.mField.get_moab(), dof->getEnt(), cell_scale);
    CHKERR VecSetValue(
        vector, index,
        state_scale * cell_scale *
            getCoefficient(ep, dof->getName(), dof->getDofCoeffIdx()),
        INSERT_VALUES);
  }
  CHKERR VecAssemblyBegin(vector);
  CHKERR VecAssemblyEnd(vector);
  CHKERR VecGhostUpdateBegin(vector, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(vector, INSERT_VALUES, SCATTER_FORWARD);
  MoFEMFunctionReturn(0);
}
 
struct OpCheckMaterialFields : public VolumeElement::UserDataOperator {
  OpCheckMaterialFields(EshelbianCore &ep,
                        boost::shared_ptr<AuxiliaryLogarithmicStressData> data,
                        MatrixPtr reconstructed, const double state_scale,
                        boost::shared_ptr<PetscInt> visited,
                        const bool compatible_state)
      : VolumeElement::UserDataOperator(NOSPACE, OPSPACE), eP(ep),
        dataPtr(std::move(data)), reconstructedPtr(std::move(reconstructed)),
        stateScale(state_scale), visitedPtr(std::move(visited)),
        compatibleState(compatible_state) {}
 
  MoFEMErrorCode doWork(int, EntityType, EntData &) override {
    MoFEMFunctionBegin;
    FTensor::Index<'A', 5> A;
    FTensor::Index<'i', 3> i;
    FTensor::Index<'j', 3> j;
    const int points = getGaussPts().size2();
    if (points <= 1)
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "Layout test requires multiple integration points");
    auto t_d = MatrixSizeHelper<GetFTensor1FromMatType<5, -1, DL>, DL>::get(
        *dataPtr->logDeviator, points)();
    auto t_theta = MatrixSizeHelper<GetFTensor1FromMatType<1, -1, DL>, DL>::get(
        *dataPtr->logJacobian, points)();
    auto t_stress =
        MatrixSizeHelper<GetFTensor1FromMatType<5, -1, DL>, DL>::get(
            *dataPtr->stress, points)();
    auto t_h =
        MatrixSizeHelper<GetFTensor2SymmetricFromMatType<3, -1, DL>, DL>::get(
            *reconstructedPtr, points)();
    double cell_scale = 1.;
    if (!compatibleState)
      CHKERR getCellScale(eP.mField.get_moab(), getFEEntityHandle(),
                          cell_scale);
    const double scale = stateScale * cell_scale;
    Coordinates t_expected_d, t_expected_stress;
    for (int aa = 0; aa != 5; ++aa) {
      t_expected_d(aa) = scale * getCoefficient(eP, eP.logDeviator, aa);
      t_expected_stress(aa) =
          scale * getCoefficient(eP, eP.auxiliaryLogStress, aa);
    }
    if (compatibleState) {
      double energy;
      CHKERR NeoHookeanLogarithmicMaterial::evaluateDeviator(
          {1.7, 8.5}, t_expected_d, energy, &t_expected_stress);
    }
    for (int gg = 0; gg != points; ++gg) {
      Coordinates t_error;
      t_error(A) = t_d(A) - t_expected_d(A);
      CHKERR checkNear(std::sqrt(t_error(A) * t_error(A)), 0., "evaluated D");
      t_error(A) = t_stress(A) - t_expected_stress(A);
      CHKERR checkNear(std::sqrt(t_error(A) * t_error(A)), 0., "evaluated Td");
      CHKERR checkNear(t_theta(0), 0.18 * scale, "evaluated theta");
      CHKERR checkNear(t_h(i, i), t_theta(0), "reconstructed trace");
      CHKERR checkNear(t_h(i, j) * t_h(i, j),
                       t_d(A) * t_d(A) + t_theta(0) * t_theta(0) / 3.,
                       "reconstructed Frobenius norm");
      const auto t_recovered = Tensor2SymmetricDeviatorBasis::getCoordinates(t_h);
      t_error(A) = t_recovered(A) - t_d(A);
      CHKERR checkNear(std::sqrt(t_error(A) * t_error(A)), 0.,
                       "reconstructed deviator");
      if (compatibleState) {
        Coordinates t_material_stress;
        t_material_stress(A) = t_stress(A);
        NeoHookeanLogarithmicMaterial::InverseState inverse;
        CHKERR NeoHookeanLogarithmicMaterial::evaluateInverse(
            {1.7, 8.5}, t_material_stress, inverse);
        t_error(A) = inverse.tMaterialDeviator(A) - t_d(A);
        CHKERR checkNear(std::sqrt(t_error(A) * t_error(A)), 0.,
                         "constant initialized material copy");
      }
      ++t_d;
      ++t_theta;
      ++t_stress;
      ++t_h;
    }
    ++*visitedPtr;
    MoFEMFunctionReturn(0);
  }
 
  EshelbianCore &eP;
  boost::shared_ptr<AuxiliaryLogarithmicStressData> dataPtr;
  MatrixPtr reconstructedPtr;
  double stateScale;
  boost::shared_ptr<PetscInt> visitedPtr;
  bool compatibleState;
};
 
MoFEMErrorCode checkStateEvaluation(EshelbianCore &ep,
                                    SmartPetscObj<Vec> current,
                                    SmartPetscObj<Vec> previous,
                                    const double current_scale,
                                    const double previous_scale,
                                    const bool compatible_state = false) {
  MoFEMFunctionBegin;
  if (compatible_state) {
    CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, current, INSERT_VALUES,
                                   SCATTER_FORWARD);
    CHKERR VecCopy(current, previous);
    CHKERR VecGhostUpdateBegin(previous, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(previous, INSERT_VALUES, SCATTER_FORWARD);
  } else {
    CHKERR setMaterialVector(ep, current, current_scale);
    CHKERR setMaterialVector(ep, previous, previous_scale);
  }
  CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, current, INSERT_VALUES,
                                 SCATTER_REVERSE);
  auto round_trip = createDMVector(ep.dmElastic);
  CHKERR VecZeroEntries(round_trip);
  CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, round_trip, INSERT_VALUES,
                                 SCATTER_FORWARD);
  CHKERR VecAXPY(round_trip, -1., current);
  PetscReal error;
  CHKERR VecNorm(round_trip, NORM_INFINITY, &error);
  CHKERR checkNear(error, 0., "current vector-mesh-vector transfer");
 
  auto fe = boost::make_shared<VolumeElement>(ep.mField);
  fe->getUserPolynomialBase() =
      boost::make_shared<CGGUserPolynomialBase>(nullptr, true);
  fe->getRuleHook = [](int, int, int) { return 4; };
  MoFEM::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      fe->getOpPtrVector(), {L2}, ep.materialH1Positions);
  auto visited = boost::make_shared<PetscInt>(0);
  for (const bool previous_state : {false, true}) {
    auto data = boost::make_shared<AuxiliaryLogarithmicStressData>();
    auto reconstructed = boost::make_shared<MatrixDouble>();
    CHKERR pushAuxiliaryLogarithmicFields(
        ep, fe->getOpPtrVector(), data, reconstructed,
        previous_state ? previous : SmartPetscObj<Vec>());
    fe->getOpPtrVector().push_back(new OpCheckMaterialFields(
        ep, data, reconstructed,
        previous_state ? previous_scale : current_scale, visited,
        compatible_state));
  }
  CHKERR DMoFEMLoopFiniteElements(ep.dmElastic, ep.elementVolumeName, fe);
  PetscInt global_visited;
  CHKERR MPI_Allreduce(visited.get(), &global_visited, 1, MPIU_INT, MPI_SUM,
                       ep.mField.get_comm());
  if (!global_visited)
    SETERRQ(ep.mField.get_comm(), MOFEM_ATOM_TEST_INVALID,
            "No auxiliary field evaluations were executed");
  MoFEMFunctionReturn(0);
}
 
/** Set the fixture's constant material state, including local ghost copies. */
MoFEMErrorCode setConstantMaterialState(EshelbianCore &ep,
                                      const Coordinates &t_deviator,
                                      const double theta) {
  MoFEMFunctionBegin;
  Coordinates t_stress;
  double energy;
  CHKERR NeoHookeanLogarithmicMaterial::evaluateDeviator(
      {1.7, 8.5}, t_deviator, energy, &t_stress);
  auto field_blas = ep.mField.getInterface<FieldBlas>();
  auto set_field = [&]<int Dim>(const std::string &field,
                               const FTensor::Tensor1<double, Dim> &t_values) {
    MoFEMFunctionBegin;
    CHKERR field_blas->setField(0., field);
    auto set_constant = [&](boost::shared_ptr<FieldEntity> entity_ptr) {
      MoFEMFunctionBegin;
      auto data = entity_ptr->getEntFieldData();
      if (data.size() < Dim)
        SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                "Field %s has no constant mode", field.c_str());
      FTensor::Index<'i', Dim> i;
      auto t_dof = getFTensor1FromPtr<Dim>(&data[0]);
      t_dof(i) = t_values(i);
      MoFEMFunctionReturn(0);
    };
    CHKERR field_blas->fieldLambdaOnEntities(set_constant, field);
    MoFEMFunctionReturn(0);
  };
  CHKERR set_field(ep.logDeviator, t_deviator);
  const FTensor::Tensor1<double, 1> t_theta(theta);
  CHKERR set_field(ep.logJacobian, t_theta);
  CHKERR set_field(ep.auxiliaryLogStress, t_stress);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode runLayoutAtom(EshelbianCore &ep) {
  MoFEMFunctionBegin;
  CHKERR checkFieldLayout(ep);
  auto current = createDMVector(ep.dmElastic);
  auto previous = createDMVector(ep.dmElastic);
  FTensor::Index<'A', 5> A;
  Coordinates t_initial;
  t_initial(A) = 0.;
  CHKERR setConstantMaterialState(ep, t_initial, 0.);
  CHKERR checkStateEvaluation(ep, current, previous, 0., 0., true);
  for (const auto &name : ep.physicalEquations->getMaterialFields(ep))
    CHKERR ep.mField.getInterface<FieldBlas>()->setField(0.4, name);
  for (int aa = 0; aa != 5; ++aa)
    t_initial(aa) = getCoefficient(ep, ep.logDeviator, aa);
  CHKERR setConstantMaterialState(ep, t_initial, 0.18);
  CHKERR checkStateEvaluation(ep, current, previous, 1., 1., true);
  CHKERR checkStateEvaluation(ep, current, previous, 1., -0.5);
  CHKERR checkStateEvaluation(ep, current, previous, -0.25, 1.5);
  CHKERR PetscPrintf(PETSC_COMM_WORLD,
                     "Auxiliary logarithmic stress layout atom passed\n");
  MoFEMFunctionReturn(0);
}
 
} // namespace
 
static char help[] = "Verify the auxiliary logarithmic stress field layout.\n";
 
int main(int argc, char *argv[]) {
  MoFEM::Core::Initialize(&argc, &argv, nullptr, help);
  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 runAuxiliaryLogarithmicStressTest(runLayoutAtom);
  }
  CATCH_ERRORS;
  MoFEM::Core::Finalize();
  return 0;
}