EshelbianADOL-C.cpp

Go to the documentation of this file.

/**
 * \file EshelbianPlasticity.cpp
 * \brief Implementation of automatic differentiation
 */
 
#include <MoFEM.hpp>
using namespace MoFEM;
 
#include <BasicFiniteElements.hpp>
 
#include <EshelbianPlasticity.hpp>
#include <boost/math/constants/constants.hpp>
 
constexpr double third = boost::math::constants::third<double>();
 
namespace EshelbianPlasticity {
 
struct OpHMHH : public OpJacobian {
  OpHMHH(const int tag, const bool eval_rhs, const bool eval_lhs,
         boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
         boost::shared_ptr<PhysicalEquations> &physics_ptr)
      : OpJacobian(tag, eval_rhs, eval_lhs, data_ptr, physics_ptr) {}
 
  MoFEMErrorCode evaluateRhs(EntData &data);
  MoFEMErrorCode evaluateLhs(EntData &data);
};
 
MoFEMErrorCode OpHMHH::evaluateRhs(EntData &ent_data) {
  MoFEMFunctionBegin;
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'l', 3> l;
  FTensor::Index<'m', 3> m;
  FTensor::Index<'n', 3> n;
 
  const auto nb_integration_pts = getGaussPts().size2();
  auto iu = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
  auto t_grad_h1 = getFTensor2FromMat<3, 3>(dataAtPts->wGradH1AtPts);
 
  auto create_data_vec = [nb_integration_pts](auto &v) {
    v.resize(nb_integration_pts, false);
  };
  auto create_data_mat = [nb_integration_pts](auto &m) {
    m.resize(9, nb_integration_pts, false);
  };
 
  create_data_mat(dataAtPts->PAtPts);
 
  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
  auto r_P = getFTensor2FromMat<3, 3>(dataAtPts->PAtPts);
  for (unsigned int gg = 0; gg != nb_integration_pts; ++gg) {
 
    FTensor::Tensor2<double, 3, 3> t_h1;
 
    switch (EshelbianCore::gradApperoximator) {
    case LARGE_ROT:
    case MODERATE_ROT:
      t_h1(i, j) = t_grad_h1(i, j) + t_kd(i, j);
      physicsPtr->get_h()(i, j) = iu(i, k) * t_h1(k, j);
      break;
    case SMALL_ROT:
      physicsPtr->get_h()(i, j) = iu(i, j);
      break;
    };
 
    int r = ::function(tAg, physicsPtr->dependentVariablesPiola.size(),
                       physicsPtr->activeVariables.size(),
                       &physicsPtr->activeVariables[0],
                       &physicsPtr->dependentVariablesPiola[0]);
    if (r < 0) { // function is locally analytic
      SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
               "ADOL-C function evaluation with error r = %d", r);
    }
 
    switch (EshelbianCore::gradApperoximator) {
    case LARGE_ROT:
    case MODERATE_ROT: {
      FTensor::Tensor4<double, 3, 3, 3, 3> t_h1_du;
      t_h1_du(i, j, m, n) = t_kd(i, m) * (t_kd(k, n) * t_h1(k, j));
      r_P(k, l) = physicsPtr->get_P()(i, j) * t_h1_du(i, j, k, l);
    }; break;
    case SMALL_ROT:
      r_P(i, j) = physicsPtr->get_P()(i, j);
      break;
    };
 
    ++iu;
    ++t_grad_h1;
    ++r_P;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode OpHMHH::evaluateLhs(EntData &ent_data) {
  MoFEMFunctionBegin;
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'l', 3> l;
  FTensor::Index<'m', 3> m;
  FTensor::Index<'n', 3> n;
  FTensor::Index<'o', 3> o;
  FTensor::Index<'p', 3> p;
 
  const int number_of_active_variables = physicsPtr->activeVariables.size();
  const int number_of_dependent_variables =
      physicsPtr->dependentVariablesPiola.size();
  std::vector<double *> jac_ptr(number_of_dependent_variables);
  for (unsigned int n = 0; n != number_of_dependent_variables; ++n) {
    jac_ptr[n] =
        &(physicsPtr
              ->dependentVariablesPiolaDirevatives[n *
                                                   number_of_active_variables]);
  }
 
  const auto nb_integration_pts = getGaussPts().size2();
 
  auto create_data_mat = [nb_integration_pts](auto &m) {
    m.resize(9, nb_integration_pts, false);
  };
 
  dataAtPts->P_du.resize(81, nb_integration_pts, false);
 
  auto iu = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
  auto t_grad_h1 = getFTensor2FromMat<3, 3>(dataAtPts->wGradH1AtPts);
  auto r_P_du = getFTensor4FromMat<3, 3, 3, 3>(dataAtPts->P_du);
 
  FTensor::Number<0> N0;
  FTensor::Number<1> N1;
  FTensor::Number<2> N2;
 
  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
 
  for (unsigned int gg = 0; gg != nb_integration_pts; ++gg) {
 
    FTensor::Tensor2<double, 3, 3> t_h1;
 
    switch (EshelbianCore::gradApperoximator) {
    case LARGE_ROT:
    case MODERATE_ROT:
      t_h1(i, j) = t_grad_h1(i, j) + t_kd(i, j);
      physicsPtr->get_h()(i, j) = iu(i, k) * t_h1(k, j);
      break;
    case SMALL_ROT:
      physicsPtr->get_h()(i, j) = iu(i, j);
      break;
    };
 
    // play recorder for jacobians
    int r = ::jacobian(tAg, number_of_dependent_variables,
                       number_of_active_variables,
                       &physicsPtr->activeVariables[0], &jac_ptr[0]);
    if (r < 0) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
              "ADOL-C function evaluation with error");
    }
 
    FTensor::Tensor4<double, 3, 3, 3, 3> t_P_dh_tmp;
    t_P_dh_tmp(i, j, N0, k) = physicsPtr->get_P_dh0()(i, j, k);
    t_P_dh_tmp(i, j, N1, k) = physicsPtr->get_P_dh1()(i, j, k);
    t_P_dh_tmp(i, j, N2, k) = physicsPtr->get_P_dh2()(i, j, k);
 
    switch (EshelbianCore::gradApperoximator) {
    case LARGE_ROT:
    case MODERATE_ROT: {
      FTensor::Tensor4<double, 3, 3, 3, 3> t_h1_du;
      t_h1_du(i, j, m, n) = t_kd(i, m) * (t_kd(k, n) * t_h1(k, j));
      r_P_du(k, l, m, n) =
          (t_P_dh_tmp(i, j, o, p) * t_h1_du(o, p, m, n)) * t_h1_du(i, j, k, l);
    } break;
    case SMALL_ROT:
      r_P_du(i, j, m, n) = t_P_dh_tmp(i, j, m, n);
      break;
    };
 
    ++iu;
    ++t_grad_h1;
    ++r_P_du;
  }
  MoFEMFunctionReturn(0);
}
 
struct HMHStVenantKirchhoff : public PhysicalEquations {
 
  static constexpr int numberOfActiveVariables = 9;
  static constexpr int numberOfDependentVariables = 9;
 
  HMHStVenantKirchhoff(const double lambda, const double mu)
      : PhysicalEquations(numberOfActiveVariables, numberOfDependentVariables),
        lambda(lambda), mu(mu) {}
 
  MoFEMErrorCode getOptions() {
    MoFEMFunctionBegin;
 
    double E = 1;
    double nu = 0;
 
    CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "stvenant_", "", "none");
 
    CHKERR PetscOptionsScalar("-young_modulus", "Young modulus", "", E, &E,
                              PETSC_NULL);
    CHKERR PetscOptionsScalar("-poisson_ratio", "poisson ratio", "", nu, &nu,
                              PETSC_NULL);
 
    ierr = PetscOptionsEnd();
    CHKERRG(ierr);
 
    MOFEM_LOG("EP", Sev::inform) << "St Venant Kirchhoff model parameters: "
                                 << "E = " << E << ", nu = " << nu;
 
    lambda = LAMBDA(E, nu);
    mu = MU(E, nu);
 
    MoFEMFunctionReturn(0);
  }
 
  virtual OpJacobian *
  returnOpJacobian(const int tag, const bool eval_rhs, const bool eval_lhs,
                   boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
                   boost::shared_ptr<PhysicalEquations> &physics_ptr) {
    return (new OpHMHH(tag, eval_rhs, eval_lhs, data_ptr, physics_ptr));
  }
 
  double lambda;
  double mu;
  double sigmaY;
 
  ATensor2 th;
  ATensor2 tH;
  ATensor2 tP;
  ATensor2 tSigma;
  adouble phi;
  adouble s2;
  adouble f;
  adouble energy;
 
  adouble detH;
  adouble detF;
  adouble trE;
  ATensor2 tInvH;
  ATensor2 tF;
  ATensor2 tInvF;
  ATensor2 tC;
  ATensor2 tE;
  ATensor2 tS;
 
  ATensor2 tPulledP;
 
  MoFEMErrorCode recordTape(const int tape, DTensor2Ptr *t_h_ptr) {
    MoFEMFunctionBegin;
    FTensor::Index<'i', 3> i;
    FTensor::Index<'j', 3> j;
    FTensor::Index<'I', 3> I;
    FTensor::Index<'J', 3> J;
    FTensor::Number<0> N0;
    FTensor::Number<1> N1;
    FTensor::Number<2> N2;
 
    CHKERR getOptions();
 
    auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
 
    auto ih = get_h();
    // auto iH = get_H();
    if (t_h_ptr)
      ih(i, j) = (*t_h_ptr)(i, j);
    else {
      ih(i, j) = t_kd(i, j);
    }
 
    auto r_P = get_P();
 
    enableMinMaxUsingAbs();
    trace_on(tape);
 
    // Set active variables to ADOL-C
    th(i, j) <<= get_h()(i, j);
 
    tF(i, I) = th(i, I);
    CHKERR determinantTensor3by3(tF, detF);
    CHKERR invertTensor3by3(tF, detF, tInvF);
 
    // Deformation and strain
    tC(I, J) = tF(i, I) * tF(i, J);
    tE(I, J) = tC(I, J);
    tE(N0, N0) -= 1;
    tE(N1, N1) -= 1;
    tE(N2, N2) -= 1;
    tE(I, J) *= 0.5;
 
    // Energy
    trE = tE(I, I);
    energy = 0.5 * lambda * trE * trE + mu * (tE(I, J) * tE(I, J));
 
    // Stress Piola II
    tS(I, J) = tE(I, J);
    tS(I, J) *= 2 * mu;
    tS(N0, N0) += lambda * trE;
    tS(N1, N1) += lambda * trE;
    tS(N2, N2) += lambda * trE;
    // Stress Piola I
    tP(i, J) = tF(i, I) * tS(I, J);
 
    // Set dependent variables to ADOL-C
    tP(i, j) >>= r_P(i, j);
 
    trace_off();
 
    MoFEMFunctionReturn(0);
  }
};
 
