29 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
30 boost::shared_ptr<PhysicalEquations> physics_ptr) {
31 return (
new OpJacobian(tag, eval_rhs, eval_lhs, data_ptr, physics_ptr));
36 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"neo_hookean_",
"",
"none");
39 CHKERR PetscOptionsScalar(
"-c10",
"C10",
"",
c10, &
c10, PETSC_NULL);
41 CHKERR PetscOptionsScalar(
"-K",
"Bulk modulus K",
"",
K, &
K, PETSC_NULL);
43 ierr = PetscOptionsEnd();
60 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
61 const double alpha_u);
71 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
72 const double alpha_u) {
78 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
87 std::string row_field, std::string col_field,
88 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
const double alpha) {
96 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
const double alpha_u)
102 auto neohookean_ptr =
103 boost::dynamic_pointer_cast<HMHNeohookean>(dataAtPts->physicsPtr);
104 if (!neohookean_ptr) {
106 "Pointer to HMHNeohookean is null");
112 "Stretch selector is not equal to LOG");
116 "Exponent base is not equal to exp(1)");
121 const auto c10 = neohookean_ptr->c10;
122 const auto K = neohookean_ptr->K;
131 int nb_integration_pts = data.
getN().size1();
132 auto v = getVolume();
133 auto t_w = getFTensor0IntegrationWeight();
134 auto t_approx_P_adjont_log_du =
135 getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
136 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
137 auto t_grad_h1 = getFTensor2FromMat<3, 3>(dataAtPts->wGradH1AtPts);
139 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
141 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
143 getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
145 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
147 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretch2H1AtPts);
158 auto get_ftensor2 = [](
auto &
v) {
160 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
163 int nb_base_functions = data.
getN().size2();
165 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
176 t_Ldiff_u(
i,
j,
L) = t_diff_u(
i,
j,
k,
l) * t_L(
k,
l,
L);
182 t_Ldiff_u(
i,
j,
L) = (t_diff_u(
i,
m,
k,
l) * t_h1(
m,
j)) * t_L(
k,
l,
L);
187 t_Ldiff_u(
i,
j,
L) = t_diff_u(
i,
j,
k,
l) * t_L(
k,
l,
L);
194 t_Sigma_u(
i,
j) = 2.0 *
c10 * (t_u(
i,
m) * t_h1(
m,
j));
195 const double tr = t_total_log_u(
i,
j) * t_kd_sym(
i,
j);
197 const double Simga_J = -2 *
c10 +
K * (
J - 1) *
J;
206 t_residual(
L) = t_approx_P_adjont_log_du(
L);
207 t_residual(
L) -= t_Ldiff_u(
i,
j,
L) * t_Sigma_u(
i,
j);
208 t_residual(
L) -= (t_L(
i,
j,
L) *
t_kd(
i,
j)) * Simga_J;
209 t_residual(
L) -= t_L(
i,
j,
L) * t_viscous_P(
i,
j);
212 ++t_approx_P_adjont_log_du;
217 auto t_nf = getFTensor1FromPtr<size_symm>(&*nF.data().begin());
219 for (; bb != nb_dofs /
size_symm; ++bb) {
220 t_nf(
L) += t_row_base_fun * t_residual(
L);
224 for (; bb != nb_base_functions; ++bb)
231 std::string row_field, std::string col_field,
232 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
const double alpha)
243 auto neohookean_ptr =
244 boost::dynamic_pointer_cast<HMHNeohookean>(dataAtPts->physicsPtr);
245 if (!neohookean_ptr) {
247 "Pointer to HMHNeohookean is null");
253 "Stretch selector is not equal to LOG");
257 "Exponent base is not equal to exp(1)");
262 const auto c10 = neohookean_ptr->c10;
263 const auto lambda = neohookean_ptr->K;
274 int nb_integration_pts = row_data.
getN().size1();
275 int row_nb_dofs = row_data.
getIndices().size();
276 int col_nb_dofs = col_data.
getIndices().size();
282 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2), &
m(
r + 0,
c + 3),
283 &
m(
r + 0,
c + 4), &
m(
r + 0,
c + 5),
285 &
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2), &
m(
r + 1,
c + 3),
286 &
m(
r + 1,
c + 4), &
m(
r + 1,
c + 5),
288 &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2), &
m(
r + 2,
c + 3),
289 &
m(
r + 2,
c + 4), &
m(
r + 2,
c + 5),
291 &
m(
r + 3,
c + 0), &
m(
r + 3,
c + 1), &
m(
r + 3,
c + 2), &
m(
r + 3,
c + 3),
292 &
m(
r + 3,
c + 4), &
m(
r + 3,
c + 5),
294 &
m(
r + 4,
c + 0), &
m(
r + 4,
c + 1), &
m(
r + 4,
c + 2), &
m(
r + 4,
c + 3),
295 &
m(
r + 4,
c + 4), &
m(
r + 4,
c + 5),
297 &
m(
r + 5,
c + 0), &
m(
r + 5,
c + 1), &
m(
r + 5,
c + 2), &
m(
r + 5,
c + 3),
298 &
m(
r + 5,
c + 4), &
m(
r + 5,
c + 5)
310 auto v = getVolume();
311 auto ts_a = getTSa();
312 auto t_w = getFTensor0IntegrationWeight();
314 int row_nb_base_functions = row_data.
getN().size2();
317 auto t_grad_h1 = getFTensor2FromMat<3, 3>(dataAtPts->wGradH1AtPts);
319 getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
321 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
322 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
323 auto t_approx_P_adjont_dstretch =
324 getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
325 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
326 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
327 auto &nbUniq = dataAtPts->nbUniq;
329 auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->PAtPts);
330 auto r_P_du = getFTensor4FromMat<3, 3, 3, 3>(dataAtPts->P_du);
332 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
342 t_Ldiff_u(
i,
j,
L) = t_diff_u(
i,
j,
k,
l) * t_L(
k,
l,
L);
346 t_Ldiff_u(
i,
j,
L) = (t_diff_u(
i,
m,
k,
l) * t_h1(
m,
j)) * t_L(
k,
l,
L);
349 t_Ldiff_u(
i,
j,
L) = t_diff_u(
i,
j,
k,
l) * t_L(
k,
l,
L);
355 const double tr = t_total_log_u(
i,
j) * t_kd_sym(
i,
j);
363 t_dP(
L,
J) = (-2.0 *
c10) * (t_Ldiff_u(
i,
j,
L) * t_Ldiff_u(
i,
j,
J));
364 t_dP(
L,
J) -= (t_L(
i,
j,
L) *
t_kd(
i,
j)) * Simga_J_dtr(
J);
366 (alphaU * ts_a) * (t_L(
i,
j,
L) * (t_diff(
i,
j,
k,
l) * t_L(
k,
l,
J)));
369 t_Sigma_u(
i,
j) = 2.0 *
c10 * (t_u(
i,
m) * t_h1(
m,
j));
374 t_approx_P_adjont_dstretch(
i,
j) - t_h1(
j,
n) * t_Sigma_u(
i,
n);
375 ++t_approx_P_adjont_dstretch;
379 t_deltaP_sym(
i,
j) = (t_deltaP(
i,
j) || t_deltaP(
j,
i));
380 t_deltaP_sym(
i,
j) /= 2.0;
384 t_dP(
L,
J) += t_L(
i,
j,
L) * (t_diff2_uP2(
i,
j,
k,
l) * t_L(
k,
l,
J));
390 for (; rr != row_nb_dofs /
size_symm; ++rr) {
392 auto t_m = get_ftensor2(
K, 6 * rr, 0);
393 for (
int cc = 0; cc != col_nb_dofs /
size_symm; ++cc) {
394 double b =
a * t_row_base_fun * t_col_base_fun;
395 t_m(
L,
J) += b * t_dP(
L,
J);
402 for (; rr != row_nb_base_functions; ++rr) {