HMHHStVenantKirchhoff.cpp

Go to the documentation of this file.

/**
 * @file HMHHStVenantKirchhoff.cpp
 * @brief Implementation StVenantKirchhoff
 * @version 0.1
 * @date 2024-08-31
 *
 * @copyright Copyright (c) 2024
 *
 */
 
namespace EshelbianPlasticity {
 
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;
 
    PetscOptionsBegin(PETSC_COMM_WORLD, "stvenant_", "", "none");
 
    CHKERR PetscOptionsScalar("-young_modulus", "Young modulus", "", E, &E,
                              PETSC_NULLPTR);
    CHKERR PetscOptionsScalar("-poisson_ratio", "poisson ratio", "", nu, &nu,
                              PETSC_NULLPTR);
 
    PetscOptionsEnd();
 
    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);
  }
};
 
 
 
} // namespace EshelbianPlasticity