v0.15.0
Loading...
Searching...
No Matches
HMHHencky.cpp
Go to the documentation of this file.
1/**
2 * @file Hencky.cpp
3 * @brief Implementation of Hencky material
4 * @date 2024-08-31
5 *
6 * @copyright Copyright (c) 2024
7 *
8 */
9
10namespace EshelbianPlasticity {
11
13
14 HMHHencky(MoFEM::Interface &m_field, const double E, const double nu)
15 : PhysicalEquations(0, 0), mField(m_field), E(E), nu(nu) {}
16
17 MoFEMErrorCode recordTape(const int tag, DTensor2Ptr *t_h) { return 0; }
18
19 struct OpHenckyJacobian : public OpJacobian {
20 OpHenckyJacobian(boost::shared_ptr<DataAtIntegrationPts> data_ptr,
21 boost::shared_ptr<HMHHencky> hencky_ptr)
22 : OpJacobian(), dataAtGaussPts(data_ptr), henckyPtr(hencky_ptr) {
24 "getOptions failed");
25 CHK_THROW_MESSAGE(henckyPtr->extractBlockData(Sev::verbose),
26 "Can not get data from block");
27 }
28
29 MoFEMErrorCode doWork(int side, EntityType type,
30 EntitiesFieldData::EntData &data) {
32
33 for (auto &b : henckyPtr->blockData) {
34 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
35 dataAtGaussPts->mu = b.shearModulusG;
36 dataAtGaussPts->lambda = b.bulkModulusK - 2 * b.shearModulusG / 3;
37 CHKERR henckyPtr->getMatDPtr(dataAtGaussPts->getMatDPtr(),
38 dataAtGaussPts->getMatAxiatorDPtr(),
39 dataAtGaussPts->getMatDeviatorDPtr(),
40 b.bulkModulusK, b.shearModulusG);
42 }
43 }
44 const auto E = henckyPtr->E;
45 const auto nu = henckyPtr->nu;
46
47 double bulk_modulus_K = E / (3 * (1 - 2 * nu));
48 double shear_modulus_G = E / (2 * (1 + nu));
49
52
53 CHKERR henckyPtr->getMatDPtr(dataAtGaussPts->getMatDPtr(),
54 dataAtGaussPts->getMatAxiatorDPtr(),
55 dataAtGaussPts->getMatDeviatorDPtr(),
57
59 }
60
61 MoFEMErrorCode evaluateRhs(EntData &data) { return 0; }
62 MoFEMErrorCode evaluateLhs(EntData &data) { return 0; }
63
64 private:
65 boost::shared_ptr<DataAtIntegrationPts> dataAtGaussPts;
66 boost::shared_ptr<HMHHencky> henckyPtr;
67 };
68
69 virtual UserDataOperator *
70 returnOpJacobian(const int tag, const bool eval_rhs, const bool eval_lhs,
71 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
72 boost::shared_ptr<PhysicalEquations> physics_ptr) {
73 return (new OpHenckyJacobian(
74 data_ptr, boost::dynamic_pointer_cast<HMHHencky>(physics_ptr)));
75 }
76
78
79 OpSpatialPhysical(const std::string &field_name,
80 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
81 const double alpha_u);
82
83 MoFEMErrorCode integrate(EntData &data);
84
85 MoFEMErrorCode integrateHencky(EntData &data);
86
87 MoFEMErrorCode integratePolyconvexHencky(EntData &data);
88
89 private:
90 const double alphaU;
91 PetscBool polyConvex = PETSC_FALSE;
92 };
93
94 virtual VolUserDataOperator *
96 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
97 const double alpha_u) {
98 return new OpSpatialPhysical(field_name, data_ptr, alpha_u);
99 }
100
103 const std::string &field_name,
104 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
105 boost::shared_ptr<ExternalStrainVec> &external_strain_vec_ptr,
106 std::map<std::string, boost::shared_ptr<ScalingMethod>> smv);
107
108 MoFEMErrorCode integrate(EntData &data);
109
110 private:
111 boost::shared_ptr<ExternalStrainVec> externalStrainVecPtr;
112 std::map<std::string, boost::shared_ptr<ScalingMethod>> scalingMethodsMap;
113 };
114
116 const std::string &field_name,
117 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
118 boost::shared_ptr<ExternalStrainVec> external_strain_vec_ptr,
119 std::map<std::string, boost::shared_ptr<ScalingMethod>> smv)
120 {
122 field_name, data_ptr, external_strain_vec_ptr, smv);
123 }
124
126 const double alphaU;
127 OpSpatialPhysical_du_du(std::string row_field, std::string col_field,
128 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
129 const double alpha);
130 MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
131 MoFEMErrorCode integrateHencky(EntData &row_data, EntData &col_data);
132 MoFEMErrorCode integratePolyconvexHencky(EntData &row_data,
133 EntData &col_data);
134
135 private:
136 PetscBool polyConvex = PETSC_FALSE;
137
138 MatrixDouble dP;
139 };
140
142 std::string row_field, std::string col_field,
143 boost::shared_ptr<DataAtIntegrationPts> data_ptr, const double alpha) {
144 return new OpSpatialPhysical_du_du(row_field, col_field, data_ptr, alpha);
145 }
146
147 /**
148 * @brief Calculate energy density for Hencky material model
149 *
150 *
151 * \f[
152 *
153 * \Psi(\log{\mathbf{U}}) = \frac{1}{2} U_{IJ} D_{IJKL} U_{KL} = \frac{1}{2}
154 * U_{IJ} T_{IJ}
155 *
156 * \f]
157 * where \f$T_{IJ} = D_{IJKL} U_{KL}\f$ is a a Hencky stress.
158 *
159 */
161
162 OpCalculateEnergy(boost::shared_ptr<DataAtIntegrationPts> data_ptr,
163 boost::shared_ptr<double> total_energy_ptr);
164 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
165
166 private:
167 boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
168 boost::shared_ptr<double> totalEnergyPtr;
169 };
170
172 returnOpCalculateEnergy(boost::shared_ptr<DataAtIntegrationPts> data_ptr,
173 boost::shared_ptr<double> total_energy_ptr) {
174 return new OpCalculateEnergy(data_ptr, total_energy_ptr);
175 }
176
179 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
180 boost::shared_ptr<MatrixDouble> strain_ptr,
181 boost::shared_ptr<MatrixDouble> stress_ptr,
182 boost::shared_ptr<HMHHencky> hencky_ptr);
183 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
184
185 private:
186 boost::shared_ptr<DataAtIntegrationPts>
187 dataAtPts; ///< data at integration pts
188 boost::shared_ptr<MatrixDouble> strainPtr;
189 boost::shared_ptr<MatrixDouble> stressPtr;
190 boost::shared_ptr<HMHHencky> henckyPtr;
191 };
192
194 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
195 boost::shared_ptr<PhysicalEquations> physics_ptr) {
197 data_ptr, data_ptr->getLogStretchTensorAtPts(),
198 data_ptr->getApproxPAtPts(),
199 boost::dynamic_pointer_cast<HMHHencky>(physics_ptr));
200 }
201
203 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
204 boost::shared_ptr<PhysicalEquations> physics_ptr) {
206 data_ptr, data_ptr->getVarLogStreachPts(), data_ptr->getVarPiolaPts(),
207 boost::dynamic_pointer_cast<HMHHencky>(physics_ptr));
208 }
209
210 MoFEMErrorCode getOptions(boost::shared_ptr<DataAtIntegrationPts> data_ptr) {
212 PetscOptionsBegin(PETSC_COMM_WORLD, "hencky_", "", "none");
213
214 CHKERR PetscOptionsScalar("-young_modulus", "Young modulus", "", E, &E,
215 PETSC_NULLPTR);
216 CHKERR PetscOptionsScalar("-poisson_ratio", "poisson ratio", "", nu, &nu,
217 PETSC_NULLPTR);
218
219 PetscOptionsEnd();
220
222 << "Hencky: E = " << E << " nu = " << nu;
223 getOptionsSeverityLevels = Sev::verbose;
224
225 CHKERRG(ierr);
226
228 }
229
230 MoFEMErrorCode extractBlockData(Sev sev) {
231 return extractBlockData(
232
233 mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
234
235 (boost::format("%s(.*)") % "MAT_ELASTIC").str()
236
237 )),
238
239 sev);
240 }
241
242 MoFEMErrorCode
243 extractBlockData(std::vector<const CubitMeshSets *> meshset_vec_ptr,
244 Sev sev) {
246
247 for (auto m : meshset_vec_ptr) {
248 MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock") << *m;
249 std::vector<double> block_data;
250 CHKERR m->getAttributes(block_data);
251 if (block_data.size() < 2) {
252 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
253 "Expected that block has atleast two attributes");
254 }
255 auto get_block_ents = [&]() {
256 Range ents;
257 CHKERR mField.get_moab().get_entities_by_handle(m->meshset, ents, true);
258 return ents;
259 };
260
261 double young_modulus = block_data[0];
262 double poisson_ratio = block_data[1];
263 double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
264 double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
265
266 MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock")
267 << "E = " << young_modulus << " nu = " << poisson_ratio;
268
269 blockData.push_back({bulk_modulus_K, shear_modulus_G, get_block_ents()});
270 }
271 MOFEM_LOG_CHANNEL("WORLD");
273 }
274
275 MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
276 boost::shared_ptr<MatrixDouble> mat_axiator_D_ptr,
277 boost::shared_ptr<MatrixDouble> mat_deviator_D_ptr,
278 double bulk_modulus_K, double shear_modulus_G) {
280 //! [Calculate elasticity tensor]
281 auto set_material_stiffness = [&]() {
287 auto t_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(
288 &*(mat_D_ptr->data().begin()));
289 auto t_axiator_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(
290 &*mat_axiator_D_ptr->data().begin());
291 auto t_deviator_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(
292 &*mat_deviator_D_ptr->data().begin());
293 t_axiator_D(i, j, k, l) = (bulk_modulus_K - (2. / 3.) * shear_modulus_G) *
294 t_kd(i, j) * t_kd(k, l);
295 t_deviator_D(i, j, k, l) =
296 2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.);
297 t_D(i, j, k, l) = t_axiator_D(i, j, k, l) + t_deviator_D(i, j, k, l);
298 };
299 //! [Calculate elasticity tensor]
300 constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
301 mat_D_ptr->resize(size_symm, size_symm, false);
302 mat_axiator_D_ptr->resize(size_symm, size_symm, false);
303 mat_deviator_D_ptr->resize(size_symm, size_symm, false);
304 set_material_stiffness();
306 }
307
308 MoFEMErrorCode
309 getInvMatDPtr(boost::shared_ptr<MatrixDouble> mat_inv_D_ptr,
310 double bulk_modulus_K, double shear_modulus_G) {
312 //! [Calculate elasticity tensor]
313 auto set_material_compilance = [&]() {
319 const double A = 1. / (2. * shear_modulus_G);
320 const double B =
321 (1. / (9. * bulk_modulus_K)) - (1. / (6. * shear_modulus_G));
322 auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(
323 &*(mat_inv_D_ptr->data().begin()));
324 t_inv_D(i, j, k, l) =
325 A * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) + B * t_kd(i, j) * t_kd(k, l);
326 };
327 //! [Calculate elasticity tensor]
328 constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
329 mat_inv_D_ptr->resize(size_symm, size_symm, false);
330 set_material_compilance();
332 }
333
334private:
336
342 std::vector<BlockData> blockData;
343
344 double E;
345 double nu;
346
347 // Set verbosity level, it verbile can changes that informatno is pronated
348 // only once at particular level
349 Sev getOptionsSeverityLevels = Sev::inform;
350};
351
353 const std::string &field_name,
354 boost::shared_ptr<DataAtIntegrationPts> data_ptr, const double alpha_u)
355 : OpAssembleVolume(field_name, data_ptr, OPROW), alphaU(alpha_u) {
356
357 CHK_MOAB_THROW(PetscOptionsGetBool(PETSC_NULLPTR, "", "-poly_convex",
358 &polyConvex, PETSC_NULLPTR),
359 "get polyconvex option failed");
360}
361
364 if (polyConvex) {
365 CHKERR integratePolyconvexHencky(data);
366 } else {
367 CHKERR integrateHencky(data);
368 }
370}
371
374
376 auto t_L = symm_L_tensor();
377
378 int nb_dofs = data.getIndices().size();
379 int nb_integration_pts = data.getN().size1();
380 auto v = getVolume();
381 auto t_w = getFTensor0IntegrationWeight();
382 auto t_approx_P_adjoint_log_du =
383 getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
384 auto t_log_stretch_h1 =
385 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
386 auto t_dot_log_u =
387 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
388
389 auto t_D = getFTensor4DdgFromPtr<3, 3, 0>(&*dataAtPts->matD.data().begin());
390
391 FTensor::Index<'i', 3> i;
392 FTensor::Index<'j', 3> j;
393 FTensor::Index<'k', 3> k;
394 FTensor::Index<'l', 3> l;
395 auto get_ftensor2 = [](auto &v) {
397 &v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
398 };
399
400 int nb_base_functions = data.getN().size2();
401 auto t_row_base_fun = data.getFTensor0N();
402 for (int gg = 0; gg != nb_integration_pts; ++gg) {
403 double a = v * t_w;
404 auto t_nf = get_ftensor2(nF);
405
407 t_T(i, j) =
408 t_D(i, j, k, l) * (t_log_stretch_h1(k, l) + alphaU * t_dot_log_u(k, l));
410 t_residual(L) =
411 a * (t_approx_P_adjoint_log_du(L) - t_L(i, j, L) * t_T(i, j));
412
413 int bb = 0;
414 for (; bb != nb_dofs / 6; ++bb) {
415 t_nf(L) -= t_row_base_fun * t_residual(L);
416 ++t_nf;
417 ++t_row_base_fun;
418 }
419 for (; bb != nb_base_functions; ++bb)
420 ++t_row_base_fun;
421
422 ++t_w;
423 ++t_approx_P_adjoint_log_du;
424 ++t_dot_log_u;
425 ++t_log_stretch_h1;
426 }
428}
429
430MoFEMErrorCode
433
435 auto t_L = symm_L_tensor();
436
437 int nb_dofs = data.getIndices().size();
438 int nb_integration_pts = data.getN().size1();
439 auto v = getVolume();
440 auto t_w = getFTensor0IntegrationWeight();
441 auto t_approx_P_adjoint_log_du =
442 getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
443 auto t_log_stretch_h1 =
444 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
445 auto t_dot_log_u =
446 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
447
448 auto t_D = getFTensor4DdgFromPtr<3, 3, 0>(&*dataAtPts->matD.data().begin());
449
450 FTensor::Index<'i', 3> i;
451 FTensor::Index<'j', 3> j;
452 FTensor::Index<'k', 3> k;
453 FTensor::Index<'l', 3> l;
454 auto get_ftensor2 = [](auto &v) {
456 &v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
457 };
458
459 constexpr double nohat_k = 1. / 4;
460 constexpr double hat_k = 1. / 8;
461 double mu = dataAtPts->mu;
462 double lambda = dataAtPts->lambda;
463
464 constexpr double third = boost::math::constants::third<double>();
466 auto t_diff_deviator = diff_deviator(diff_tensor());
467
468 int nb_base_functions = data.getN().size2();
469 auto t_row_base_fun = data.getFTensor0N();
470 for (int gg = 0; gg != nb_integration_pts; ++gg) {
471 double a = v * t_w;
472 auto t_nf = get_ftensor2(nF);
473
474 double log_det = t_log_stretch_h1(i, i);
475 double log_det2 = log_det * log_det;
477 t_dev(i, j) = t_log_stretch_h1(i, j) - t_kd(i, j) * (third * log_det);
478 double dev_norm2 = t_dev(i, j) * t_dev(i, j);
479
481 auto A = 2 * mu * std::exp(nohat_k * dev_norm2);
482 auto B = lambda * std::exp(hat_k * log_det2) * log_det;
483 t_T(i, j) =
484
485 A * (t_dev(k, l) * t_diff_deviator(k, l, i, j))
486
487 +
488
489 B * t_kd(i, j)
490
491 +
492
493 alphaU * t_D(i, j, k, l) * t_dot_log_u(k, l);
494
496 t_residual(L) =
497 a * (t_approx_P_adjoint_log_du(L) - t_L(i, j, L) * t_T(i, j));
498
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;
503 ++t_row_base_fun;
504 }
505 for (; bb != nb_base_functions; ++bb)
506 ++t_row_base_fun;
507
508 ++t_w;
509 ++t_approx_P_adjoint_log_du;
510 ++t_dot_log_u;
511 ++t_log_stretch_h1;
512 }
514}
515
517 std::string row_field, std::string col_field,
518 boost::shared_ptr<DataAtIntegrationPts> data_ptr, const double alpha)
519 : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL, false),
520 alphaU(alpha) {
521 sYmm = false;
522
523 CHK_MOAB_THROW(PetscOptionsGetBool(PETSC_NULLPTR, "", "-poly_convex",
524 &polyConvex, PETSC_NULLPTR),
525 "get polyconvex option failed");
526}
527
529 const std::string &field_name,
530 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
531 boost::shared_ptr<ExternalStrainVec> &external_strain_vec_ptr,
532 std::map<std::string, boost::shared_ptr<ScalingMethod>> smv)
533 : OpAssembleVolume(field_name, data_ptr, OPROW),
534 externalStrainVecPtr(external_strain_vec_ptr), scalingMethodsMap(smv) {}
535
536MoFEMErrorCode
539
541
542 double time = OpAssembleVolume::getFEMethod()->ts_t;
545 }
546 // get entity of tet
548 // iterate over all block data
549
550 for (auto &ext_strain_block : (*externalStrainVecPtr)) {
551 // check if finite element entity is part of the EXTERNALSTRAIN block
552 if (ext_strain_block.ents.find(fe_ent) != ext_strain_block.ents.end()) {
553
554 double scale = 1;
555 if (scalingMethodsMap.find(ext_strain_block.blockName) !=
556 scalingMethodsMap.end()) {
557 scale *=
558 scalingMethodsMap.at(ext_strain_block.blockName)->getScale(time);
559 } else {
560 MOFEM_LOG("SELF", Sev::warning)
561 << "No scaling method found for " << ext_strain_block.blockName;
562 }
563
564 // get ExternalStrain data
565 double external_strain_val = scale * ext_strain_block.val;
566 double bulk_modulus_K = ext_strain_block.bulkModulusK;
567 auto t_L = symm_L_tensor();
568
569 int nb_dofs = data.getIndices().size();
570 int nb_integration_pts = data.getN().size1();
571 auto v = getVolume();
572 auto t_w = getFTensor0IntegrationWeight();
573
575
576 FTensor::Index<'i', 3> i;
577 FTensor::Index<'j', 3> j;
578 FTensor::Index<'k', 3> k;
579 FTensor::Index<'l', 3> l;
580 auto get_ftensor2 = [](auto &v) {
582 &v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
583 };
584
585 int nb_base_functions = data.getN().size2();
586 auto t_row_base_fun = data.getFTensor0N();
587 for (int gg = 0; gg != nb_integration_pts; ++gg) {
588 auto tr = 3.0 * external_strain_val;
589 double a = v * t_w;
590 auto t_nf = get_ftensor2(nF);
591
593
594 t_T(i, j) = -bulk_modulus_K * tr * t_kd(i, j);
595
597 t_residual(L) = a * (t_L(i, j, L) * t_T(i, j));
598
599 int bb = 0;
600 for (; bb != nb_dofs / 6; ++bb) {
601 t_nf(L) += t_row_base_fun * t_residual(L);
602 ++t_nf;
603 ++t_row_base_fun;
604 }
605 for (; bb != nb_base_functions; ++bb)
606 ++t_row_base_fun;
607
608 ++t_w;
609 }
610 }
611 }
612
614}
615
616MoFEMErrorCode
618 EntData &col_data) {
620
621 if (polyConvex) {
622 CHKERR integratePolyconvexHencky(row_data, col_data);
623 } else {
624 CHKERR integrateHencky(row_data, col_data);
625 }
627}
628
629MoFEMErrorCode
631 EntData &col_data) {
633
636 auto t_L = symm_L_tensor();
637 auto t_diff = diff_tensor();
638
639 int nb_integration_pts = row_data.getN().size1();
640 int row_nb_dofs = row_data.getIndices().size();
641 int col_nb_dofs = col_data.getIndices().size();
642
643 auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
645 size_symm>(
646
647 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 0, c + 3),
648 &m(r + 0, c + 4), &m(r + 0, c + 5),
649
650 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 1, c + 3),
651 &m(r + 1, c + 4), &m(r + 1, c + 5),
652
653 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2), &m(r + 2, c + 3),
654 &m(r + 2, c + 4), &m(r + 2, c + 5),
655
656 &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2), &m(r + 3, c + 3),
657 &m(r + 3, c + 4), &m(r + 3, c + 5),
658
659 &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2), &m(r + 4, c + 3),
660 &m(r + 4, c + 4), &m(r + 4, c + 5),
661
662 &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2), &m(r + 5, c + 3),
663 &m(r + 5, c + 4), &m(r + 5, c + 5)
664
665 );
666 };
667 FTensor::Index<'i', 3> i;
668 FTensor::Index<'j', 3> j;
669 FTensor::Index<'k', 3> k;
670 FTensor::Index<'l', 3> l;
671 FTensor::Index<'m', 3> m;
672 FTensor::Index<'n', 3> n;
673
674 auto v = getVolume();
675 auto t_w = getFTensor0IntegrationWeight();
676
677 auto t_approx_P_adjoint__dstretch =
678 getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
679 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
680 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
681
682 int row_nb_base_functions = row_data.getN().size2();
683 auto t_row_base_fun = row_data.getFTensor0N();
684
685 auto get_dP = [&]() {
686 dP.resize(size_symm * size_symm, nb_integration_pts, false);
687 auto ts_a = getTSa();
688
689 auto t_D = getFTensor4DdgFromPtr<3, 3, 0>(&*dataAtPts->matD.data().begin());
691 t_dP_tmp(L, J) = -(1 + alphaU * ts_a) *
692 (t_L(i, j, L) *
693 ((t_D(i, j, m, n) * t_diff(m, n, k, l)) * t_L(k, l, J)));
694
697 auto t_approx_P_adjoint__dstretch =
698 getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
699 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
700 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
701 auto &nbUniq = dataAtPts->nbUniq;
702
703 auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
704 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
705
706 // Work of symmetric tensor on undefined tensor is equal to the work
707 // of the symmetric part of it
709 t_sym(i, j) = (t_approx_P_adjoint__dstretch(i, j) ||
710 t_approx_P_adjoint__dstretch(j, i));
711 t_sym(i, j) /= 2.0;
712 auto t_diff2_uP2 = EigenMatrix::getDiffDiffMat(
713 t_eigen_vals, t_eigen_vecs, EshelbianCore::f, EshelbianCore::d_f,
714 EshelbianCore::dd_f, t_sym, nbUniq[gg]);
715 t_dP(L, J) = t_L(i, j, L) *
716 ((t_diff2_uP2(i, j, k, l) + t_diff2_uP2(k, l, i, j)) *
717 t_L(k, l, J)) /
718 2. +
719 t_dP_tmp(L, J);
720
721 ++t_dP;
722 ++t_approx_P_adjoint__dstretch;
723 ++t_eigen_vals;
724 ++t_eigen_vecs;
725 }
726 } else {
727 auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
728 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
729 t_dP(L, J) = t_dP_tmp(L, J);
730 ++t_dP;
731 }
732 }
733
734 return getFTensor2FromMat<size_symm, size_symm>(dP);
735 };
736
737 auto t_dP = get_dP();
738 for (int gg = 0; gg != nb_integration_pts; ++gg) {
739 double a = v * t_w;
740
741 int rr = 0;
742 for (; rr != row_nb_dofs / 6; ++rr) {
743 auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
744 auto t_m = get_ftensor2(K, 6 * rr, 0);
745 for (int cc = 0; cc != col_nb_dofs / 6; ++cc) {
746 const double b = a * t_row_base_fun * t_col_base_fun;
747 t_m(L, J) -= b * t_dP(L, J);
748 ++t_m;
749 ++t_col_base_fun;
750 }
751 ++t_row_base_fun;
752 }
753
754 for (; rr != row_nb_base_functions; ++rr) {
755 ++t_row_base_fun;
756 }
757
758 ++t_w;
759 ++t_dP;
760 }
762}
763
765 EntData &row_data, EntData &col_data) {
767
770 auto t_L = symm_L_tensor();
771 auto t_diff = diff_tensor();
772
773 int nb_integration_pts = row_data.getN().size1();
774 int row_nb_dofs = row_data.getIndices().size();
775 int col_nb_dofs = col_data.getIndices().size();
776
777 auto get_ftensor2 = [](MatrixDouble &m, const int r, const int c) {
779 size_symm>(
780
781 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 0, c + 3),
782 &m(r + 0, c + 4), &m(r + 0, c + 5),
783
784 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 1, c + 3),
785 &m(r + 1, c + 4), &m(r + 1, c + 5),
786
787 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2), &m(r + 2, c + 3),
788 &m(r + 2, c + 4), &m(r + 2, c + 5),
789
790 &m(r + 3, c + 0), &m(r + 3, c + 1), &m(r + 3, c + 2), &m(r + 3, c + 3),
791 &m(r + 3, c + 4), &m(r + 3, c + 5),
792
793 &m(r + 4, c + 0), &m(r + 4, c + 1), &m(r + 4, c + 2), &m(r + 4, c + 3),
794 &m(r + 4, c + 4), &m(r + 4, c + 5),
795
796 &m(r + 5, c + 0), &m(r + 5, c + 1), &m(r + 5, c + 2), &m(r + 5, c + 3),
797 &m(r + 5, c + 4), &m(r + 5, c + 5)
798
799 );
800 };
801 FTensor::Index<'i', 3> i;
802 FTensor::Index<'j', 3> j;
803 FTensor::Index<'k', 3> k;
804 FTensor::Index<'l', 3> l;
805 FTensor::Index<'m', 3> m;
806 FTensor::Index<'n', 3> n;
807
808 auto v = getVolume();
809 auto t_w = getFTensor0IntegrationWeight();
810
811 int row_nb_base_functions = row_data.getN().size2();
812 auto t_row_base_fun = row_data.getFTensor0N();
813
814 auto get_dP = [&]() {
815 dP.resize(size_symm * size_symm, nb_integration_pts, false);
816 auto ts_a = getTSa();
817
818 auto t_D = getFTensor4DdgFromPtr<3, 3, 0>(&*dataAtPts->matD.data().begin());
819
820 constexpr double nohat_k = 1. / 4;
821 constexpr double hat_k = 1. / 8;
822 double mu = dataAtPts->mu;
823 double lambda = dataAtPts->lambda;
824
825 constexpr double third = boost::math::constants::third<double>();
827 auto t_diff_deviator = diff_deviator(diff_tensor());
828
829 auto t_approx_P_adjoint__dstretch =
830 getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
831 auto t_log_stretch_h1 =
832 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
833 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
834 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
835 auto &nbUniq = dataAtPts->nbUniq;
836
837 auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
838 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
839
840 double log_det = t_log_stretch_h1(i, i);
841 double log_det2 = log_det * log_det;
843 t_dev(i, j) = t_log_stretch_h1(i, j) - t_kd(i, j) * (third * log_det);
844 double dev_norm2 = t_dev(i, j) * t_dev(i, j);
845
846 auto A = 2 * mu * std::exp(nohat_k * dev_norm2);
847 auto B = lambda * std::exp(hat_k * log_det2) * log_det;
848
849 FTensor::Tensor2_symmetric<double, 3> t_A_diff, t_B_diff;
850 t_A_diff(i, j) =
851 (A * 2 * nohat_k) * (t_dev(k, l) * t_diff_deviator(k, l, i, j));
852 t_B_diff(i, j) = (B * 2 * hat_k) * log_det * t_kd(i, j) +
853 lambda * std::exp(hat_k * log_det2) * t_kd(i, j);
855 t_dT(i, j, k, l) =
856 t_A_diff(i, j) * (t_dev(m, n) * t_diff_deviator(m, n, k, l))
857
858 +
859
860 A * t_diff_deviator(m, n, i, j) * t_diff_deviator(m, n, k, l)
861
862 +
863
864 t_B_diff(i, j) * t_kd(k, l);
865
866 t_dP(L, J) = -t_L(i, j, L) *
867 ((
868
869 t_dT(i, j, k, l)
870
871 +
872
873 (alphaU * ts_a) * (t_D(i, j, m, n) * t_diff(m, n, k, l)
874
875 )) *
876 t_L(k, l, J));
877
878 // Work of symmetric tensor on undefined tensor is equal to the work
879 // of the symmetric part of it
883 t_sym(i, j) = (t_approx_P_adjoint__dstretch(i, j) ||
884 t_approx_P_adjoint__dstretch(j, i));
885 t_sym(i, j) /= 2.0;
886 auto t_diff2_uP2 = EigenMatrix::getDiffDiffMat(
887 t_eigen_vals, t_eigen_vecs, EshelbianCore::f, EshelbianCore::d_f,
888 EshelbianCore::dd_f, t_sym, nbUniq[gg]);
889 t_dP(L, J) += t_L(i, j, L) *
890 ((t_diff2_uP2(i, j, k, l) + t_diff2_uP2(k, l, i, j)) *
891 t_L(k, l, J)) /
892 2.;
893 }
894
895 ++t_dP;
896 ++t_approx_P_adjoint__dstretch;
897 ++t_log_stretch_h1;
898 ++t_eigen_vals;
899 ++t_eigen_vecs;
900 }
901
902 return getFTensor2FromMat<size_symm, size_symm>(dP);
903 };
904
905 auto t_dP = get_dP();
906 for (int gg = 0; gg != nb_integration_pts; ++gg) {
907 double a = v * t_w;
908
909 int rr = 0;
910 for (; rr != row_nb_dofs / 6; ++rr) {
911 auto t_col_base_fun = col_data.getFTensor0N(gg, 0);
912 auto t_m = get_ftensor2(K, 6 * rr, 0);
913 for (int cc = 0; cc != col_nb_dofs / 6; ++cc) {
914 const double b = a * t_row_base_fun * t_col_base_fun;
915 t_m(L, J) -= b * t_dP(L, J);
916 ++t_m;
917 ++t_col_base_fun;
918 }
919 ++t_row_base_fun;
920 }
921
922 for (; rr != row_nb_base_functions; ++rr) {
923 ++t_row_base_fun;
924 }
925
926 ++t_w;
927 ++t_dP;
928 }
930}
931
933 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
934 boost::shared_ptr<double> total_energy_ptr)
935 : VolUserDataOperator(NOSPACE, OPSPACE), dataAtPts(data_ptr),
936 totalEnergyPtr(total_energy_ptr) {
937
938 if (!dataAtPts) {
940 "dataAtPts is not allocated. Please set it before "
941 "using this operator.");
942 }
943
944}
945
946MoFEMErrorCode HMHHencky::OpCalculateEnergy::doWork(int side, EntityType type,
947 EntData &data) {
949
954
955 int nb_integration_pts = getGaussPts().size2();
956 auto t_log_u =
957 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
958
959
960 auto &mat_d = dataAtPts->matD;
961 if (mat_d.size1() != size_symm || mat_d.size2() != size_symm) {
962 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
963 "wrong matD size, should be 6 by 6 but is %d by %d", size_symm,
964 size_symm);
965 }
966
967 auto t_D = getFTensor4DdgFromPtr<3, 3, 0>(mat_d.data().data());
968
969 dataAtPts->energyAtPts.resize(nb_integration_pts, false);
970 auto t_energy = getFTensor0FromVec(dataAtPts->energyAtPts);
971
972 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
973
974 t_energy = 0.5 * (t_log_u(i, j) * (t_D(i, j, k, l) * t_log_u(k, l)));
975
976 ++t_log_u;
977 ++t_energy;
978 }
979
980 if (totalEnergyPtr) {
981 auto t_w = getFTensor0IntegrationWeight();
982 auto t_energy = getFTensor0FromVec(dataAtPts->energyAtPts);
983 double loc_energy = 0;
984 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
985 loc_energy += t_energy * t_w;
986 ++t_w;
987 ++t_energy;
988 }
989 *totalEnergyPtr += getMeasure() * loc_energy;
990 }
991
993}
994
996 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
997 boost::shared_ptr<MatrixDouble> strain_ptr,
998 boost::shared_ptr<MatrixDouble> stress_ptr,
999 boost::shared_ptr<HMHHencky> hencky_ptr)
1000 : VolUserDataOperator(NOSPACE, OPSPACE), dataAtPts(data_ptr),
1001 strainPtr(strain_ptr), stressPtr(stress_ptr), henckyPtr(hencky_ptr) {}
1002
1004 EntityType type,
1005 EntData &data) {
1007
1014
1015 auto nb_integration_pts = stressPtr->size2();
1016#ifndef NDEBUG
1017 if (nb_integration_pts != getGaussPts().size2()) {
1018 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1019 "inconsistent number of integration points");
1020 }
1021#endif // NDEBUG
1022
1023 auto get_D = [&]() {
1025 for (auto &b : henckyPtr->blockData) {
1026
1027 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
1028 dataAtPts->mu = b.shearModulusG;
1029 dataAtPts->lambda = b.bulkModulusK - 2 * b.shearModulusG / 3;
1030 CHKERR henckyPtr->getMatDPtr(
1031 dataAtPts->getMatDPtr(), dataAtPts->getMatAxiatorDPtr(),
1032 dataAtPts->getMatDeviatorDPtr(), b.bulkModulusK, b.shearModulusG);
1033 CHKERR henckyPtr->getInvMatDPtr(dataAtPts->getMatInvDPtr(),
1034 b.bulkModulusK, b.shearModulusG);
1036 }
1037 }
1038
1039 const auto E = henckyPtr->E;
1040 const auto nu = henckyPtr->nu;
1041
1042 double bulk_modulus_K = E / (3 * (1 - 2 * nu));
1043 double shear_modulus_G = E / (2 * (1 + nu));
1044
1045 dataAtPts->mu = shear_modulus_G;
1046 dataAtPts->lambda = bulk_modulus_K - 2 * shear_modulus_G / 3;
1047
1048 CHKERR henckyPtr->getMatDPtr(
1049 dataAtPts->getMatDPtr(), dataAtPts->getMatAxiatorDPtr(),
1050 dataAtPts->getMatDeviatorDPtr(), bulk_modulus_K, shear_modulus_G);
1051 CHKERR henckyPtr->getInvMatDPtr(dataAtPts->getMatInvDPtr(), bulk_modulus_K,
1054 };
1055
1056 CHKERR get_D();
1057
1058 strainPtr->resize(size_symm, nb_integration_pts, false);
1059 auto t_strain = getFTensor2SymmetricFromMat<3>(*strainPtr);
1060 auto t_stress = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*stressPtr);
1061 auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(
1062 &*dataAtPts->matInvD.data().begin());
1063#ifndef NDEBUG
1064 auto t_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(
1065 &*dataAtPts->matD.data().begin());
1066#endif
1067
1068 const double lambda = dataAtPts->lambda;
1069 const double mu = dataAtPts->mu;
1071
1072 // note: add rotation, so we can extract rigid body motion, work then with
1073 // symmetric part.
1074 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
1075 t_strain(i, j) = t_inv_D(i, j, k, l) * t_stress(k, l);
1076
1077#ifndef NDEBUG
1078 FTensor::Tensor2_symmetric<double, 3> t_stress_symm_debug;
1079 t_stress_symm_debug(i, j) = (t_stress(i, j) || t_stress(j, i)) / 2;
1080 FTensor::Tensor2_symmetric<double, 3> t_stress_symm_debug_diff;
1081 t_stress_symm_debug_diff(i, j) =
1082 t_D(i, j, k, l) * t_strain(k, l) - t_stress_symm_debug(i, j);
1083 double nrm =
1084 t_stress_symm_debug_diff(i, j) * t_stress_symm_debug_diff(i, j);
1085 double nrm0 = t_stress_symm_debug(i, j) * t_stress_symm_debug(i, j) +
1086 std::numeric_limits<double>::epsilon();
1087 constexpr double eps = 1e-10;
1088 if (std::fabs(std::sqrt(nrm / nrm0)) > eps) {
1089 MOFEM_LOG("SELF", Sev::error)
1090 << "Stress symmetry check failed: " << std::endl
1091 << t_stress_symm_debug_diff << std::endl
1092 << t_stress;
1094 "Norm is too big: " + std::to_string(nrm / nrm0));
1095 }
1096 ++t_D;
1097#endif
1098
1099 ++t_strain;
1100 ++t_stress;
1101 ++t_inv_D;
1102 }
1103
1105}
1106
1107} // namespace EshelbianPlasticity
constexpr double third
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
#define FTENSOR_INDEX(DIM, I)
constexpr double a
static PetscErrorCode ierr
static const double eps
constexpr int SPACE_DIM
Kronecker Delta class symmetric.
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ NOSPACE
Definition definitions.h:83
#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.
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
double bulk_modulus_K
double shear_modulus_G
constexpr auto t_kd
auto get_D
#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 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.
auto diff_deviator(FTensor::Ddg< double, 3, 3 > &&t_diff_stress)
static constexpr auto size_symm
double young_modulus
Young modulus.
Definition plastic.cpp:125
double poisson_ratio
Poisson ratio.
Definition plastic.cpp:126
double scale
Definition plastic.cpp:123
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
Calculate energy density for Hencky material model.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpCalculateEnergy(boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< double > total_energy_ptr)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpCalculateStretchFromStress(boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< MatrixDouble > strain_ptr, boost::shared_ptr< MatrixDouble > stress_ptr, boost::shared_ptr< HMHHencky > hencky_ptr)
MoFEMErrorCode evaluateRhs(EntData &data)
Definition HMHHencky.cpp:61
boost::shared_ptr< HMHHencky > henckyPtr
Definition HMHHencky.cpp:66
MoFEMErrorCode evaluateLhs(EntData &data)
Definition HMHHencky.cpp:62
OpHenckyJacobian(boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< HMHHencky > hencky_ptr)
Definition HMHHencky.cpp:20
boost::shared_ptr< DataAtIntegrationPts > dataAtGaussPts
Definition HMHHencky.cpp:65
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition HMHHencky.cpp:29
std::map< std::string, boost::shared_ptr< ScalingMethod > > scalingMethodsMap
boost::shared_ptr< ExternalStrainVec > externalStrainVecPtr
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 integrateHencky(EntData &row_data, EntData &col_data)
MoFEMErrorCode integratePolyconvexHencky(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integratePolyconvexHencky(EntData &data)
OpSpatialPhysical(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const double alpha_u)
MoFEMErrorCode integrateHencky(EntData &data)
virtual VolUserDataOperator * returnOpSpatialPhysical(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const double alpha_u)
Definition HMHHencky.cpp:95
VolUserDataOperator * returnOpCalculateStretchFromStress(boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< PhysicalEquations > physics_ptr)
MoFEMErrorCode extractBlockData(Sev sev)
std::vector< BlockData > blockData
MoFEMErrorCode extractBlockData(std::vector< const CubitMeshSets * > meshset_vec_ptr, Sev sev)
HMHHencky(MoFEM::Interface &m_field, const double E, const double nu)
Definition HMHHencky.cpp:14
VolUserDataOperator * returnOpCalculateVarStretchFromStress(boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< PhysicalEquations > physics_ptr)
VolUserDataOperator * returnOpCalculateEnergy(boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< double > total_energy_ptr)
MoFEMErrorCode getMatDPtr(boost::shared_ptr< MatrixDouble > mat_D_ptr, boost::shared_ptr< MatrixDouble > mat_axiator_D_ptr, boost::shared_ptr< MatrixDouble > mat_deviator_D_ptr, double bulk_modulus_K, double shear_modulus_G)
MoFEMErrorCode getOptions(boost::shared_ptr< DataAtIntegrationPts > data_ptr)
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)
Definition HMHHencky.cpp:70
VolUserDataOperator * returnOpSpatialPhysical_du_du(std::string row_field, std::string col_field, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const double alpha)
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)
MoFEMErrorCode recordTape(const int tag, DTensor2Ptr *t_h)
Definition HMHHencky.cpp:17
MoFEMErrorCode getInvMatDPtr(boost::shared_ptr< MatrixDouble > mat_inv_D_ptr, double bulk_modulus_K, double shear_modulus_G)
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::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.