struct HMHPMooneyRivlinWriggersEq63 : public PhysicalEquations {
 
  static constexpr int numberOfActiveVariables = 9;
  static constexpr int numberOfDependentVariables = 9;
 
  HMHPMooneyRivlinWriggersEq63(const double alpha, const double beta,
                               const double lambda)
      : PhysicalEquations(numberOfActiveVariables, numberOfDependentVariables),
        alpha(alpha), beta(beta), lambda(lambda), epsilon(0) {}
 
  virtual OpJacobian *
  returnOpJacobian(const int tag, const bool eval_rhs, const bool eval_lhs,
                   boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
                   boost::shared_ptr<PhysicalEquations> &physics_ptr) {
    return (new OpHMHH(tag, eval_rhs, eval_lhs, data_ptr, physics_ptr));
  }
 
  MoFEMErrorCode getOptions() {
    MoFEMFunctionBegin;
    CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "mooneyrivlin_", "", "none");
 
    alpha = 1;
    CHKERR PetscOptionsScalar("-alpha", "Alpha", "", alpha, &alpha, PETSC_NULL);
    beta = 1;
    CHKERR PetscOptionsScalar("-beta", "Beta", "", beta, &beta, PETSC_NULL);
 
    lambda = 1;
    CHKERR PetscOptionsScalar("-lambda", "Lambda", "", lambda, &lambda,
                              PETSC_NULL);
 
    epsilon = 0;
    CHKERR PetscOptionsScalar("-epsilon", "Epsilon", "", epsilon, &epsilon,
                              PETSC_NULL);
 
