666 {
668
669 auto neohookean_ptr =
670 boost::dynamic_pointer_cast<HMHNeohookean>(
dataAtPts->physicsPtr);
671 if (!neohookean_ptr) {
673 "Pointer to HMHNeohookean is null");
674 }
675 auto [def_c10, def_K] =
677
678 double c10 = def_c10 / neohookean_ptr->eqScaling;
679 double alpha_u =
alphaU / neohookean_ptr->eqScaling;
680 double lambda = def_K / neohookean_ptr->eqScaling;
681
682 double alpha_grad_u =
683 neohookean_ptr->alphaGradU / neohookean_ptr->eqScaling;
684
687
690
693
694 int nb_integration_pts = row_data.
getN().size1();
695 int row_nb_dofs = row_data.
getIndices().size();
696 int col_nb_dofs = col_data.
getIndices().size();
697
701
702 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 0,
c + 3),
703 &
m(r + 0,
c + 4), &
m(r + 0,
c + 5),
704
705 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 1,
c + 3),
706 &
m(r + 1,
c + 4), &
m(r + 1,
c + 5),
707
708 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2), &
m(r + 2,
c + 3),
709 &
m(r + 2,
c + 4), &
m(r + 2,
c + 5),
710
711 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2), &
m(r + 3,
c + 3),
712 &
m(r + 3,
c + 4), &
m(r + 3,
c + 5),
713
714 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2), &
m(r + 4,
c + 3),
715 &
m(r + 4,
c + 4), &
m(r + 4,
c + 5),
716
717 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2), &
m(r + 5,
c + 3),
718 &
m(r + 5,
c + 4), &
m(r + 5,
c + 5)
719
720 );
721 };
722
729
733
734 int row_nb_base_functions = row_data.
getN().size2();
737
738 auto t_grad_h1 = getFTensor2FromMat<3, 3>(
dataAtPts->wGradH1AtPts);
739 auto t_diff_u =
740 getFTensor4DdgFromMat<3, 3, 1>(
dataAtPts->diffStretchTensorAtPts);
741 auto t_log_u =
742 getFTensor2SymmetricFromMat<3>(
dataAtPts->logStretchTensorAtPts);
743 auto t_log_u2_h1 =
744 getFTensor2SymmetricFromMat<3>(
dataAtPts->logStretch2H1AtPts);
745 auto t_u = getFTensor2SymmetricFromMat<3>(
dataAtPts->stretchTensorAtPts);
746 auto t_approx_P_adjoint__dstretch =
747 getFTensor2FromMat<3, 3>(
dataAtPts->adjointPdstretchAtPts);
748 auto t_eigen_vals = getFTensor1FromMat<3>(
dataAtPts->eigenVals);
749 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(
dataAtPts->eigenVecs);
751 auto t_nb_uniq =
753
754 auto no_h1 = [&]() {
756
757 for (int gg = 0; gg != nb_integration_pts; ++gg) {
759 ++t_w;
760
762 auto d_neohookean = [c10,
lambda](
auto v) {
764 };
765
767 t_eigen_vals, t_eigen_vecs, neohookean, d_neohookean, t_nb_uniq);
768
769 const auto tr = t_log_u(
i,
i);
771 t_dP(L,
J) = -t_L(
i,
j, L) * ((t_diff_neohookean(
i,
j,
k,
l) +
773 t_kd_sym(
i,
j) * t_kd_sym(
k,
l)) *
775 t_dP(L,
J) -= (alpha_u * ts_a) *
776 (t_L(
i,
j, L) * (t_diff(
i,
j,
k,
l) * t_L(
k,
l,
J)));
777
778 if constexpr (1) {
780 t_deltaP(
i,
j) = (t_approx_P_adjoint__dstretch(
i,
j) ||
781 t_approx_P_adjoint__dstretch(
j,
i)) /
782 2.;
786 t_dP(L,
J) += t_L(
i,
j, L) * (t_diff2_uP(
i,
j,
k,
l) * t_L(
k,
l,
J));
787 }
788 ++t_approx_P_adjoint__dstretch;
789 ++t_log_u;
790 ++t_eigen_vals;
791 ++t_eigen_vecs;
792 ++t_nb_uniq;
793
794 int rr = 0;
795 for (; rr != row_nb_dofs /
size_symm; ++rr) {
798
799 auto t_m = get_ftensor2(
K, 6 * rr, 0);
800 for (
int cc = 0; cc != col_nb_dofs /
size_symm; ++cc) {
801 double b =
a * t_row_base_fun * t_col_base_fun;
802 double c = (
a * alpha_grad_u * ts_a) *
803 (t_row_grad_fun(
i) * t_col_grad_fun(
i));
804 t_m(L,
J) -= b * t_dP(L,
J);
805 t_m(L,
J) +=
c * t_kd_sym(L,
J);
806
807 ++t_m;
808 ++t_col_base_fun;
809 ++t_col_grad_fun;
810 }
811 ++t_row_base_fun;
812 ++t_row_grad_fun;
813 }
814
815 for (; rr != row_nb_base_functions; ++rr) {
816 ++t_row_base_fun;
817 ++t_row_grad_fun;
818 }
819
820 }
822 };
823
824 auto large = [&]() {
827 "Not implemented for Neo-Hookean (used ADOL-C)");
829 };
830
834 break;
838 break;
839 default:
841 "gradApproximator not handled");
842 };
843
845}
#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