651 {
653
654 auto neohookean_ptr =
655 boost::dynamic_pointer_cast<HMHNeohookean>(
dataAtPts->physicsPtr);
656 if (!neohookean_ptr) {
658 "Pointer to HMHNeohookean is null");
659 }
660 auto [def_c10, def_K] =
662
663 double c10 = def_c10 / neohookean_ptr->eqScaling;
664 double alpha_u =
alphaU / neohookean_ptr->eqScaling;
665 double lambda = def_K / neohookean_ptr->eqScaling;
666
668 double alpha_grad_u =
669 neohookean_ptr->alphaGradU / neohookean_ptr->eqScaling;
670
673
676
679
680 int nb_integration_pts = row_data.
getN().size1();
681 int row_nb_dofs = row_data.
getIndices().size();
682 int col_nb_dofs = col_data.
getIndices().size();
683
687
688 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 0,
c + 3),
689 &
m(r + 0,
c + 4), &
m(r + 0,
c + 5),
690
691 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 1,
c + 3),
692 &
m(r + 1,
c + 4), &
m(r + 1,
c + 5),
693
694 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2), &
m(r + 2,
c + 3),
695 &
m(r + 2,
c + 4), &
m(r + 2,
c + 5),
696
697 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2), &
m(r + 3,
c + 3),
698 &
m(r + 3,
c + 4), &
m(r + 3,
c + 5),
699
700 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2), &
m(r + 4,
c + 3),
701 &
m(r + 4,
c + 4), &
m(r + 4,
c + 5),
702
703 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2), &
m(r + 5,
c + 3),
704 &
m(r + 5,
c + 4), &
m(r + 5,
c + 5)
705
706 );
707 };
708
715
719
720 int row_nb_base_functions = row_data.
getN().size2();
723
724 auto t_grad_h1 = getFTensor2FromMat<3, 3>(
dataAtPts->wGradH1AtPts);
725 auto t_diff_u =
726 getFTensor4DdgFromMat<3, 3, 1>(
dataAtPts->diffStretchTensorAtPts);
727 auto t_log_u =
728 getFTensor2SymmetricFromMat<3>(
dataAtPts->logStretchTensorAtPts);
729 auto t_log_u2_h1 =
730 getFTensor2SymmetricFromMat<3>(
dataAtPts->logStretch2H1AtPts);
731 auto t_u = getFTensor2SymmetricFromMat<3>(
dataAtPts->stretchTensorAtPts);
732 auto t_approx_P_adjoint__dstretch =
733 getFTensor2FromMat<3, 3>(
dataAtPts->adjointPdstretchAtPts);
734 auto t_eigen_vals = getFTensor1FromMat<3>(
dataAtPts->eigenVals);
735 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(
dataAtPts->eigenVecs);
737 auto t_nb_uniq =
739
740 auto no_h1 = [&]() {
742
743 for (int gg = 0; gg != nb_integration_pts; ++gg) {
745 ++t_w;
746
748 auto d_neohookean = [c10,
lambda](
auto v) {
750 };
751
753 t_eigen_vals, t_eigen_vecs, neohookean, d_neohookean, t_nb_uniq);
754
755 const auto tr = t_log_u(
i,
i);
757 t_dP(L,
J) = -t_L(
i,
j, L) * ((t_diff_neohookean(
i,
j,
k,
l) +
759 t_kd_sym(
i,
j) * t_kd_sym(
k,
l)) *
761 t_dP(L,
J) -= (alpha_u * ts_a) *
762 (t_L(
i,
j, L) * (t_diff(
i,
j,
k,
l) * t_L(
k,
l,
J)));
763
764 if constexpr (1) {
766 t_deltaP(
i,
j) = (t_approx_P_adjoint__dstretch(
i,
j) ||
767 t_approx_P_adjoint__dstretch(
j,
i)) /
768 2.;
772 t_dP(L,
J) += t_L(
i,
j, L) * (t_diff2_uP(
i,
j,
k,
l) * t_L(
k,
l,
J));
773 }
774 ++t_approx_P_adjoint__dstretch;
775 ++t_log_u;
776 ++t_eigen_vals;
777 ++t_eigen_vecs;
778 ++t_nb_uniq;
779
780 int rr = 0;
781 for (; rr != row_nb_dofs /
size_symm; ++rr) {
784
785 auto t_m = get_ftensor2(
K, 6 * rr, 0);
786 for (
int cc = 0; cc != col_nb_dofs /
size_symm; ++cc) {
787 double b =
a * t_row_base_fun * t_col_base_fun;
788 double c = (
a * alpha_grad_u * ts_a) *
789 (t_row_grad_fun(
i) * t_col_grad_fun(
i));
790 t_m(L,
J) -= b * t_dP(L,
J);
791 t_m(L,
J) +=
c * t_kd_sym(L,
J);
792
793 ++t_m;
794 ++t_col_base_fun;
795 ++t_col_grad_fun;
796 }
797 ++t_row_base_fun;
798 ++t_row_grad_fun;
799 }
800
801 for (; rr != row_nb_base_functions; ++rr) {
802 ++t_row_base_fun;
803 ++t_row_grad_fun;
804 }
805
806 }
808 };
809
810 auto large = [&]() {
813 "Not implemented for Neo-Hookean (used ADOL-C)");
815 };
816
820 break;
824 break;
825 default:
827 "gradApproximator not handled");
828 };
829
831}
#define FTENSOR_INDEX(DIM, I)
constexpr int SPACE_DIM
[Define dimension]
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