    ierr = PetscOptionsEnd();
    CHKERRG(ierr);
    MoFEMFunctionReturn(0);
  }
 
  double alpha;
  double beta;
  double lambda;
  double epsilon;
 
  ATensor2 th;
  ATensor2 tH;
  ATensor2 tF;
 
  adouble detH;
  adouble detF;
  ATensor2 tInvH;
  ATensor2 tInvF;
 
  ATensor2 tP;
  ATensor2 tSigma;
 
  ATensor2 tCof;
  // ATensor3 tCofT1;
  // ATensor3 tCofT2;
  ATensor2 tBF;
  ATensor2 tBCof;
  adouble tBj;
 
  adouble energy;
  adouble phi;
  adouble A;
  adouble B;
 
  ATensor2 tPulledP;
 
  MoFEMErrorCode recordTape(const int tape, DTensor2Ptr *t_h_ptr) {
    MoFEMFunctionBegin;
 
    FTensor::Index<'i', 3> i;
    FTensor::Index<'j', 3> j;
    FTensor::Index<'k', 3> k;
    FTensor::Index<'I', 3> I;
    FTensor::Index<'J', 3> J;
    FTensor::Index<'K', 3> K;
    FTensor::Number<0> N0;
    FTensor::Number<1> N1;
    FTensor::Number<2> N2;
 
    CHKERR getOptions();
 
    auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
 
    auto ih = get_h();
    // auto iH = get_H();
    if (t_h_ptr)
      ih(i, j) = (*t_h_ptr)(i, j);
    else {
      ih(i, j) = t_kd(i, j);
    }
 
    auto r_P = get_P();
 
    enableMinMaxUsingAbs();
    trace_on(tape);
 
    // Set active variables to ADOL-C
    th(i, j) <<= get_h()(i, j);
 
    tF(i, I) = th(i, I);
    CHKERR determinantTensor3by3(tF, detF);
    CHKERR invertTensor3by3(tF, detF, tInvF);
 
    tCof(i, I) = detF * tInvF(I, i);
 
    A = tF(k, K) * tF(k, K);
    B = tCof(k, K) * tCof(k, K);
 
    tBF(i, I) = 4 * alpha * (A * tF(i, I));
    tBCof(i, I) = 4 * beta * (B * tCof(i, I));
    tBj = (-12 * alpha - 24 * beta) / detF +
          0.5 * (lambda / epsilon) *
              (pow(detF, epsilon - 1) - pow(detF, -epsilon - 1));
 
    tP(i, I) = tBF(i, I);
    tP(i, I) += (levi_civita(i, j, k) * tBCof(j, J)) *
                (levi_civita(I, J, K) * tF(k, K));
    tP(i, I) += tCof(i, I) * tBj;
 
    // Set dependent variables to ADOL-C
    tP(i, j) >>= r_P(i, j);
 
    trace_off();
 
    MoFEMFunctionReturn(0);
  }
};
 
MoFEMErrorCode EshelbianCore::addMaterial_HMHHStVenantKirchhoff(
    const int tape, const double lambda, const double mu,
    const double sigma_y) {
  MoFEMFunctionBegin;
  physicalEquations = boost::make_shared<HMHStVenantKirchhoff>(lambda, mu);
  CHKERR physicalEquations->recordTape(tape, nullptr);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::addMaterial_HMHMooneyRivlin(
    const int tape, const double alpha, const double beta, const double lambda,
    const double sigma_y) {
  MoFEMFunctionBegin;
  physicalEquations =
      boost::make_shared<HMHPMooneyRivlinWriggersEq63>(alpha, beta, lambda);
  CHKERR physicalEquations->recordTape(tape, nullptr);
  MoFEMFunctionReturn(0);
}
 
struct HMHHencky : public PhysicalEquations {
 
  HMHHencky(MoFEM::Interface &m_field, const double E, const double nu)
      : PhysicalEquations(0, 0), mField(m_field), E(E), nu(nu) {}
 
  MoFEMErrorCode recordTape(const int tag, DTensor2Ptr *t_h) { return 0; }
 
  struct OpHenckyJacobian : public OpJacobian {
    OpHenckyJacobian(boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                     boost::shared_ptr<HMHHencky> hencky_ptr)
        : OpJacobian(), dataAtGaussPts(data_ptr), henckyPtr(hencky_ptr) {
      CHK_THROW_MESSAGE(henckyPtr->getOptions(dataAtGaussPts),
                        "getOptions failed");
      CHK_THROW_MESSAGE(henckyPtr->extractBlockData(Sev::inform),
                        "Can not get data from block");
    }
 
    MoFEMErrorCode doWork(int side, EntityType type,
                          EntitiesFieldData::EntData &data) {
      MoFEMFunctionBegin;
 
      for (auto &b : henckyPtr->blockData) {
 
        if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
          CHKERR henckyPtr->getMatDPtr(dataAtGaussPts->getMatDPtr(),
                                       b.bulkModulusK, b.shearModulusG);
          MoFEMFunctionReturnHot(0);
        }
      }
 
      const auto E = henckyPtr->E;
      const auto nu = henckyPtr->nu;
 
      double bulk_modulus_K = E / (3 * (1 - 2 * nu));
      double shear_modulus_G = E / (2 * (1 + nu));
 
      CHKERR henckyPtr->getMatDPtr(dataAtGaussPts->getMatDPtr(), bulk_modulus_K,
                                   shear_modulus_G);
 
      MoFEMFunctionReturn(0);
    }
 
    MoFEMErrorCode evaluateRhs(EntData &data) { return 0; }
    MoFEMErrorCode evaluateLhs(EntData &data) { return 0; }
 
  private:
    boost::shared_ptr<DataAtIntegrationPts> dataAtGaussPts;
    boost::shared_ptr<HMHHencky> henckyPtr;
  };
 
  virtual OpJacobian *
  returnOpJacobian(const int tag, const bool eval_rhs, const bool eval_lhs,
                   boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
                   boost::shared_ptr<PhysicalEquations> &physics_ptr) {
    return (new OpHenckyJacobian(
        data_ptr, boost::dynamic_pointer_cast<HMHHencky>(physics_ptr)));
  }
 
  MoFEMErrorCode getOptions(boost::shared_ptr<DataAtIntegrationPts> data_ptr) {
    MoFEMFunctionBegin;
    CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "hencky_", "", "none");
 
    CHKERR PetscOptionsScalar("-young_modulus", "Young modulus", "", E, &E,
                              PETSC_NULL);
    CHKERR PetscOptionsScalar("-poisson_ratio", "poisson ratio", "", nu, &nu,
                              PETSC_NULL);
 
    ierr = PetscOptionsEnd();
    CHKERRG(ierr);
 
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode extractBlockData(Sev sev) {
    return extractBlockData(
 
        mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
 
            (boost::format("%s(.*)") % "MAT_ELASTIC").str()
 
                )),
 
        sev);
  }
 
  MoFEMErrorCode
  extractBlockData(std::vector<const CubitMeshSets *> meshset_vec_ptr,
                   Sev sev) {
    MoFEMFunctionBegin;
 
    for (auto m : meshset_vec_ptr) {
      MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock") << *m;
      std::vector<double> block_data;
      CHKERR m->getAttributes(block_data);
      if (block_data.size() != 2) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Expected that block has two attribute");
      }
      auto get_block_ents = [&]() {
        Range ents;
        CHKERR mField.get_moab().get_entities_by_handle(m->meshset, ents, true);
        return ents;
      };
 
      double young_modulus = block_data[0];
      double poisson_ratio = block_data[1];
      double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
      double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
 
      MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock")
          << "E = " << young_modulus << " nu = " << poisson_ratio;
 
      blockData.push_back({bulk_modulus_K, shear_modulus_G, get_block_ents()});
    }
    MOFEM_LOG_CHANNEL("WORLD");
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
                            double bulk_modulus_K, double shear_modulus_G) {
    MoFEMFunctionBegin;
    //! [Calculate elasticity tensor]
    auto set_material_stiffness = [&]() {
      FTensor::Index<'i', SPACE_DIM> i;
      FTensor::Index<'j', SPACE_DIM> j;
      FTensor::Index<'k', SPACE_DIM> k;
      FTensor::Index<'l', SPACE_DIM> l;
      constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
      double A = (SPACE_DIM == 2)
                     ? 2 * shear_modulus_G /
                           (bulk_modulus_K + (4. / 3.) * shear_modulus_G)
                     : 1;
      auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mat_D_ptr);
      t_D(i, j, k, l) = 2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
                        A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) *
                            t_kd(i, j) * t_kd(k, l);
    };
    //! [Calculate elasticity tensor]
    constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
    mat_D_ptr->resize(size_symm * size_symm, 1);
    set_material_stiffness();
    MoFEMFunctionReturn(0);
  }
 
private:
  MoFEM::Interface &mField;
 
  struct BlockData {
    double bulkModulusK;
    double shearModulusG;
    Range blockEnts;
  };
  std::vector<BlockData> blockData;
 
  double E;
  double nu;
};
 
MoFEMErrorCode EshelbianCore::addMaterial_Hencky(double E, double nu) {
  MoFEMFunctionBegin;
  physicalEquations = boost::make_shared<HMHHencky>(mField, E, nu);
  MoFEMFunctionReturn(0);
}
 
}; // namespace EshelbianPlasticity