656 {
658
659 auto neohookean_ptr =
660 boost::dynamic_pointer_cast<HMHNeohookean>(
dataAtPts->physicsPtr);
661 if (!neohookean_ptr) {
663 "Pointer to HMHNeohookean is null");
664 }
665 auto [def_c10, def_K] =
667
668 double c10 = def_c10 / neohookean_ptr->eqScaling;
669 double alpha_u =
alphaU / neohookean_ptr->eqScaling;
670 double lambda = def_K / neohookean_ptr->eqScaling;
671
672 double alpha_grad_u =
673 neohookean_ptr->alphaGradU / neohookean_ptr->eqScaling;
674
677
680
683
684 int nb_integration_pts = row_data.
getN().size1();
685 int row_nb_dofs = row_data.
getIndices().size();
686 int col_nb_dofs = col_data.
getIndices().size();
687
691
692 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 0,
c + 3),
693 &
m(r + 0,
c + 4), &
m(r + 0,
c + 5),
694
695 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 1,
c + 3),
696 &
m(r + 1,
c + 4), &
m(r + 1,
c + 5),
697
698 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2), &
m(r + 2,
c + 3),
699 &
m(r + 2,
c + 4), &
m(r + 2,
c + 5),
700
701 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2), &
m(r + 3,
c + 3),
702 &
m(r + 3,
c + 4), &
m(r + 3,
c + 5),
703
704 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2), &
m(r + 4,
c + 3),
705 &
m(r + 4,
c + 4), &
m(r + 4,
c + 5),
706
707 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2), &
m(r + 5,
c + 3),
708 &
m(r + 5,
c + 4), &
m(r + 5,
c + 5)
709
710 );
711 };
712
719
723
724 int row_nb_base_functions = row_data.
getN().size2();
727
728 auto t_grad_h1 = getFTensor2FromMat<3, 3>(
dataAtPts->wGradH1AtPts);
729 auto t_diff_u =
730 getFTensor4DdgFromMat<3, 3, 1>(
dataAtPts->diffStretchTensorAtPts);
731 auto t_log_u =
732 getFTensor2SymmetricFromMat<3>(
dataAtPts->logStretchTensorAtPts);
733 auto t_log_u2_h1 =
734 getFTensor2SymmetricFromMat<3>(
dataAtPts->logStretch2H1AtPts);
735 auto t_u = getFTensor2SymmetricFromMat<3>(
dataAtPts->stretchTensorAtPts);
736 auto t_approx_P_adjoint__dstretch =
737 getFTensor2FromMat<3, 3>(
dataAtPts->adjointPdstretchAtPts);
738 auto t_eigen_vals = getFTensor1FromMat<3>(
dataAtPts->eigenVals);
739 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(
dataAtPts->eigenVecs);
741 auto t_nb_uniq =
743
744 auto no_h1 = [&]() {
746
747 for (int gg = 0; gg != nb_integration_pts; ++gg) {
749 ++t_w;
750
752 auto d_neohookean = [c10,
lambda](
auto v) {
754 };
755
757 t_eigen_vals, t_eigen_vecs, neohookean, d_neohookean, t_nb_uniq);
758
759 const auto tr = t_log_u(
i,
i);
761 t_dP(L,
J) = -t_L(
i,
j, L) * ((t_diff_neohookean(
i,
j,
k,
l) +
763 t_kd_sym(
i,
j) * t_kd_sym(
k,
l)) *
765 t_dP(L,
J) -= (alpha_u * ts_a) *
766 (t_L(
i,
j, L) * (t_diff(
i,
j,
k,
l) * t_L(
k,
l,
J)));
767
768 if constexpr (1) {
770 t_deltaP(
i,
j) = (t_approx_P_adjoint__dstretch(
i,
j) ||
771 t_approx_P_adjoint__dstretch(
j,
i)) /
772 2.;
776 t_dP(L,
J) += t_L(
i,
j, L) * (t_diff2_uP(
i,
j,
k,
l) * t_L(
k,
l,
J));
777 }
778 ++t_approx_P_adjoint__dstretch;
779 ++t_log_u;
780 ++t_eigen_vals;
781 ++t_eigen_vecs;
782 ++t_nb_uniq;
783
784 int rr = 0;
785 for (; rr != row_nb_dofs /
size_symm; ++rr) {
788
789 auto t_m = get_ftensor2(
K, 6 * rr, 0);
790 for (
int cc = 0; cc != col_nb_dofs /
size_symm; ++cc) {
791 double b =
a * t_row_base_fun * t_col_base_fun;
792 double c = (
a * alpha_grad_u * ts_a) *
793 (t_row_grad_fun(
i) * t_col_grad_fun(
i));
794 t_m(L,
J) -= b * t_dP(L,
J);
795 t_m(L,
J) +=
c * t_kd_sym(L,
J);
796
797 ++t_m;
798 ++t_col_base_fun;
799 ++t_col_grad_fun;
800 }
801 ++t_row_base_fun;
802 ++t_row_grad_fun;
803 }
804
805 for (; rr != row_nb_base_functions; ++rr) {
806 ++t_row_base_fun;
807 ++t_row_grad_fun;
808 }
809
810 }
812 };
813
814 auto large = [&]() {
817 "Not implemented for Neo-Hookean (used ADOL-C)");
819 };
820
824 break;
828 break;
829 default:
831 "gradApproximator not handled");
832 };
833
835}
#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 getVolume() const
element volume (linear geometry)
MatrixDouble K
local tangent matrix
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts