70                      "extract block data failed");
 
   72#ifdef NEOHOOKEAN_SCALING 
   87                          "Stretch selector is not equal to LOG");
 
   91                            "Exponent base is not equal to exp(1)");
 
 
   99      if (b.blockEnts.find(ent) != b.blockEnts.end()) {
 
  100        return std::make_pair(b.c10, b.K);
 
  105                        "Block not found for entity handle. If you mat set " 
  106                        "block, set it to all elements");
 
 
  112               boost::shared_ptr<PhysicalEquations> physics_ptr)
 
 
  116    MoFEMErrorCode 
doWork(
int side, EntityType type,
 
  117                          EntitiesFieldData::EntData &data) {
 
  119      auto neohookean_ptr =
 
  120          boost::dynamic_pointer_cast<HMHNeohookean>(
physicsPtr);
 
  121      *(
scalePtr) = neohookean_ptr->eqScaling;
 
 
 
  132                   boost::shared_ptr<PhysicalEquations> physics_ptr) {
 
  134    return new OpGetScale(scale_ptr, physics_ptr);
 
 
  138    using EshelbianPlasticity::OpJacobian::OpJacobian;
 
 
  145                   boost::shared_ptr<DataAtIntegrationPts> data_ptr,
 
  146                   boost::shared_ptr<PhysicalEquations> physics_ptr) {
 
  149                                                 data_ptr, physics_ptr);
 
  151      return (
new OpJacobian(tag, eval_rhs, eval_lhs, data_ptr, physics_ptr));
 
 
  157    PetscOptionsBegin(PETSC_COMM_WORLD, 
"neo_hookean_", 
"", 
"none");
 
  167    CHKERR PetscOptionsScalar(
"-viscosity_alpha_grad_u", 
"viscosity", 
"",
 
 
  191            (boost::format(
"%s(.*)") % 
"MAT_NEOHOOKEAN").str()
 
 
  203    for (
auto m : meshset_vec_ptr) {
 
  205      std::vector<double> block_data;
 
  206      CHKERR m->getAttributes(block_data);
 
  207      if (block_data.size() < 2) {
 
  209                "Expected that block has atleast two attributes");
 
  211      auto get_block_ents = [&]() {
 
  217      double c10 = block_data[0];
 
  218      double K = block_data[1];
 
  220      blockData.push_back({c10, K, get_block_ents()});
 
  222      MOFEM_LOG(
"EP", sev) << 
"MatBlock Neo-Hookean c10 = " 
  224                           << 
" K = " << 
blockData.back().K << 
" nb ents. = " 
 
  245        ih(
i, 
j) = (*t_h_ptr)(
i, 
j);
 
  252      enableMinMaxUsingAbs();
 
  283              "gradApproximator not handled");
 
 
  292                      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
 
  293                      const double alpha_u);
 
 
  303                          boost::shared_ptr<DataAtIntegrationPts> data_ptr,
 
  304                          const double alpha_u) {
 
 
  317        boost::shared_ptr<DataAtIntegrationPts> data_ptr,
 
  318        boost::shared_ptr<ExternalStrainVec> external_strain_vec_ptr,
 
  319        std::map<std::string, boost::shared_ptr<ScalingMethod>> smv);
 
 
  330      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
 
  331      boost::shared_ptr<ExternalStrainVec> external_strain_vec_ptr,
 
  332      std::map<std::string, boost::shared_ptr<ScalingMethod>> smv) {
 
  334                                               external_strain_vec_ptr, smv);
 
 
  339                            boost::shared_ptr<DataAtIntegrationPts> data_ptr,
 
 
  348      std::string row_field, std::string col_field,
 
  349      boost::shared_ptr<DataAtIntegrationPts> data_ptr, 
const double alpha) {
 
  352          row_field, col_field, data_ptr, alpha);
 
 
 
  386    boost::shared_ptr<DataAtIntegrationPts> data_ptr, 
const double alpha_u)
 
 
  391                        double circumcenter[3], 
double *xi, 
double *
eta,
 
  398  auto neohookean_ptr =
 
  399      boost::dynamic_pointer_cast<HMHNeohookean>(dataAtPts->physicsPtr);
 
  400  if (!neohookean_ptr) {
 
  402            "Pointer to HMHNeohookean is null");
 
  404  auto [def_c10, def_K] =
 
  405      neohookean_ptr->getMaterialParameters(getFEEntityHandle());
 
  407  double c10 = def_c10 / neohookean_ptr->eqScaling;
 
  408  double alpha_u = alphaU / neohookean_ptr->eqScaling;
 
  409  double K = def_K / neohookean_ptr->eqScaling;
 
  411  auto time_step = getTStimeStep();
 
  412  double alpha_grad_u =
 
  413      neohookean_ptr->alphaGradU / neohookean_ptr->eqScaling;
 
  422  int nb_integration_pts = data.
getN().size1();
 
  423  auto v = getVolume();
 
  424  auto t_w = getFTensor0IntegrationWeight();
 
  425  auto t_approx_P_adjoint_log_du =
 
  426      getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
 
  427  auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
 
  428  auto t_grad_h1 = getFTensor2FromMat<3, 3>(dataAtPts->wGradH1AtPts);
 
  430      getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
 
  432      getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
 
  434      getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
 
  436      getFTensor2FromMat<6, 3>(dataAtPts->gradLogStretchDotTensorAtPts);
 
  438      getFTensor2SymmetricFromMat<3>(dataAtPts->logStretch2H1AtPts);
 
  440  auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
 
  441  auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
 
  442  auto &nbUniq = dataAtPts->nbUniq;
 
  455  auto get_ftensor2 = [](
auto &
v) {
 
  457        &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
 
  460  int nb_base_functions = data.
getN().size2();
 
  467    for (
int gg = 0; gg != nb_integration_pts; ++gg) {
 
  472      auto t_neohookean_hencky =
 
  478      const double tr = t_log_u(
i, 
i);
 
  480      t_P(L) = t_L(
i, 
j, L) * (t_neohookean_hencky(
i, 
j) +
 
  483      t_viscous_P(L) = alpha_u * (t_L(
i, 
j, L) * t_dot_log_u(
i, 
j));
 
  486      t_residual(L) = t_approx_P_adjoint_log_du(L) - t_P(L) - t_viscous_P(L);
 
  490      t_grad_residual(L, 
i) = alpha_grad_u * t_grad_log_u(L, 
i);
 
  491      t_grad_residual(L, 
i) *= 
a;
 
  493      ++t_approx_P_adjoint_log_du;
 
  498      auto t_nf = getFTensor1FromPtr<size_symm>(&*nF.data().begin());
 
  500      for (; bb != nb_dofs / 
size_symm; ++bb) {
 
  501        t_nf(L) -= t_row_base_fun * t_residual(L);
 
  502        t_nf(L) += t_grad_base_fun(
i) * t_grad_residual(L, 
i);
 
  507      for (; bb != nb_base_functions; ++bb) {
 
  519            "Not implemented for Neo-Hookean (used ADOL-C)");
 
  533            "gradApproximator not handled");
 
 
  541    boost::shared_ptr<DataAtIntegrationPts> data_ptr,
 
  542    boost::shared_ptr<ExternalStrainVec> external_strain_vec_ptr,
 
  543    std::map<std::string, boost::shared_ptr<ScalingMethod>> smv)
 
  545      externalStrainVecPtr(external_strain_vec_ptr), scalingMethodsMap{smv} {}
 
 
  551  auto neohookean_ptr =
 
  552      boost::dynamic_pointer_cast<HMHNeohookean>(dataAtPts->physicsPtr);
 
  553  if (!neohookean_ptr) {
 
  555            "Pointer to HMHNeohookean is null");
 
  566  for (
auto &ext_strain_block : (*externalStrainVecPtr)) {
 
  568    if (ext_strain_block.ents.find(fe_ent) != ext_strain_block.ents.end()) {
 
  571      if (scalingMethodsMap.find(ext_strain_block.blockName) !=
 
  572          scalingMethodsMap.end()) {
 
  574            scalingMethodsMap.at(ext_strain_block.blockName)->getScale(time);
 
  577            << 
"No scaling method found for " << ext_strain_block.blockName;
 
  580      auto [def_c10, def_K] =
 
  581          neohookean_ptr->getMaterialParameters(getFEEntityHandle());
 
  582      double c10 = def_c10 / neohookean_ptr->eqScaling;
 
  584      double external_strain_val = 
scale * ext_strain_block.val;
 
  585      double K = ext_strain_block.bulkModulusK;
 
  592      int nb_integration_pts = data.
getN().size1();
 
  593      auto v = getVolume();
 
  594      auto t_w = getFTensor0IntegrationWeight();
 
  603      int nb_base_functions = data.
getN().size2();
 
  606      const double tr = 3.0 * external_strain_val;
 
  607      const double sigma_J = K * tr;
 
  609      for (
int gg = 0; gg != nb_integration_pts; ++gg) {
 
  615        t_residual(L) += (t_L(
i, 
j, L) * 
t_kd(
i, 
j)) * sigma_J;
 
  618        auto t_nf = getFTensor1FromPtr<size_symm>(&*nF.data().begin());
 
  620        for (; bb != nb_dofs / 
size_symm; ++bb) {
 
  621          t_nf(L) -= t_row_base_fun * t_residual(L);
 
  625        for (; bb != nb_base_functions; ++bb) {
 
 
  635    std::string row_field, std::string col_field,
 
  636    boost::shared_ptr<DataAtIntegrationPts> data_ptr, 
const double alpha)
 
 
  648  auto neohookean_ptr =
 
  649      boost::dynamic_pointer_cast<HMHNeohookean>(dataAtPts->physicsPtr);
 
  650  if (!neohookean_ptr) {
 
  652            "Pointer to HMHNeohookean is null");
 
  654  auto [def_c10, def_K] =
 
  655      neohookean_ptr->getMaterialParameters(getFEEntityHandle());
 
  657  double c10 = def_c10 / neohookean_ptr->eqScaling;
 
  658  double alpha_u = alphaU / neohookean_ptr->eqScaling;
 
  659  double lambda = def_K / neohookean_ptr->eqScaling;
 
  661  auto time_step = getTStimeStep();
 
  662  double alpha_grad_u =
 
  663      neohookean_ptr->alphaGradU / neohookean_ptr->eqScaling;
 
  674  int nb_integration_pts = row_data.
getN().size1();
 
  675  int row_nb_dofs = row_data.
getIndices().size();
 
  676  int col_nb_dofs = col_data.
getIndices().size();
 
  678  auto get_ftensor2 = [](MatrixDouble &
m, 
const int r, 
const int c) {
 
  682        &
m(r + 0, 
c + 0), &
m(r + 0, 
c + 1), &
m(r + 0, 
c + 2), &
m(r + 0, 
c + 3),
 
  683        &
m(r + 0, 
c + 4), &
m(r + 0, 
c + 5),
 
  685        &
m(r + 1, 
c + 0), &
m(r + 1, 
c + 1), &
m(r + 1, 
c + 2), &
m(r + 1, 
c + 3),
 
  686        &
m(r + 1, 
c + 4), &
m(r + 1, 
c + 5),
 
  688        &
m(r + 2, 
c + 0), &
m(r + 2, 
c + 1), &
m(r + 2, 
c + 2), &
m(r + 2, 
c + 3),
 
  689        &
m(r + 2, 
c + 4), &
m(r + 2, 
c + 5),
 
  691        &
m(r + 3, 
c + 0), &
m(r + 3, 
c + 1), &
m(r + 3, 
c + 2), &
m(r + 3, 
c + 3),
 
  692        &
m(r + 3, 
c + 4), &
m(r + 3, 
c + 5),
 
  694        &
m(r + 4, 
c + 0), &
m(r + 4, 
c + 1), &
m(r + 4, 
c + 2), &
m(r + 4, 
c + 3),
 
  695        &
m(r + 4, 
c + 4), &
m(r + 4, 
c + 5),
 
  697        &
m(r + 5, 
c + 0), &
m(r + 5, 
c + 1), &
m(r + 5, 
c + 2), &
m(r + 5, 
c + 3),
 
  698        &
m(r + 5, 
c + 4), &
m(r + 5, 
c + 5)
 
  710  auto v = getVolume();
 
  711  auto ts_a = getTSa();
 
  712  auto t_w = getFTensor0IntegrationWeight();
 
  714  int row_nb_base_functions = row_data.
getN().size2();
 
  718  auto t_grad_h1 = getFTensor2FromMat<3, 3>(dataAtPts->wGradH1AtPts);
 
  720      getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
 
  722      getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
 
  724      getFTensor2SymmetricFromMat<3>(dataAtPts->logStretch2H1AtPts);
 
  725  auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
 
  726  auto t_approx_P_adjoint__dstretch =
 
  727      getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
 
  728  auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
 
  729  auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
 
  730  auto &nbUniq = dataAtPts->nbUniq;
 
  737    for (
int gg = 0; gg != nb_integration_pts; ++gg) {
 
  742      auto d_neohookean = [c10, 
lambda](
auto v) {
 
  747          t_eigen_vals, t_eigen_vecs, neohookean, d_neohookean, t_nb_uniq);
 
  749      const auto tr = t_log_u(
i, 
i);
 
  751      t_dP(L, 
J) = -t_L(
i, 
j, L) * ((t_diff_neohookean(
i, 
j, 
k, 
l) +
 
  753                                         t_kd_sym(
i, 
j) * t_kd_sym(
k, 
l)) *
 
  755      t_dP(L, 
J) -= (alpha_u * ts_a) *
 
  756                    (t_L(
i, 
j, L) * (t_diff(
i, 
j, 
k, 
l) * t_L(
k, 
l, 
J)));
 
  760        t_deltaP(
i, 
j) = (t_approx_P_adjoint__dstretch(
i, 
j) ||
 
  761                          t_approx_P_adjoint__dstretch(
j, 
i)) /
 
  766        t_dP(L, 
J) += t_L(
i, 
j, L) * (t_diff2_uP(
i, 
j, 
k, 
l) * t_L(
k, 
l, 
J));
 
  768      ++t_approx_P_adjoint__dstretch;
 
  775      for (; rr != row_nb_dofs / 
size_symm; ++rr) {
 
  779        auto t_m = get_ftensor2(K, 6 * rr, 0);
 
  780        for (
int cc = 0; cc != col_nb_dofs / 
size_symm; ++cc) {
 
  781          double b = 
a * t_row_base_fun * t_col_base_fun;
 
  782          double c = (
a * alpha_grad_u * ts_a) *
 
  783                     (t_row_grad_fun(
i) * t_col_grad_fun(
i));
 
  784          t_m(L, 
J) -= b * t_dP(L, 
J);
 
  785          t_m(L, 
J) += 
c * t_kd_sym(L, 
J);
 
  795      for (; rr != row_nb_base_functions; ++rr) {
 
  807            "Not implemented for Neo-Hookean (used ADOL-C)");
 
  821              "gradApproximator not handled");
 
 
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
#define FTENSOR_INDEX(DIM, I)
Kronecker Delta class symmetric.
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
const double v
phase velocity of light in medium (cm/ns)
const double n
refractive index of diffusive medium
FTensor::Index< 'J', DIM1 > J
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
auto getMat(A &&t_val, B &&t_vec, Fun< double > f)
Get the Mat object.
auto getDiffMat(A &&t_val, B &&t_vec, Fun< double > f, Fun< double > d_f, const int nb)
Get the Diff Mat object.
auto getDiffDiffMat(A &&t_val, B &&t_vec, Fun< double > f, Fun< double > d_f, Fun< double > dd_f, C &&t_S, const int nb)
Get the Diff Diff Mat object.
void tetcircumcenter_tp(double a[3], double b[3], double c[3], double d[3], double circumcenter[3], double *xi, double *eta, double *zeta)
ForcesAndSourcesCore::UserDataOperator UserDataOperator
static constexpr auto size_symm
constexpr IntegrationType I
double zeta
Viscous hardening.
constexpr auto field_name
FTensor::Index< 'm', 3 > m
static enum StretchSelector stretchSelector
static double dynamicTime
static enum RotSelector gradApproximator
static PetscBool dynamicRelaxation
static double exponentBase
static boost::function< double(const double)> f
static boost::function< double(const double)> dd_f
static boost::function< double(const double)> d_f
OpGetScale(boost::shared_ptr< double > scale_ptr, boost::shared_ptr< PhysicalEquations > physics_ptr)
boost::shared_ptr< double > scalePtr
boost::shared_ptr< PhysicalEquations > physicsPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
MoFEMErrorCode evaluateRhs(EntData &data)
MoFEMErrorCode evaluateLhs(EntData &data)
MoFEMErrorCode integrate(EntData &data)
boost::shared_ptr< ExternalStrainVec > externalStrainVecPtr
std::map< std::string, boost::shared_ptr< ScalingMethod > > scalingMethodsMap
OpSpatialPhysicalExternalStrain(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< ExternalStrainVec > external_strain_vec_ptr, std::map< std::string, boost::shared_ptr< ScalingMethod > > smv)
OpSpatialPhysical_du_du(std::string row_field, std::string col_field, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const double alpha)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
OpSpatialPhysical(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const double alpha_u)
MoFEMErrorCode integrate(EntData &data)
static double fun_diff_neohookean_bulk(double K, double tr)
Definition of derivative of axiator of Neo-hookean function.
static double fun_d_neohookean(double c10, double v)
Definition of derivative of Neo-hookean function.
static double fun_neohookean_bulk(double K, double tr)
Definition of axiator of Neo-hookean function.
virtual VolUserDataOperator * returnOpSpatialPhysical(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const double alpha_u)
std::vector< BlockData > blockData
MoFEMErrorCode extractBlockData(std::vector< const CubitMeshSets * > meshset_vec_ptr, Sev sev)
MoFEMErrorCode extractBlockData(Sev sev)
static constexpr int numberOfDependentVariables
VolUserDataOperator * returnOpSetScale(boost::shared_ptr< double > scale_ptr, boost::shared_ptr< PhysicalEquations > physics_ptr)
static double fun_neohookean(double c10, double v)
Definition of Neo-hookean function.
MoFEMErrorCode getOptions()
MoFEM::Interface & mField
auto getMaterialParameters(EntityHandle ent)
static constexpr int numberOfActiveVariables
virtual VolUserDataOperator * returnOpSpatialPhysicalExternalStrain(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< ExternalStrainVec > external_strain_vec_ptr, std::map< std::string, boost::shared_ptr< ScalingMethod > > smv)
HMHNeohookean(MoFEM::Interface &m_field, const double c10, const double K)
VolUserDataOperator * returnOpSpatialPhysical_du_du(std::string row_field, std::string col_field, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const double alpha)
MoFEMErrorCode recordTape(const int tape, DTensor2Ptr *t_h_ptr)
UserDataOperator * returnOpJacobian(const int tag, const bool eval_rhs, const bool eval_lhs, boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< PhysicalEquations > physics_ptr)
virtual VolUserDataOperator * returnOpSpatialPhysical(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const double alpha_u)
virtual VolUserDataOperator * returnOpSpatialPhysical_du_du(std::string row_field, std::string col_field, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const double alpha)
virtual UserDataOperator * returnOpJacobian(const int tag, const bool eval_rhs, const bool eval_lhs, boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< PhysicalEquations > physics_ptr)
virtual moab::Interface & get_moab()=0
bool sYmm
If true assume that matrix is symmetric structure.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorInt & getIndices() const
Get global indices of degrees of freedom on entity.
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
PetscReal ts_t
Current time value.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.