v0.15.0
Loading...
Searching...
No Matches
HMHNeohookean.cpp
Go to the documentation of this file.
1/**
2 * @file HMHNeohookean.cpp
3 * @brief Implementation of NeoHookean material
4 * @date 2024-08-31
5 *
6 * @copyright Copyright (c) 2024
7 *
8 */
9
10namespace EshelbianPlasticity {
11
12// #define NEOHOOKEAN_SCALING
14
15 /**
16 * @brief Definition of Neo-hookean function
17 *
18 * @param c10
19 * @param v
20 * @return double
21 */
22 static inline double fun_neohookean(double c10, double v) {
23 return (2.0 * c10) * (-1. + EshelbianCore::f(2 * v));
24 }
25
26 /**
27 * @brief Definition of derivative of Neo-hookean function
28 *
29 * @param c10
30 * @param v
31 * @return * double
32 */
33 static inline double fun_d_neohookean(double c10, double v) {
34 return
35
36 (4.0 * c10) * EshelbianCore::d_f(2 * v);
37 }
38
39 /**
40 * @brief Definition of axiator of Neo-hookean function
41 *
42 * @param c10
43 * @param v
44 * @return double
45 */
46 static inline double fun_neohookean_bulk(double K, double tr) {
47 return K * tr; // K*ln(J)^2
48 }
49
50 /**
51 * @brief Definition of derivative of axiator of Neo-hookean function
52 *
53 * @param K
54 * @param tr
55 * @return double
56 */
57 static inline double fun_diff_neohookean_bulk(double K, double tr) {
58 return K;
59 }
60
61 static constexpr int numberOfActiveVariables = 9;
62 static constexpr int numberOfDependentVariables = 9;
63
64 HMHNeohookean(MoFEM::Interface &m_field, const double c10, const double K)
66 mField(m_field), c10_default(c10), K_default(K) {
67
68 CHK_THROW_MESSAGE(getOptions(), "get options failed");
70 "extract block data failed");
71
72#ifdef NEOHOOKEAN_SCALING
73 if (blockData.size()) {
74 double ave_K = 0;
75 for (auto b : blockData) {
76 ave_K += b.K;
77 }
78 eqScaling = ave_K / blockData.size();
79 }
80#endif
81
82 MOFEM_LOG("EP", Sev::inform) << "Neo-Hookean scaling = " << eqScaling;
83
84 if (adolCOn == PETSC_FALSE) {
87 "Stretch selector is not equal to LOG");
88 } else {
89 if (EshelbianCore::exponentBase != exp(1)) {
91 "Exponent base is not equal to exp(1)");
92 }
93 }
94 }
95 }
96
98 for (auto &b : blockData) {
99 if (b.blockEnts.find(ent) != b.blockEnts.end()) {
100 return std::make_pair(b.c10, b.K);
101 }
102 }
103 if (blockData.size() != 0)
105 "Block not found for entity handle. If you mat set "
106 "block, set it to all elements");
107 return std::make_pair(c10_default, K_default);
108 }
109
111 OpGetScale(boost::shared_ptr<double> scale_ptr,
112 boost::shared_ptr<PhysicalEquations> physics_ptr)
114 scalePtr(scale_ptr), physicsPtr(physics_ptr) {}
115
116 MoFEMErrorCode doWork(int side, EntityType type,
117 EntitiesFieldData::EntData &data) {
119 auto neohookean_ptr =
120 boost::dynamic_pointer_cast<HMHNeohookean>(physicsPtr);
121 *(scalePtr) = neohookean_ptr->eqScaling;
123 }
124
125 private:
126 boost::shared_ptr<double> scalePtr;
127 boost::shared_ptr<PhysicalEquations> physicsPtr;
128 };
129
131 returnOpSetScale(boost::shared_ptr<double> scale_ptr,
132 boost::shared_ptr<PhysicalEquations> physics_ptr) {
133
134 return new OpGetScale(scale_ptr, physics_ptr);
135 };
136
139 MoFEMErrorCode evaluateRhs(EntData &data) { return 0; }
140 MoFEMErrorCode evaluateLhs(EntData &data) { return 0; }
141 };
142
144 returnOpJacobian(const int tag, const bool eval_rhs, const bool eval_lhs,
145 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
146 boost::shared_ptr<PhysicalEquations> physics_ptr) {
147 if (adolCOn) {
148 return PhysicalEquations::returnOpJacobian(tag, eval_rhs, eval_lhs,
149 data_ptr, physics_ptr);
150 } else {
151 return (new OpJacobian(tag, eval_rhs, eval_lhs, data_ptr, physics_ptr));
152 }
153 }
154
155 MoFEMErrorCode getOptions() {
157 PetscOptionsBegin(PETSC_COMM_WORLD, "neo_hookean_", "", "none");
158
159 c10_default = 1;
160 CHKERR PetscOptionsScalar("-c10", "C10", "", c10_default, &c10_default,
161 PETSC_NULLPTR);
162 K_default = 1;
163 CHKERR PetscOptionsScalar("-K", "Bulk modulus K", "", K_default, &K_default,
164 PETSC_NULLPTR);
165
166 alphaGradU = 0;
167 CHKERR PetscOptionsScalar("-viscosity_alpha_grad_u", "viscosity", "",
168 alphaGradU, &alphaGradU, PETSC_NULLPTR);
169
170 CHKERR PetscOptionsBool("-adolc_on", "adolc ON", "", adolCOn, &adolCOn,
171 PETSC_NULLPTR);
172
173 PetscOptionsEnd();
174
175 MOFEM_LOG_CHANNEL("WORLD");
176 MOFEM_TAG_AND_LOG("WORLD", Sev::inform, "MatBlock Neo-Hookean (default)")
177 << "c10 = " << c10_default << " K = " << K_default
178 << " grad alpha u = " << alphaGradU;
179 MOFEM_TAG_AND_LOG("WORLD", Sev::inform, "MatBlock Neo-Hookean (default)")
180 << "ADOL-C On " << (adolCOn)
181 ? "Yes"
182 : "No";
184 }
185
186 MoFEMErrorCode extractBlockData(Sev sev) {
187 return extractBlockData(
188
189 mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
190
191 (boost::format("%s(.*)") % "MAT_NEOHOOKEAN").str()
192
193 )),
194
195 sev);
196 }
197
198 MoFEMErrorCode
199 extractBlockData(std::vector<const CubitMeshSets *> meshset_vec_ptr,
200 Sev sev) {
202
203 for (auto m : meshset_vec_ptr) {
204 MOFEM_LOG("EP", sev) << *m;
205 std::vector<double> block_data;
206 CHKERR m->getAttributes(block_data);
207 if (block_data.size() < 2) {
208 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
209 "Expected that block has atleast two attributes");
210 }
211 auto get_block_ents = [&]() {
212 Range ents;
213 CHKERR mField.get_moab().get_entities_by_handle(m->meshset, ents, true);
214 return ents;
215 };
216
217 double c10 = block_data[0];
218 double K = block_data[1];
219
220 blockData.push_back({c10, K, get_block_ents()});
221
222 MOFEM_LOG("EP", sev) << "MatBlock Neo-Hookean c10 = "
223 << blockData.back().c10
224 << " K = " << blockData.back().K << " nb ents. = "
225 << blockData.back().blockEnts.size();
226 }
228 }
229
230 MoFEMErrorCode recordTape(const int tape, DTensor2Ptr *t_h_ptr) {
232
233 FTENSOR_INDEX(3, i);
234 FTENSOR_INDEX(3, j);
235 FTENSOR_INDEX(3, I);
236 FTENSOR_INDEX(3, J);
237
238 auto large = [&]() {
240
241 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
242
243 auto ih = get_h();
244 if (t_h_ptr)
245 ih(i, j) = (*t_h_ptr)(i, j);
246 else {
247 ih(i, j) = t_kd(i, j);
248 }
249
250 auto r_P = get_P();
251
252 enableMinMaxUsingAbs();
253 trace_on(tape);
254
255 // Set active variables to ADOL-C
256 th(i, j) <<= get_h()(i, j);
257
258 tF(i, I) = th(i, I);
259 CHKERR determinantTensor3by3(tF, detF);
260 CHKERR invertTensor3by3(tF, detF, tInvF);
261
262 tCof(i, I) = detF * tInvF(I, i);
263
264 tP(i, I) = 2. * c10_default * (tF(i, I) - tInvF(i, I)) +
265 K_default * log(detF) * tCof(i, I);
266
267 // Set dependent variables to ADOL-C
268 tP(i, j) >>= r_P(i, j);
269
270 trace_off();
271
273 };
274
277 case LARGE_ROT:
278 case MODERATE_ROT:
279 CHKERR large();
280 break;
281 default:
282 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
283 "gradApproximator not handled");
284 };
285
287 }
288
290
291 OpSpatialPhysical(const std::string &field_name,
292 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
293 const double alpha_u);
294
295 MoFEMErrorCode integrate(EntData &data);
296
297 private:
298 const double alphaU;
299 };
300
301 virtual VolUserDataOperator *
303 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
304 const double alpha_u) {
305 if (adolCOn) {
307 alpha_u);
308 } else {
309 return new OpSpatialPhysical(field_name, data_ptr, alpha_u);
310 }
311 }
312
314
316 const std::string &field_name,
317 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
318 boost::shared_ptr<ExternalStrainVec> external_strain_vec_ptr,
319 std::map<std::string, boost::shared_ptr<ScalingMethod>> smv);
320
321 MoFEMErrorCode integrate(EntData &data);
322
323 private:
324 boost::shared_ptr<ExternalStrainVec> externalStrainVecPtr;
325 std::map<std::string, boost::shared_ptr<ScalingMethod>> scalingMethodsMap;
326 };
327
329 const std::string &field_name,
330 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
331 boost::shared_ptr<ExternalStrainVec> external_strain_vec_ptr,
332 std::map<std::string, boost::shared_ptr<ScalingMethod>> smv) {
333 return new OpSpatialPhysicalExternalStrain(field_name, data_ptr,
334 external_strain_vec_ptr, smv);
335 }
336
338 OpSpatialPhysical_du_du(std::string row_field, std::string col_field,
339 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
340 const double alpha);
341 MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
342
343 private:
344 const double alphaU;
345 };
346
348 std::string row_field, std::string col_field,
349 boost::shared_ptr<DataAtIntegrationPts> data_ptr, const double alpha) {
350 if (adolCOn) {
352 row_field, col_field, data_ptr, alpha);
353 } else {
354 return new OpSpatialPhysical_du_du(row_field, col_field, data_ptr, alpha);
355 }
356 }
357
358private:
360
362 double K_default;
364 struct BlockData {
365 double c10;
366 double K;
368 };
369 std::vector<BlockData> blockData;
370
371 double eqScaling = 1.;
372
373 // ADOL-C
374 PetscBool adolCOn = PETSC_FALSE;
375
382};
383
385 const std::string &field_name,
386 boost::shared_ptr<DataAtIntegrationPts> data_ptr, const double alpha_u)
387 : OpAssembleVolume(field_name, data_ptr, OPROW), alphaU(alpha_u) {}
388
389extern "C" {
390void tetcircumcenter_tp(double a[3], double b[3], double c[3], double d[3],
391 double circumcenter[3], double *xi, double *eta,
392 double *zeta);
393}
394
397
398 auto neohookean_ptr =
399 boost::dynamic_pointer_cast<HMHNeohookean>(dataAtPts->physicsPtr);
400 if (!neohookean_ptr) {
401 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
402 "Pointer to HMHNeohookean is null");
403 }
404 auto [def_c10, def_K] =
405 neohookean_ptr->getMaterialParameters(getFEEntityHandle());
406
407 double c10 = def_c10 / neohookean_ptr->eqScaling;
408 double alpha_u = alphaU / neohookean_ptr->eqScaling;
409 double K = def_K / neohookean_ptr->eqScaling;
410
411 auto time_step = getTStimeStep();
412 double alpha_grad_u =
413 neohookean_ptr->alphaGradU / neohookean_ptr->eqScaling;// / time_step;
414
416 auto t_L = symm_L_tensor();
417
418 constexpr auto t_kd_sym = FTensor::Kronecker_Delta_symmetric<int>();
419 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
420
421 int nb_dofs = data.getIndices().size();
422 int nb_integration_pts = data.getN().size1();
423 auto v = getVolume();
424 auto t_w = getFTensor0IntegrationWeight();
425 auto t_approx_P_adjoint_log_du =
426 getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
427 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
428 auto t_grad_h1 = getFTensor2FromMat<3, 3>(dataAtPts->wGradH1AtPts);
429 auto t_dot_log_u =
430 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
431 auto t_diff_u =
432 getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
433 auto t_log_u =
434 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
435 auto t_grad_log_u =
436 getFTensor2FromMat<6, 3>(dataAtPts->gradLogStretchDotTensorAtPts);
437 auto t_log_u2_h1 =
438 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretch2H1AtPts);
439
440 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
441 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
442 auto &nbUniq = dataAtPts->nbUniq;
443 auto t_nb_uniq =
444 FTensor::Tensor0<FTensor::PackPtr<int *, 1>>(nbUniq.data().data());
445
446 auto t_diff = diff_tensor();
447
454
455 auto get_ftensor2 = [](auto &v) {
457 &v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
458 };
459
460 int nb_base_functions = data.getN().size2();
461 auto t_row_base_fun = data.getFTensor0N();
462 auto t_grad_base_fun = data.getFTensor1DiffN<3>();
463
464 auto no_h1 = [&]() {
466
467 for (int gg = 0; gg != nb_integration_pts; ++gg) {
468 double a = v * t_w;
469 ++t_w;
470
471 auto neohookean = [c10](auto v) { return fun_neohookean(c10, v); };
472 auto t_neohookean_hencky =
473 EigenMatrix::getMat(t_eigen_vals, t_eigen_vecs, neohookean);
474 ++t_eigen_vals;
475 ++t_eigen_vecs;
476 ++t_nb_uniq;
477
478 const double tr = t_log_u(i, i);
480 t_P(L) = t_L(i, j, L) * (t_neohookean_hencky(i, j) +
481 t_kd(i, j) * fun_neohookean_bulk(K, tr));
483 t_viscous_P(L) = alpha_u * (t_L(i, j, L) * t_dot_log_u(i, j));
484
486 t_residual(L) = t_approx_P_adjoint_log_du(L) - t_P(L) - t_viscous_P(L);
487 t_residual(L) *= a;
488
490 t_grad_residual(L, i) = alpha_grad_u * t_grad_log_u(L, i);
491 t_grad_residual(L, i) *= a;
492
493 ++t_approx_P_adjoint_log_du;
494 ++t_log_u;
495 ++t_dot_log_u;
496 ++t_grad_log_u;
497
498 auto t_nf = getFTensor1FromPtr<size_symm>(&*nF.data().begin());
499 int bb = 0;
500 for (; bb != nb_dofs / size_symm; ++bb) {
501 t_nf(L) -= t_row_base_fun * t_residual(L);
502 t_nf(L) += t_grad_base_fun(i) * t_grad_residual(L, i);
503 ++t_nf;
504 ++t_row_base_fun;
505 ++t_grad_base_fun;
506 }
507 for (; bb != nb_base_functions; ++bb) {
508 ++t_row_base_fun;
509 ++t_grad_base_fun;
510 }
511 }
512
514 };
515
516 auto large = [&]() {
518 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
519 "Not implemented for Neo-Hookean (used ADOL-C)");
521 };
522
525 CHKERR no_h1();
526 break;
527 case LARGE_ROT:
528 case MODERATE_ROT:
529 CHKERR large();
530 break;
531 default:
532 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
533 "gradApproximator not handled");
534 };
535
537}
538
540 const std::string &field_name,
541 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
542 boost::shared_ptr<ExternalStrainVec> external_strain_vec_ptr,
543 std::map<std::string, boost::shared_ptr<ScalingMethod>> smv)
544 : OpAssembleVolume(field_name, data_ptr, OPROW),
545 externalStrainVecPtr(external_strain_vec_ptr), scalingMethodsMap{smv} {}
546
547MoFEMErrorCode
550
551 auto neohookean_ptr =
552 boost::dynamic_pointer_cast<HMHNeohookean>(dataAtPts->physicsPtr);
553 if (!neohookean_ptr) {
554 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
555 "Pointer to HMHNeohookean is null");
556 }
557
558 double time = OpAssembleVolume::getFEMethod()->ts_t;
561 }
562 // get entity of tet
564 // iterate over all block data
565
566 for (auto &ext_strain_block : (*externalStrainVecPtr)) {
567 // check if finite element entity is part of the EXTERNALSTRAIN block
568 if (ext_strain_block.ents.find(fe_ent) != ext_strain_block.ents.end()) {
569
570 double scale = 1;
571 if (scalingMethodsMap.find(ext_strain_block.blockName) !=
572 scalingMethodsMap.end()) {
573 scale *=
574 scalingMethodsMap.at(ext_strain_block.blockName)->getScale(time);
575 } else {
576 MOFEM_LOG("SELF", Sev::warning)
577 << "No scaling method found for " << ext_strain_block.blockName;
578 }
579
580 auto [def_c10, def_K] =
581 neohookean_ptr->getMaterialParameters(getFEEntityHandle());
582 double c10 = def_c10 / neohookean_ptr->eqScaling;
583 // get ExternalStrain block data
584 double external_strain_val = scale * ext_strain_block.val;
585 double K = ext_strain_block.bulkModulusK;
586
588 auto t_L = symm_L_tensor();
589 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
590
591 int nb_dofs = data.getIndices().size();
592 int nb_integration_pts = data.getN().size1();
593 auto v = getVolume();
594 auto t_w = getFTensor0IntegrationWeight();
595
602
603 int nb_base_functions = data.getN().size2();
604 auto t_row_base_fun = data.getFTensor0N();
605
606 const double tr = 3.0 * external_strain_val;
607 const double sigma_J = K * tr;
608
609 for (int gg = 0; gg != nb_integration_pts; ++gg) {
610 double a = v * t_w;
611 ++t_w;
612
614 t_residual(L) = 0.0;
615 t_residual(L) += (t_L(i, j, L) * t_kd(i, j)) * sigma_J;
616 t_residual(L) *= a;
617
618 auto t_nf = getFTensor1FromPtr<size_symm>(&*nF.data().begin());
619 int bb = 0;
620 for (; bb != nb_dofs / size_symm; ++bb) {
621 t_nf(L) -= t_row_base_fun * t_residual(L);
622 ++t_nf;
623 ++t_row_base_fun;
624 }
625 for (; bb != nb_base_functions; ++bb) {
626 ++t_row_base_fun;
627 }
628 }
629 }
630 }
632}
633
635 std::string row_field, std::string col_field,
636 boost::shared_ptr<DataAtIntegrationPts> data_ptr, const double alpha)
637 : OpAssembleVolumePositiveDefine(row_field, col_field, data_ptr, OPROWCOL,
638 false),
639 alphaU(alpha) {
640 sYmm = false;
641}
642
643MoFEMErrorCode
645 EntData &col_data) {
647
648 auto neohookean_ptr =
649 boost::dynamic_pointer_cast<HMHNeohookean>(dataAtPts->physicsPtr);
650 if (!neohookean_ptr) {
651 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
652 "Pointer to HMHNeohookean is null");
653 }
654 auto [def_c10, def_K] =
655 neohookean_ptr->getMaterialParameters(getFEEntityHandle());
656
657 double c10 = def_c10 / neohookean_ptr->eqScaling;
658 double alpha_u = alphaU / neohookean_ptr->eqScaling;
659 double lambda = def_K / neohookean_ptr->eqScaling;
660
661 auto time_step = getTStimeStep();
662 double alpha_grad_u =
663 neohookean_ptr->alphaGradU / neohookean_ptr->eqScaling;// / time_step;
664
667
668 constexpr auto t_kd_sym = FTensor::Kronecker_Delta_symmetric<int>();
669 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
670
671 auto t_L = symm_L_tensor();
672 auto t_diff = diff_tensor();
673
674 int nb_integration_pts = row_data.getN().size1();
675 int row_nb_dofs = row_data.getIndices().size();
676 int col_nb_dofs = col_data.getIndices().size();
677
678 auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
680 size_symm>(
681
682 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 0, c + 3),
683 &m(r + 0, c + 4), &m(r + 0, c + 5),
684
685 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 1, c + 3),
686 &m(r + 1, c + 4), &m(r + 1, c + 5),
687
688 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2), &m(r + 2, c + 3),
689 &m(r + 2, c + 4), &m(r + 2, c + 5),
690
691 &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2), &m(r + 3, c + 3),
692 &m(r + 3, c + 4), &m(r + 3, c + 5),
693
694 &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2), &m(r + 4, c + 3),
695 &m(r + 4, c + 4), &m(r + 4, c + 5),
696
697 &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2), &m(r + 5, c + 3),
698 &m(r + 5, c + 4), &m(r + 5, c + 5)
699
700 );
701 };
702
709
710 auto v = getVolume();
711 auto ts_a = getTSa();
712 auto t_w = getFTensor0IntegrationWeight();
713
714 int row_nb_base_functions = row_data.getN().size2();
715 auto t_row_base_fun = row_data.getFTensor0N();
716 auto t_row_grad_fun = row_data.getFTensor1DiffN<3>();
717
718 auto t_grad_h1 = getFTensor2FromMat<3, 3>(dataAtPts->wGradH1AtPts);
719 auto t_diff_u =
720 getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
721 auto t_log_u =
722 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
723 auto t_log_u2_h1 =
724 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretch2H1AtPts);
725 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
726 auto t_approx_P_adjoint__dstretch =
727 getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
728 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
729 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
730 auto &nbUniq = dataAtPts->nbUniq;
731 auto t_nb_uniq =
732 FTensor::Tensor0<FTensor::PackPtr<int *, 1>>(nbUniq.data().data());
733
734 auto no_h1 = [&]() {
736
737 for (int gg = 0; gg != nb_integration_pts; ++gg) {
738 double a = v * t_w;
739 ++t_w;
740
741 auto neohookean = [c10](auto v) { return fun_neohookean(c10, v); };
742 auto d_neohookean = [c10, lambda](auto v) {
743 return fun_d_neohookean(c10, v);
744 };
745
746 auto t_diff_neohookean = EigenMatrix::getDiffMat(
747 t_eigen_vals, t_eigen_vecs, neohookean, d_neohookean, t_nb_uniq);
748
749 const auto tr = t_log_u(i, i);
751 t_dP(L, J) = -t_L(i, j, L) * ((t_diff_neohookean(i, j, k, l) +
753 t_kd_sym(i, j) * t_kd_sym(k, l)) *
754 t_L(k, l, J));
755 t_dP(L, J) -= (alpha_u * ts_a) *
756 (t_L(i, j, L) * (t_diff(i, j, k, l) * t_L(k, l, J)));
757
758 if constexpr (1) {
760 t_deltaP(i, j) = (t_approx_P_adjoint__dstretch(i, j) ||
761 t_approx_P_adjoint__dstretch(j, i)) /
762 2.;
763 auto t_diff2_uP = EigenMatrix::getDiffDiffMat(
764 t_eigen_vals, t_eigen_vecs, EshelbianCore::f, EshelbianCore::d_f,
765 EshelbianCore::dd_f, t_deltaP, t_nb_uniq);
766 t_dP(L, J) += t_L(i, j, L) * (t_diff2_uP(i, j, k, l) * t_L(k, l, J));
767 }
768 ++t_approx_P_adjoint__dstretch;
769 ++t_log_u;
770 ++t_eigen_vals;
771 ++t_eigen_vecs;
772 ++t_nb_uniq;
773
774 int rr = 0;
775 for (; rr != row_nb_dofs / size_symm; ++rr) {
776 auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
777 auto t_col_grad_fun = col_data.getFTensor1DiffN<3>(gg, 0);
778
779 auto t_m = get_ftensor2(K, 6 * rr, 0);
780 for (int cc = 0; cc != col_nb_dofs / size_symm; ++cc) {
781 double b = a * t_row_base_fun * t_col_base_fun;
782 double c = (a * alpha_grad_u * ts_a) *
783 (t_row_grad_fun(i) * t_col_grad_fun(i));
784 t_m(L, J) -= b * t_dP(L, J);
785 t_m(L, J) += c * t_kd_sym(L, J);
786
787 ++t_m;
788 ++t_col_base_fun;
789 ++t_col_grad_fun;
790 }
791 ++t_row_base_fun;
792 ++t_row_grad_fun;
793 }
794
795 for (; rr != row_nb_base_functions; ++rr) {
796 ++t_row_base_fun;
797 ++t_row_grad_fun;
798 }
799
800 }
802 };
803
804 auto large = [&]() {
806 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
807 "Not implemented for Neo-Hookean (used ADOL-C)");
809 };
810
813 CHKERR no_h1();
814 break;
815 case LARGE_ROT:
816 case MODERATE_ROT:
817 CHKERR large();
818 break;
819 default:
820 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
821 "gradApproximator not handled");
822 };
823
825}
826
827} // namespace EshelbianPlasticity
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
#define FTENSOR_INDEX(DIM, I)
constexpr double a
constexpr int SPACE_DIM
Kronecker Delta class symmetric.
Kronecker Delta class.
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
@ NOSPACE
Definition definitions.h:83
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
constexpr auto t_kd
double eta
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
FTensor::Index< 'i', SPACE_DIM > i
static double lambda
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
Definition level_set.cpp:30
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
auto getMat(A &&t_val, B &&t_vec, Fun< double > f)
Get the Mat object.
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.
void tetcircumcenter_tp(double a[3], double b[3], double c[3], double d[3], double circumcenter[3], double *xi, double *eta, double *zeta)
static constexpr auto size_symm
constexpr IntegrationType I
double scale
Definition plastic.cpp:123
double zeta
Viscous hardening.
Definition plastic.cpp:130
constexpr auto field_name
FTensor::Index< 'm', 3 > m
static enum StretchSelector stretchSelector
static boost::function< double(const double)> f
static boost::function< double(const double)> dd_f
static boost::function< double(const double)> d_f
OpGetScale(boost::shared_ptr< double > scale_ptr, boost::shared_ptr< PhysicalEquations > physics_ptr)
boost::shared_ptr< PhysicalEquations > physicsPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
std::map< std::string, boost::shared_ptr< ScalingMethod > > scalingMethodsMap
OpSpatialPhysicalExternalStrain(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< ExternalStrainVec > external_strain_vec_ptr, std::map< std::string, boost::shared_ptr< ScalingMethod > > smv)
OpSpatialPhysical_du_du(std::string row_field, std::string col_field, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const double alpha)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
OpSpatialPhysical(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const double alpha_u)
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_bulk(double K, double tr)
Definition of axiator of Neo-hookean function.
virtual VolUserDataOperator * returnOpSpatialPhysical(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const double alpha_u)
std::vector< BlockData > blockData
MoFEMErrorCode extractBlockData(std::vector< const CubitMeshSets * > meshset_vec_ptr, Sev sev)
MoFEMErrorCode extractBlockData(Sev sev)
static constexpr int numberOfDependentVariables
VolUserDataOperator * returnOpSetScale(boost::shared_ptr< double > scale_ptr, boost::shared_ptr< PhysicalEquations > physics_ptr)
static double fun_neohookean(double c10, double v)
Definition of Neo-hookean function.
auto getMaterialParameters(EntityHandle ent)
static constexpr int numberOfActiveVariables
virtual VolUserDataOperator * returnOpSpatialPhysicalExternalStrain(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< ExternalStrainVec > external_strain_vec_ptr, std::map< std::string, boost::shared_ptr< ScalingMethod > > smv)
HMHNeohookean(MoFEM::Interface &m_field, const double c10, const double K)
VolUserDataOperator * returnOpSpatialPhysical_du_du(std::string row_field, std::string col_field, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const double alpha)
MoFEMErrorCode recordTape(const int tape, DTensor2Ptr *t_h_ptr)
UserDataOperator * returnOpJacobian(const int tag, const bool eval_rhs, const bool eval_lhs, boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< PhysicalEquations > physics_ptr)
virtual VolUserDataOperator * returnOpSpatialPhysical(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const double alpha_u)
virtual VolUserDataOperator * returnOpSpatialPhysical_du_du(std::string row_field, std::string col_field, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const double alpha)
virtual UserDataOperator * returnOpJacobian(const int tag, const bool eval_rhs, const bool eval_lhs, boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< PhysicalEquations > physics_ptr)
virtual moab::Interface & get_moab()=0
bool sYmm
If true assume that matrix is symmetric structure.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
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 dofs on entity.
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
PetscReal ts_t
time
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.