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 <MatrixFunction.hpp>
 
#include <EshelbianPlasticity.hpp>
 
#include <EshelbianAux.hpp>
 
#include <boost/math/constants/constants.hpp>
 
constexpr double third = boost::math::constants::third<double>();
 
namespace EshelbianPlasticity {
 
struct OpSpatialPhysical : public OpAssembleVolume {
 
  OpSpatialPhysical(const std::string &field_name,
                    boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                    const double alpha_u)
      : OpAssembleVolume(field_name, data_ptr, OPROW), alphaU(alpha_u) {}
 
  MoFEMErrorCode integrate(EntData &data);
 
private:
  const double alphaU;
};
 
MoFEMErrorCode OpSpatialPhysical::integrate(EntData &data) {
  MoFEMFunctionBegin;
 
  FTensor::Index<'L', size_symm> L;
  auto t_L = symm_L_tensor();
 
  int nb_dofs = data.getIndices().size();
  int nb_integration_pts = data.getN().size1();
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
  auto t_approx_P_adjoint_log_du =
      getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
  auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->PAtPts);
  auto t_dot_log_u =
      getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
  auto t_diff_u =
      getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'l', 3> l;
  auto get_ftensor2 = [](auto &v) {
    return FTensor::Tensor1<FTensor::PackPtr<double *, 6>, 6>(
        &v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
  };
 
  int nb_base_functions = data.getN().size2();
  auto t_row_base_fun = data.getFTensor0N();
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    double a = v * t_w;
    auto t_nf = get_ftensor2(nF);
 
    FTensor::Tensor3<double, 3, 3, size_symm> t_Ldiff_u;
    t_Ldiff_u(i, j, L) = t_diff_u(i, j, k, l) * t_L(k, l, L);
 
    FTensor::Tensor2_symmetric<double, 3> t_viscous_P;
    t_viscous_P(i, j) =
        alphaU *
        t_dot_log_u(i,
                    j); // That is cheat, should be split on axial and deviator
 
    FTensor::Tensor1<double, size_symm> t_residual;
    t_residual(L) =
        a * (t_approx_P_adjoint_log_du(L) - t_Ldiff_u(i, j, L) * t_P(i, j) -
             t_L(i, j, L) * t_viscous_P(i, j));
 
    int bb = 0;
    for (; bb != nb_dofs / 6; ++bb) {
      t_nf(L) -= t_row_base_fun * t_residual(L);
      ++t_nf;
      ++t_row_base_fun;
    }
    for (; bb != nb_base_functions; ++bb)
      ++t_row_base_fun;
 
    ++t_w;
    ++t_approx_P_adjoint_log_du;
    ++t_P;
    ++t_dot_log_u;
    ++t_diff_u;
  }
  MoFEMFunctionReturn(0);
}
 
VolUserDataOperator *PhysicalEquations::returnOpSpatialPhysical(
    const std::string &field_name,
    boost::shared_ptr<DataAtIntegrationPts> data_ptr, const double alpha_u) {
  return new OpSpatialPhysical(field_name, data_ptr, alpha_u);
}
 
struct OpSpatialPhysical_du_du : public OpAssembleVolume {
  OpSpatialPhysical_du_du(std::string row_field, std::string col_field,
                          boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                          const double alpha)
      : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL, false),
        alphaU(alpha) {
    sYmm = false;
  }
  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
  MoFEMErrorCode integratePiola(EntData &row_data, EntData &col_data);
 
private:
  const double alphaU;
  MatrixDouble dP;
 
};
 
MoFEMErrorCode OpSpatialPhysical_du_du::integrate(EntData &row_data,
                                                  EntData &col_data) {
  return integratePiola(row_data, col_data);
}
 
