HMHPMooneyRivlinWriggersEq63.cpp

Go to the documentation of this file.

/**
 * @file HMHPMooneyRivlinWriggersEq63.cpp
 * @brief Implement of MooneyRivlinWriggersEq63
 *
 * @copyright Copyright (c) 2024
 *
 */
 
namespace EshelbianPlasticity {
 
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;
    PetscOptionsBegin(PETSC_COMM_WORLD, "mooneyrivlin_", "", "none");
    alpha = 1;
    CHKERR PetscOptionsScalar("-alpha", "Alpha", "", alpha, &alpha, PETSC_NULLPTR);
    beta = 1;
    CHKERR PetscOptionsScalar("-beta", "Beta", "", beta, &beta, PETSC_NULLPTR);
 
    lambda = 1;
    CHKERR PetscOptionsScalar("-lambda", "Lambda", "", lambda, &lambda,
                              PETSC_NULLPTR);
 
    epsilon = 0;
    CHKERR PetscOptionsScalar("-epsilon", "Epsilon", "", epsilon, &epsilon,
                              PETSC_NULLPTR);
 
    PetscOptionsEnd();
    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;
  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();
    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);
  }
};
 
 
 
} // namespace EshelbianPlasticity