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
678 auto get_ftensor2 = [](MatrixDouble &
m,
const int r,
const int c) {
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
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.
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
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 dofs 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)