MoFEMErrorCode OpSpatialPhysical_du_du::integratePiola(EntData &row_data,
                                                       EntData &col_data) {
  MoFEMFunctionBegin;
 
  FTensor::Index<'L', size_symm> L;
  FTensor::Index<'J', size_symm> J;
  auto t_L = symm_L_tensor();
  auto t_diff = diff_tensor();
 
  int nb_integration_pts = row_data.getN().size1();
  int row_nb_dofs = row_data.getIndices().size();
  int col_nb_dofs = col_data.getIndices().size();
 
  auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
    return FTensor::Tensor2<FTensor::PackPtr<double *, 6>, size_symm,
                            size_symm>(
 
        &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 0, c + 3),
        &m(r + 0, c + 4), &m(r + 0, c + 5),
 
        &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 1, c + 3),
        &m(r + 1, c + 4), &m(r + 1, c + 5),
 
        &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2), &m(r + 2, c + 3),
        &m(r + 2, c + 4), &m(r + 2, c + 5),
 
        &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2), &m(r + 3, c + 3),
        &m(r + 3, c + 4), &m(r + 3, c + 5),
 
        &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2), &m(r + 4, c + 3),
        &m(r + 4, c + 4), &m(r + 4, c + 5),
 
        &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2), &m(r + 5, c + 3),
        &m(r + 5, c + 4), &m(r + 5, c + 5)
 
    );
  };
  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;
 
  auto v = getVolume();
  auto t_w = getFTensor0IntegrationWeight();
 
  auto get_dP = [&]() {
    dP.resize(size_symm * size_symm, nb_integration_pts, false);
    auto ts_a = getTSa();
 
    auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->PAtPts);
    auto r_P_du = getFTensor4FromMat<3, 3, 3, 3>(dataAtPts->P_du);
    auto t_approx_P_adjoint_dstretch =
        getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
    auto t_dot_log_u =
        getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
    auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
    auto t_diff_u =
        getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
    auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
    auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
    auto &nbUniq = dataAtPts->nbUniq;
 
    auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
    for (auto gg = 0; gg != nb_integration_pts; ++gg) {
 
      FTensor::Tensor2<double, 3, 3> t_deltaP;
      t_deltaP(i, j) = t_approx_P_adjoint_dstretch(i, j) - t_P(i, j);
 
      FTensor::Tensor3<double, 3, 3, size_symm> t_Ldiff_u;
      t_Ldiff_u(i, j, L) = t_diff_u(i, j, m, n) * t_L(m, n, L);
      t_dP(L, J) =
          -t_Ldiff_u(i, j, L) * (r_P_du(i, j, k, l) * t_Ldiff_u(k, l, J));
      // viscous stress derivative
      t_dP(L, J) -= (alphaU * ts_a) *
                    (t_L(i, j, L) * (t_diff(i, j, k, l) * t_L(k, l, J)));
 
      if (EshelbianCore::stretchSelector > LINEAR &&
          EshelbianCore::gradApproximator > MODERATE_ROT) {
        FTensor::Tensor2_symmetric<double, 3> t_deltaP_sym;
        t_deltaP_sym(i, j) = (t_deltaP(i, j) || t_deltaP(j, i));
        t_deltaP_sym(i, j) /= 2.0;
        auto t_diff2_uP2 = EigenMatrix::getDiffDiffMat(
            t_eigen_vals, t_eigen_vecs, EshelbianCore::f, EshelbianCore::d_f,
            EshelbianCore::dd_f, t_deltaP_sym, nbUniq[gg]);
        t_dP(L, J) += t_L(i, j, L) * (t_diff2_uP2(i, j, k, l) * t_L(k, l, J));
      }
 
      ++t_P;
      ++t_dP;
      ++r_P_du;
      ++t_approx_P_adjoint_dstretch;
      ++t_dot_log_u;
      ++t_u;
      ++t_diff_u;
      ++t_eigen_vals;
      ++t_eigen_vecs;
    }
 
    return getFTensor2FromMat<size_symm, size_symm>(dP);
  };
 
  int row_nb_base_functions = row_data.getN().size2();
  auto t_row_base_fun = row_data.getFTensor0N();
 
  auto t_dP = get_dP();
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    double a = v * t_w;
 
    int rr = 0;
    for (; rr != row_nb_dofs / 6; ++rr) {
      auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
      auto t_m = get_ftensor2(K, 6 * rr, 0);
      for (int cc = 0; cc != col_nb_dofs / 6; ++cc) {
        const double b = a * t_row_base_fun * t_col_base_fun;
        t_m(L, J) -= b * t_dP(L, J);
        ++t_m;
        ++t_col_base_fun;
      }
      ++t_row_base_fun;
    }
 
    for (; rr != row_nb_base_functions; ++rr) {
      ++t_row_base_fun;
    }
 
    ++t_w;
    ++t_dP;
  }
  MoFEMFunctionReturn(0);
}
 
