645                                                                     {
  647 
  648  auto neohookean_ptr =
  649      boost::dynamic_pointer_cast<HMHNeohookean>(
dataAtPts->physicsPtr);
 
  650  if (!neohookean_ptr) {
  652            "Pointer to HMHNeohookean is null");
  653  }
  654  auto [def_c10, def_K] =
  656 
  657  double c10 = def_c10 / neohookean_ptr->eqScaling;
  658  double alpha_u = 
alphaU / neohookean_ptr->eqScaling;
 
  659  double lambda = def_K / neohookean_ptr->eqScaling;
 
  660 
  662  double alpha_grad_u =
  663      neohookean_ptr->alphaGradU / neohookean_ptr->eqScaling;
  664 
  667 
  670 
  673 
  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();
 
  677 
  681 
  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),
 
  684 
  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),
 
  687 
  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),
 
  690 
  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),
 
  693 
  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),
 
  696 
  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)
 
  699 
  700    );
  701  };
  702 
  709 
  713 
  714  int row_nb_base_functions = row_data.
getN().size2();
 
  717 
  718  auto t_grad_h1 = getFTensor2FromMat<3, 3>(
dataAtPts->wGradH1AtPts);
 
  719  auto t_diff_u =
  720      getFTensor4DdgFromMat<3, 3, 1>(
dataAtPts->diffStretchTensorAtPts);
 
  721  auto t_log_u =
  722      getFTensor2SymmetricFromMat<3>(
dataAtPts->logStretchTensorAtPts);
 
  723  auto t_log_u2_h1 =
  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);
 
  731  auto t_nb_uniq =
  733 
  734  auto no_h1 = [&]() {
  736 
  737    for (int gg = 0; gg != nb_integration_pts; ++gg) {
  739      ++t_w;
  740 
  742      auto d_neohookean = [c10, 
lambda](
auto v) {
 
  744      };
  745 
  747          t_eigen_vals, t_eigen_vecs, neohookean, d_neohookean, t_nb_uniq);
  748 
  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)));
 
  757 
  758      if constexpr (1) {
  760        t_deltaP(
i, 
j) = (t_approx_P_adjoint__dstretch(
i, 
j) ||
 
  761                          t_approx_P_adjoint__dstretch(
j, 
i)) /
 
  762                         2.;
  766        t_dP(L, 
J) += t_L(
i, 
j, L) * (t_diff2_uP(
i, 
j, 
k, 
l) * t_L(
k, 
l, 
J));
 
  767      }
  768      ++t_approx_P_adjoint__dstretch;
  769      ++t_log_u;
  770      ++t_eigen_vals;
  771      ++t_eigen_vecs;
  772      ++t_nb_uniq;
  773 
  774      int rr = 0;
  775      for (; rr != row_nb_dofs / 
size_symm; ++rr) {
 
  778 
  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);
 
  786 
  787          ++t_m;
  788          ++t_col_base_fun;
  789          ++t_col_grad_fun;
  790        }
  791        ++t_row_base_fun;
  792        ++t_row_grad_fun;
  793      }
  794 
  795      for (; rr != row_nb_base_functions; ++rr) {
  796        ++t_row_base_fun;
  797        ++t_row_grad_fun;
  798      }
  799 
  800    }
  802  };
  803 
  804  auto large = [&]() {
  807            "Not implemented for Neo-Hookean (used ADOL-C)");
  809  };
  810 
  814      break;
  818      break;
  819    default:
  821              "gradApproximator not handled");
  822    };
  823 
  825}
#define FTENSOR_INDEX(DIM, I)
Kronecker Delta class symmetric.
#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.
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 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.
static constexpr auto size_symm
UBlasMatrix< double > MatrixDouble
FTensor::Index< 'm', 3 > m
static enum RotSelector gradApproximator
static boost::function< double(const double)> f
static boost::function< double(const double)> dd_f
static boost::function< double(const double)> d_f
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(double c10, double v)
Definition of Neo-hookean function.
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.
auto getFTensor0IntegrationWeight()
Get integration weights.
double getTStimeStep() const
double getVolume() const
element volume (linear geometry)
MatrixDouble K
local tangent matrix
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts