9 #include <BasicFiniteElements.hpp>
12 #include <boost/math/constants/constants.hpp>
14 constexpr
double third = boost::math::constants::third<double>();
19 OpHMHH(
const std::string &field_name,
const int tag,
const bool eval_rhs,
20 const bool eval_lhs, boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
21 boost::shared_ptr<PhysicalEquations> &physics_ptr)
22 :
OpJacobian(field_name, tag, eval_rhs, eval_lhs, data_ptr, physics_ptr) {
35 const auto nb_integration_pts = ent_data.
getN().size1();
36 auto iu = getFTensor2SymmetricFromMat<3>(dataAtPts->streachTensorAtPts);
38 auto iG = getFTensor2FromMat<3, 3>(dataAtPts->GAtPts);
40 auto create_data_vec = [nb_integration_pts](
auto &v) {
41 v.resize(nb_integration_pts,
false);
43 auto create_data_mat = [nb_integration_pts](
auto &
m) {
44 m.resize(9, nb_integration_pts,
false);
47 create_data_mat(dataAtPts->PAtPts);
48 create_data_mat(dataAtPts->SigmaAtPts);
49 create_data_vec(dataAtPts->phiAtPts);
50 create_data_mat(dataAtPts->flowAtPts);
52 auto r_P = getFTensor2FromMat<3, 3>(dataAtPts->PAtPts);
53 auto r_Sigma = getFTensor2FromMat<3, 3>(dataAtPts->SigmaAtPts);
55 auto r_flow = getFTensor2FromMat<3, 3>(dataAtPts->flowAtPts);
57 for (
unsigned int gg = 0; gg != nb_integration_pts; ++gg) {
59 auto t_h = physicsPtr->get_h();
60 for (
auto ii : {0, 1, 2})
61 for (
auto jj : {0, 1, 2})
62 t_h(ii, jj) = iu(ii, jj);
64 physicsPtr->get_H()(
i,
j) = iG(
i,
j);
65 for (
int ii = 0; ii != 3; ++ii)
66 physicsPtr->get_H()(ii, ii) += 1;
69 int r = ::function(tAg, physicsPtr->dependentVariables.size(),
70 physicsPtr->activeVariables.size(),
71 &physicsPtr->activeVariables[0],
72 &physicsPtr->dependentVariables[0]);
75 "ADOL-C function evaluation with error r = %d",
r);
78 r_P(
i,
j) = physicsPtr->get_P()(
i,
j);
79 r_Sigma(
i,
j) = physicsPtr->get_Sigma()(
i,
j);
80 r_phi = physicsPtr->get_Phi();
81 r_flow(
i,
j) = physicsPtr->get_Flow()(
i,
j);
98 const int number_of_active_variables = physicsPtr->activeVariables.size();
99 const int number_of_dependent_variables =
100 physicsPtr->dependentVariables.size();
101 std::vector<double *> jac_ptr(number_of_dependent_variables);
102 for (
unsigned int n = 0;
n != number_of_dependent_variables; ++
n) {
105 ->dependentVariablesDirevatives[
n * number_of_active_variables]);
108 const auto nb_integration_pts = ent_data.
getN().size1();
110 auto create_data_mat = [nb_integration_pts](
auto &
m) {
111 m.resize(9, nb_integration_pts,
false);
114 auto create_data_ten = [nb_integration_pts](
auto &
m) {
115 m.resize(27, nb_integration_pts,
false);
118 create_data_ten(dataAtPts->P_dh0);
119 create_data_ten(dataAtPts->P_dh1);
120 create_data_ten(dataAtPts->P_dh2);
121 create_data_ten(dataAtPts->P_dH0);
122 create_data_ten(dataAtPts->P_dH1);
123 create_data_ten(dataAtPts->P_dH2);
124 create_data_ten(dataAtPts->Sigma_dh0);
125 create_data_ten(dataAtPts->Sigma_dh1);
126 create_data_ten(dataAtPts->Sigma_dh2);
127 create_data_ten(dataAtPts->Sigma_dH0);
128 create_data_ten(dataAtPts->Sigma_dH1);
129 create_data_ten(dataAtPts->Sigma_dH2);
130 create_data_mat(dataAtPts->phi_dh);
131 create_data_mat(dataAtPts->phi_dH);
132 create_data_ten(dataAtPts->Flow_dh0);
133 create_data_ten(dataAtPts->Flow_dh1);
134 create_data_ten(dataAtPts->Flow_dh2);
135 create_data_ten(dataAtPts->Flow_dH0);
136 create_data_ten(dataAtPts->Flow_dH1);
137 create_data_ten(dataAtPts->Flow_dH2);
139 auto iu = getFTensor2SymmetricFromMat<3>(dataAtPts->streachTensorAtPts);
140 auto iG = getFTensor2FromMat<3, 3>(dataAtPts->GAtPts);
154 auto r_phi_dh = getFTensor2FromMat<3, 3>(dataAtPts->phi_dh);
155 auto r_phi_dH = getFTensor2FromMat<3, 3>(dataAtPts->phi_dH);
163 for (
unsigned int gg = 0; gg != nb_integration_pts; ++gg) {
165 auto t_h = physicsPtr->get_h();
166 for(
auto ii : {0, 1, 2})
167 for(
auto jj : {0, 1, 2})
168 t_h(ii, jj) = iu(ii, jj);
170 physicsPtr->get_H()(
i,
j) = iG(
i,
j);
171 for (
int ii = 0; ii != 3; ++ii)
172 physicsPtr->get_H()(ii, ii) += 1;
176 int r = ::jacobian(tAg, number_of_dependent_variables,
177 number_of_active_variables,
178 &physicsPtr->activeVariables[0], &jac_ptr[0]);
181 "ADOL-C function evaluation with error");
184 r_P_dh0(
i,
j,
k) = physicsPtr->get_P_dh0()(
i,
j,
k);
185 r_P_dh1(
i,
j,
k) = physicsPtr->get_P_dh1()(
i,
j,
k);
186 r_P_dh2(
i,
j,
k) = physicsPtr->get_P_dh2()(
i,
j,
k);
187 r_P_dH0(
i,
j,
k) = physicsPtr->get_P_dH0()(
i,
j,
k);
188 r_P_dH1(
i,
j,
k) = physicsPtr->get_P_dH1()(
i,
j,
k);
189 r_P_dH2(
i,
j,
k) = physicsPtr->get_P_dH2()(
i,
j,
k);
190 r_Sigma_dh0(
i,
j,
k) = physicsPtr->get_Sigma_dh0()(
i,
j,
k);
191 r_Sigma_dh1(
i,
j,
k) = physicsPtr->get_Sigma_dh1()(
i,
j,
k);
192 r_Sigma_dh2(
i,
j,
k) = physicsPtr->get_Sigma_dh2()(
i,
j,
k);
193 r_Sigma_dH0(
i,
j,
k) = physicsPtr->get_Sigma_dH0()(
i,
j,
k);
194 r_Sigma_dH1(
i,
j,
k) = physicsPtr->get_Sigma_dH1()(
i,
j,
k);
195 r_Sigma_dH2(
i,
j,
k) = physicsPtr->get_Sigma_dH2()(
i,
j,
k);
196 r_phi_dh(
i,
j) = physicsPtr->get_Phi_dh()(
i,
j);
197 r_phi_dH(
i,
j) = physicsPtr->get_Phi_dH()(
i,
j);
198 r_Flow_dh0(
i,
j,
k) = physicsPtr->get_Flow_dh0()(
i,
j,
k);
199 r_Flow_dh1(
i,
j,
k) = physicsPtr->get_Flow_dh1()(
i,
j,
k);
200 r_Flow_dh2(
i,
j,
k) = physicsPtr->get_Flow_dh2()(
i,
j,
k);
201 r_Flow_dH0(
i,
j,
k) = physicsPtr->get_Flow_dH0()(
i,
j,
k);
202 r_Flow_dH1(
i,
j,
k) = physicsPtr->get_Flow_dH1()(
i,
j,
k);
203 r_Flow_dH2(
i,
j,
k) = physicsPtr->get_Flow_dH2()(
i,
j,
k);
233 static constexpr
int numberOfActiveVariables = 9 + 9;
234 static constexpr
int numberOfDependentVariables = 9 + 9 + 1 + 9;
237 const double sigma_y)
239 lambda(lambda), mu(mu),
sigmaY(sigma_y) {}
247 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"stvenant_",
"",
"none");
249 CHKERR PetscOptionsScalar(
"-young_modulus",
"Young modulus",
"",
E, &
E,
251 CHKERR PetscOptionsScalar(
"-poisson_ratio",
"poisson ratio",
"", nu, &nu,
256 ierr = PetscOptionsEnd();
267 const bool eval_rhs,
const bool eval_lhs,
268 boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
269 boost::shared_ptr<PhysicalEquations> &physics_ptr) {
271 new OpHMHH(field_name, tag, eval_rhs, eval_lhs, data_ptr, physics_ptr));
323 ih(
i,
j) = (*t_h_ptr)(
i,
j);
326 for (
auto ii : {0, 1, 2})
331 for (
auto ii : {0, 1, 2})
335 auto r_Sigma = get_Sigma();
336 double &r_phi = get_PhiRef();
337 auto r_Flow = get_Flow();
339 enableMinMaxUsingAbs();
343 th(
i,
j) <<= get_h()(
i,
j);
344 tH(
i,
j) <<= get_H()(
i,
j);
349 tF(
i,
I) = th(
i,
J) * tInvH(
J,
I);
354 tC(
I,
J) = tF(
i,
I) * tF(
i,
J);
363 energy = 0.5 * lambda * trE * trE + mu * (tE(
I,
J) * tE(
I,
J));
368 tS(
N0,
N0) += lambda * trE;
369 tS(
N1,
N1) += lambda * trE;
370 tS(
N2,
N2) += lambda * trE;
372 tP(
i,
J) = tF(
i,
I) * tS(
I,
J);
374 tCauchy(
i,
j) = tP(
i,
J) * tF(
j,
J);
375 tCauchy(
i,
j) *= (1 / detF);
378 tSigma(
I,
J) = -tF(
i,
I) * tP(
i,
J);
379 tSigma(
N0,
N0) += energy;
380 tSigma(
N1,
N1) += energy;
381 tSigma(
N2,
N2) += energy;
384 meanCauchy =
third * tCauchy(
i,
i);
385 tDevCauchy(
i,
j) = tCauchy(
i,
j);
386 tDevCauchy(
N0,
N0) -= meanCauchy;
387 tDevCauchy(
N1,
N1) -= meanCauchy;
388 tDevCauchy(
N2,
N2) -= meanCauchy;
392 s2 = tDevCauchy(
i,
j) * tDevCauchy(
i,
j);
397 tPhi_dDevCauchy(
i,
j) = tDevCauchy(
i,
j);
398 tPhi_dDevCauchy(
i,
j) *= 1.5 /
f;
399 tPhi_dCauchy(
i,
j) = tPhi_dDevCauchy(
i,
j);
400 tPhi_dCauchy(
N0,
N0) -=
third * tPhi_dDevCauchy(
i,
i);
401 tPhi_dCauchy(
N1,
N1) -=
third * tPhi_dDevCauchy(
i,
i);
402 tPhi_dCauchy(
N2,
N2) -=
third * tPhi_dDevCauchy(
i,
i);
403 tPhi_dSigma(
I,
J) = tInvF(
I,
i) * (tPhi_dCauchy(
i,
j) * tF(
j,
J));
404 tPhi_dSigma(
I,
J) *= -(1 / detF);
407 tPulledP(
i,
J) = tP(
i,
I) * tInvH(
J,
I);
408 tPulledP(
i,
J) *= detH;
409 tPulledSigma(
i,
J) = tSigma(
i,
I) * tInvH(
J,
I);
410 tPulledSigma(
i,
J) *= detH;
411 tPulledPhi_dSigma(
i,
J) = tPhi_dSigma(
i,
I) * tInvH(
J,
I);
412 tPulledPhi_dSigma(
i,
J) *= detH;
415 tPulledP(
i,
j) >>= r_P(
i,
j);
416 tPulledSigma(
i,
j) >>= r_Sigma(
i,
j);
418 tPulledPhi_dSigma(
i,
j) >>= r_Flow(
i,
j);
428 static constexpr
int numberOfActiveVariables = 9 + 9;
429 static constexpr
int numberOfDependentVariables = 9 + 9 + 1 + 9;
432 const double lambda,
const double sigma_y)
438 const bool eval_rhs,
const bool eval_lhs,
439 boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
440 boost::shared_ptr<PhysicalEquations> &physics_ptr) {
442 new OpHMHH(field_name, tag, eval_rhs, eval_lhs, data_ptr, physics_ptr));
447 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"mooneyrivlin_",
"",
"none");
450 CHKERR PetscOptionsScalar(
"-alpha",
"Alpha",
"",
alpha, &
alpha, PETSC_NULL);
452 CHKERR PetscOptionsScalar(
"-beta",
"Beta",
"", beta, &beta, PETSC_NULL);
455 CHKERR PetscOptionsScalar(
"-lambda",
"Lambda",
"", lambda, &lambda,
459 CHKERR PetscOptionsScalar(
"-epsilon",
"Epsilon",
"", epsilon, &epsilon,
465 ierr = PetscOptionsEnd();
531 ih(
i,
j) = (*t_h_ptr)(
i,
j);
534 for (
auto ii : {0, 1, 2})
539 for (
auto ii : {0, 1, 2})
543 auto r_Sigma = get_Sigma();
544 double &r_phi = get_PhiRef();
545 auto r_Flow = get_Flow();
547 enableMinMaxUsingAbs();
551 th(
i,
j) <<= get_h()(
i,
j);
552 tH(
i,
j) <<= get_H()(
i,
j);
558 tF(
i,
I) = th(
i,
J) * tInvH(
J,
I);
562 tCof(
i,
I) = detF * tInvF(
I,
i);
564 A = tF(
k,
K) * tF(
k,
K);
565 B = tCof(
k,
K) * tCof(
k,
K);
568 tBCof(
i,
I) = 4 * beta * (
B * tCof(
i,
I));
569 tBj = (-12 *
alpha - 24 * beta) / detF +
570 0.5 * (lambda / epsilon) *
571 (pow(detF, epsilon - 1) - pow(detF, -epsilon - 1));
573 tP(
i,
I) = tBF(
i,
I);
576 tP(
i,
I) += tCof(
i,
I) * tBj;
578 energy =
alpha * pow(A, 2) + beta * pow(
B, 2);
579 energy += (-12 *
alpha - 24 * beta) * log(detF) +
580 (lambda / (2 * epsilon * epsilon)) *
581 (pow(detF, epsilon) + pow(detF, -epsilon));
584 tSigma(
I,
J) = -tF(
i,
I) * tP(
i,
J);
585 tSigma(
N0,
N0) += energy;
586 tSigma(
N1,
N1) += energy;
587 tSigma(
N2,
N2) += energy;
590 meanCauchy =
third * tCauchy(
i,
i);
591 tDevCauchy(
i,
j) = tCauchy(
i,
j);
592 tDevCauchy(
N0,
N0) -= meanCauchy;
593 tDevCauchy(
N1,
N1) -= meanCauchy;
594 tDevCauchy(
N2,
N2) -= meanCauchy;
598 s2 = tDevCauchy(
i,
j) * tDevCauchy(
i,
j);
603 tPhi_dDevCauchy(
i,
j) = tDevCauchy(
i,
j);
604 tPhi_dDevCauchy(
i,
j) *= 1.5 /
f;
605 tPhi_dCauchy(
i,
j) = tPhi_dDevCauchy(
i,
j);
606 tPhi_dCauchy(
N0,
N0) -=
third * tPhi_dDevCauchy(
i,
i);
607 tPhi_dCauchy(
N1,
N1) -=
third * tPhi_dDevCauchy(
i,
i);
608 tPhi_dCauchy(
N2,
N2) -=
third * tPhi_dDevCauchy(
i,
i);
609 tPhi_dSigma(
I,
J) = tInvF(
I,
i) * (tPhi_dCauchy(
i,
j) * tF(
j,
J));
610 tPhi_dSigma(
I,
J) *= -(1 / detF);
613 tPulledP(
i,
J) = tP(
i,
I) * tInvH(
J,
I);
614 tPulledP(
i,
J) *= detH;
615 tPulledSigma(
i,
J) = tSigma(
i,
I) * tInvH(
J,
I);
616 tPulledSigma(
i,
J) *= detH;
617 tPulledPhi_dSigma(
i,
J) = tPhi_dSigma(
i,
I) * tInvH(
J,
I);
618 tPulledPhi_dSigma(
i,
J) *= detH;
621 tPulledP(
i,
j) >>= r_P(
i,
j);
622 tPulledSigma(
i,
j) >>= r_Sigma(
i,
j);
624 tPulledPhi_dSigma(
i,
j) >>= r_Flow(
i,
j);
633 const int tape,
const double lambda,
const double mu,
634 const double sigma_y) {
636 physicalEquations = boost::shared_ptr<HMHStVenantKirchhoff>(
638 CHKERR physicalEquations->recordTape(tape,
nullptr);
643 const int tape,
const double alpha,
const double beta,
const double lambda,
644 const double sigma_y) {
646 physicalEquations = boost::shared_ptr<HMHPMooneyRivlinWriggersEq63>(
648 CHKERR physicalEquations->recordTape(tape,
nullptr);
Eshelbian plasticity interface.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
@ MOFEM_OPERATION_UNSUCCESSFUL
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
FTensor::Index< 'k', 3 > k
FTensor::Tensor3< FTensor::PackPtr< double *, 1 >, 3, 3, 3 > getFTensor3FromMat(MatrixDouble &m)
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
static FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det)
Calculate determinant 3 by 3.
const double r
rate factor
virtual OpJacobian * returnOpJacobian(const std::string &field_name, const int tag, const bool eval_rhs, const bool eval_lhs, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, boost::shared_ptr< PhysicalEquations > &physics_ptr)
MoFEMErrorCode recordTape(const int tape, DTensor2Ptr *t_h_ptr)
MoFEMErrorCode getOptions()
ATensor2 tPulledPhi_dSigma
HMHPMooneyRivlinWriggersEq63(const double alpha, const double beta, const double lambda, const double sigma_y)
HMHStVenantKirchhoff(const double lambda, const double mu, const double sigma_y)
virtual OpJacobian * returnOpJacobian(const std::string &field_name, const int tag, const bool eval_rhs, const bool eval_lhs, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, boost::shared_ptr< PhysicalEquations > &physics_ptr)
MoFEMErrorCode recordTape(const int tape, DTensor2Ptr *t_h_ptr)
MoFEMErrorCode getOptions()
ATensor2 tPulledPhi_dSigma
OpHMHH(const std::string &field_name, const int tag, const bool eval_rhs, const bool eval_lhs, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, boost::shared_ptr< PhysicalEquations > &physics_ptr)
Data on single entity (This is passed as argument to DataOperator::doWork)
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....