VolUserDataOperator *PhysicalEquations::returnOpSpatialPhysical_du_du(
    std::string row_field, std::string col_field,
    boost::shared_ptr<DataAtIntegrationPts> data_ptr, const double alpha) {
  return new OpSpatialPhysical_du_du(row_field, col_field, data_ptr, alpha);
}
 
VolUserDataOperator *PhysicalEquations::returnOpCalculateEnergy(
    boost::shared_ptr<DataAtIntegrationPts> data_ptr,
    boost::shared_ptr<double> total_energy_ptr) {
  return nullptr;
}
 
VolUserDataOperator *PhysicalEquations::returnOpCalculateStretchFromStress(
    boost::shared_ptr<DataAtIntegrationPts> data_ptr,
    boost::shared_ptr<PhysicalEquations> physics_ptr) {
  return nullptr;
}
 
VolUserDataOperator *PhysicalEquations::returnOpCalculateVarStretchFromStress(
    boost::shared_ptr<DataAtIntegrationPts> data_ptr,
    boost::shared_ptr<PhysicalEquations> physics_ptr) {
  return nullptr;
};
 
struct OpGetScale : public VolUserDataOperator {
  OpGetScale(boost::shared_ptr<double> scale_ptr)
      : VolUserDataOperator(NOSPACE, VolUserDataOperator::OPSPACE),
        scalePtr(scale_ptr) {}
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data) {
    MoFEMFunctionBegin;
    *(scalePtr) = 1.;
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<double> scalePtr;
};
 
VolUserDataOperator *PhysicalEquations::returnOpSetScale(
    boost::shared_ptr<double> scale_ptr,
    boost::shared_ptr<PhysicalEquations> physics_ptr) {
  return new OpGetScale(scale_ptr);
};
 
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::gradApproximator) {
    case NO_H1_CONFIGURATION:
      t_h1(i, j) = t_kd(i, j);
      physicsPtr->get_h()(i, j) = iu(i, j);
      break;
    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
      SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
               "ADOL-C function evaluation with error r = %d", r);
    }
 
    switch (EshelbianCore::gradApproximator) {
    case NO_H1_CONFIGURATION:
      r_P(i, j) = physicsPtr->get_P()(i, j);
      break;
    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::gradApproximator) {
    case NO_H1_CONFIGURATION:
      t_h1(i, j) = t_kd(i, j);
      physicsPtr->get_h()(i, j) = iu(i, k) * t_h1(k, j);
      break;
    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::gradApproximator) {
    case NO_H1_CONFIGURATION: {
      FTensor::Tensor4<double, 3, 3, 3, 3> t_h1_du;
      t_h1_du(i, j, m, n) = t_kd(i, m) * t_kd(j, n);
      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 LARGE_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 MODERATE_ROT:
    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);
}
 
} // namespace EshelbianPlasticity
 
#include <impl/HMHHStVenantKirchhoff.cpp>
#include <impl/HMHNeohookean.cpp>
#include <impl/HMHPMooneyRivlinWriggersEq63.cpp>
#include <impl/HMHHencky.cpp>
 
namespace EshelbianPlasticity {
 
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);
}
MoFEMErrorCode EshelbianCore::addMaterial_HMHNeohookean(const int tape,
                                                        const double c10,
                                                        const double K) {
  MoFEMFunctionBegin;
  physicalEquations = boost::make_shared<HMHNeohookean>(mField, c10, K);
  CHKERR physicalEquations->recordTape(tape, nullptr);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode EshelbianCore::addMaterial_Hencky(double E, double nu) {
  MoFEMFunctionBegin;
  physicalEquations = boost::make_shared<HMHHencky>(mField, E, nu);
  MoFEMFunctionReturn(0);
}
 
}; // namespace EshelbianPlasticity