v0.14.0
Loading...
Searching...
No Matches
PlasticOperators.hpp
Go to the documentation of this file.
1/* This file is part of MoFEM.
2 * MoFEM is free software: you can redistribute it and/or modify it under
3 * the terms of the GNU Lesser General Public License as published by the
4 * Free Software Foundation, either version 3 of the License, or (at your
5 * option) any later version.
6 *
7 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
8 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
9 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
10 * License for more details.
11 *
12 * You should have received a copy of the GNU Lesser General Public
13 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
14
15/** \file OpPlasticOperators.hpp
16
17\f[
18\left\{
19\begin{array}{ll}
20\frac{\partial \sigma_{ij}}{\partial x_j} - b_i = 0 & \forall x \in \Omega \\
21\varepsilon_{ij} = \frac{1}{2}\left( \frac{\partial u_i}{\partial x_j} +
22\frac{\partial u_j}{\partial x_i} \right)\\
23\sigma_{ij} = D_{ijkl}\left(\varepsilon_{kl}-\varepsilon^p_{kl}\right) \\
24\dot{\varepsilon}^p_{kl} - \dot{\tau} \left( \left. \frac{\partial f}{\partial
25\sigma_{kl}} \right|_{(\sigma,\tau) } \right) = 0 \\
26f(\sigma, \tau) \leq 0,\; \dot{\tau} \geq 0,\;\dot{\tau}f(\sigma, \tau)=0\\
27u_i = \overline{u}_i & \forall x \in \partial\Omega_u \\
28\sigma_{ij}n_j = \overline{t}_i & \forall x \in \partial\Omega_\sigma \\
29\Omega_u \cup \Omega_\sigma = \Omega \\
30\Omega_u \cap \Omega_\sigma = \emptyset
31\end{array}
32\right.
33\f]
34
35\f[
36\left\{
37\begin{array}{ll}
38\left(\frac{\partial \delta u_i}{\partial x_j},\sigma_{ij}\right)_\Omega-(\delta
39u_i,b_i)_\Omega -(\delta u_i,\overline{t}_i)_{\partial\Omega_\sigma}=0 & \forall
40\delta u_i \in H^1(\Omega)\\ \left(\delta\varepsilon^p_{kl} ,D_{ijkl}\left(
41\dot{\varepsilon}^p_{kl} - \dot{\tau} A_{kl} \right)\right) = 0
42& \forall \delta\varepsilon^p_{ij} \in L^2(\Omega) \cap \mathcal{S} \\
43\left(\delta\tau,c_n\dot{\tau} - \frac{1}{2}\left\{c_n \dot{\tau} +
44(f(\pmb\sigma,\tau) - \sigma_y) +
45\| c_n \dot{\tau} + (f(\pmb\sigma,\tau) - \sigma_y) \|\right\}\right) = 0 &
46\forall \delta\tau \in L^2(\Omega) \end{array} \right. \f]
47
48*/
49
50namespace OpPlasticTools {
51
52//! [Common data]
54 boost::shared_ptr<VectorDouble> plasticSurfacePtr;
55 boost::shared_ptr<MatrixDouble> plasticFlowPtr;
56 boost::shared_ptr<VectorDouble> plasticTauPtr;
57 boost::shared_ptr<VectorDouble> plasticTauDotPtr;
58 boost::shared_ptr<MatrixDouble> plasticStrainPtr;
59 boost::shared_ptr<MatrixDouble> plasticStrainDotPtr;
60
61 boost::shared_ptr<MatrixDouble> plasticGradTauPtr;
62 boost::shared_ptr<MatrixDouble> plasticGradStrainPtr;
63
64 boost::shared_ptr<VectorDouble> plasticTauJumpPtr;
65 boost::shared_ptr<MatrixDouble> plasticStrainJumpPtr;
66 boost::shared_ptr<MatrixDouble> guidingVelocityPtr;
67
68 Ddg<double, 3, 3> diffDeviator;
69
70 MatrixDouble dualBaseMat;
72
73 // data for skeleton computation
74 map<int, MatrixDouble> plasticStrainSideMap;
75 map<int, VectorDouble> plasticTauSideMap;
76 map<int, MatrixDouble> velocityVecSideMap;
77
78 MoFEMErrorCode calculateDiffDeviator();
79};
80
83 const std::string field_name,
84 boost::shared_ptr<CommonData> common_data_ptr);
85 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
86
87private:
88 boost::shared_ptr<CommonData> commonDataPtr;
89};
90
92 OpPlasticStress(const std::string field_name,
93 boost::shared_ptr<CommonData> common_data_ptr,
94 boost::shared_ptr<MatrixDouble> mat_d_ptr);
95 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
96
97private:
98 boost::shared_ptr<CommonData> commonDataPtr;
99 boost::shared_ptr<MatrixDouble> matDPtr;
100};
101
103 OpCalculatePlasticFlowRhs(const std::string field_name,
104 boost::shared_ptr<CommonData> common_data_ptr);
105 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
106
107private:
108 boost::shared_ptr<CommonData> commonDataPtr;
109};
110
112 OpCalculateConstraintRhs(const std::string field_name,
113 boost::shared_ptr<CommonData> common_data_ptr);
114 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
115
116private:
117 boost::shared_ptr<CommonData> commonDataPtr;
118};
119
120template <bool LOGSTRAIN>
123 const std::string row_field_name, const std::string col_field_name,
124 boost::shared_ptr<CommonData> common_data_ptr,
125 boost::shared_ptr<MatrixDouble> mat_D_ptr);
126 MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
127
128private:
129 boost::shared_ptr<CommonData> commonDataPtr;
130 boost::shared_ptr<MatrixDouble> matDPtr;
131};
132
133template <bool LOGSTRAIN>
135 OpCalculatePlasticFlowLhs_dU(const std::string row_field_name,
136 const std::string col_field_name,
137 boost::shared_ptr<CommonData> common_data_ptr);
138 MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
139
140private:
141 boost::shared_ptr<CommonData> commonDataPtr;
142};
143
145 OpCalculatePlasticFlowLhs_dEP(const std::string row_field_name,
146 const std::string col_field_name,
147 boost::shared_ptr<CommonData> common_data_ptr);
148 MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
149
150private:
151 boost::shared_ptr<CommonData> commonDataPtr;
152};
153
155 OpCalculatePlasticFlowLhs_dTAU(const std::string row_field_name,
156 const std::string col_field_name,
157 boost::shared_ptr<CommonData> common_data_ptr);
158 MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
159
160private:
161 boost::shared_ptr<CommonData> commonDataPtr;
162};
163
165 OpCalculatePlasticInternalForceLhs_dTAU(const std::string row_field_name,
166 const std::string col_field_name,
167 boost::shared_ptr<CommonData> common_data_ptr);
168 MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
169
170private:
171 boost::shared_ptr<CommonData> commonDataPtr;
172};
173
174template <bool LOGSTRAIN>
176 OpCalculateConstraintLhs_dU(const std::string row_field_name,
177 const std::string col_field_name,
178 boost::shared_ptr<CommonData> common_data_ptr);
179 MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
180
181private:
182 boost::shared_ptr<CommonData> commonDataPtr;
183};
184
186 OpCalculateConstraintLhs_dEP(const std::string row_field_name,
187 const std::string col_field_name,
188 boost::shared_ptr<CommonData> common_data_ptr);
189 MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
190
191private:
192 boost::shared_ptr<CommonData> commonDataPtr;
193};
194
196 OpCalculateConstraintLhs_dTAU(const std::string row_field_name,
197 const std::string col_field_name,
198 boost::shared_ptr<CommonData> common_data_ptr);
199 MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
200
201private:
202 boost::shared_ptr<CommonData> commonDataPtr;
203};
204
206 OpPostProcPlastic(const std::string field_name,
207 moab::Interface &post_proc_mesh,
208 std::vector<EntityHandle> &map_gauss_pts,
209 boost::shared_ptr<CommonData> common_data_ptr);
210 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
211
212private:
213 moab::Interface &postProcMesh;
214 std::vector<EntityHandle> &mapGaussPts;
215 boost::shared_ptr<CommonData> commonDataPtr;
216};
217//! [Operators definitions]
218
219//! [Lambda functions]
220inline auto diff_tensor() {
221 Ddg<double, 3, 3> t_diff;
222 constexpr auto t_kd = Kronecker_Delta_symmetric<int>();
223 t_diff(i, j, k, l) = (t_kd(i, k) ^ t_kd(j, l)) / 4.;
224 return t_diff;
225};
226
227inline auto symm_L_tensor() {
228 Dg<double, 3, 6> L;
229 Index<'L', 6> Li;
230 // dg(T 000, T 001, T 002, T 010, T 011, T 012, T 020, T 021, T 022,
231 // T 110, T 111, T 112, T 120, T 121, T 122, T 220, T 221, T 222);
232
233 // Tensor2_symmetric_constructor(T data[(3 * 4) / 2], T d00, T d01, T d02, T
234 // d11,
235 // T d12, T d22);
236 // Tensor2<double, 3, 3> test_tensor{1., 2., 3., 2., 5., 6., 3., 6., 9.};
237 // Tensor1<double, 6> test_vec;
238
239 L(i, j, Li) = 0;
240 L(0, 0, 0) = 1;
241 L(0, 1, 1) = 1;
242 L(0, 2, 2) = 1;
243 L(1, 0, 1) = 1;
244 L(1, 1, 3) = 1;
245 L(1, 2, 4) = 1;
246 L(2, 0, 2) = 1;
247 L(2, 1, 4) = 1;
248 L(2, 2, 5) = 1;
249
250 // test_vec(Li) = test_tensor(i, j) * L(i, j, Li);
251
252 // (00, 01, 02, 11, 12, 22)
253 return L;
254}
255
256inline auto diff_symmetrize() { return (*cache).Is; };
257
258template <typename T1>
259inline auto deviator(Tensor2_symmetric<T1, 3> &t_stress) {
260 Tensor2_symmetric<double, 3> t_dev;
261 t_dev(I, J) = 0;
262 constexpr double third = boost::math::constants::third<double>();
263 const double trace = t_stress(i, i);
264
265 t_dev(i, j) = t_stress(i, j);
266 t_dev(0, 0) -= trace * third;
267 t_dev(1, 1) -= trace * third;
268 t_dev(2, 2) -= trace * third;
269 return t_dev;
270};
271
272inline auto diff_deviator(Ddg<double, 3, 3> &&t_diff_stress) {
273 Ddg<double, 3, 3> t_diff_deviator;
274 t_diff_deviator(I, J, K, L) = 0;
275 for (int ii = 0; ii != 3; ++ii)
276 for (int jj = ii; jj != 3; ++jj)
277 for (int kk = 0; kk != 3; ++kk)
278 for (int ll = kk; ll != 3; ++ll)
279 t_diff_deviator(ii, jj, kk, ll) = t_diff_stress(ii, jj, kk, ll);
280
281 constexpr double third = boost::math::constants::third<double>();
282
283 t_diff_deviator(0, 0, 0, 0) = 1. - third;
284 t_diff_deviator(0, 0, 1, 1) = -third;
285 t_diff_deviator(0, 0, 2, 2) = -third;
286
287 t_diff_deviator(1, 1, 0, 0) = -third;
288 t_diff_deviator(1, 1, 1, 1) = 1. - third;
289 t_diff_deviator(1, 1, 2, 2) = -third;
290
291 t_diff_deviator(2, 2, 0, 0) = -third;
292 t_diff_deviator(2, 2, 1, 1) = -third;
293 t_diff_deviator(2, 2, 2, 2) = 1. - third;
294
295 return t_diff_deviator;
296};
297
298template <typename T>
299inline auto get_effective_stress(Tensor2_symmetric<T, 3> &t_stress,
300 Tensor2_symmetric<T, 3> &t_plastic_strain) {
301
302 Tensor2_symmetric<double, 3> stress_tmp;
303
304 stress_tmp(i, j) =
305 t_stress(i, j) - 1.5 * (*cache).C1_k * t_plastic_strain(i, j);
306
307 return stress_tmp;
308};
309
310constexpr double sqrt2by3 = 0.816496580927726;
311inline auto hardening(long double tau) {
312 return (*cache).H * tau + (*cache).Q_inf * (1. - exp(-(*cache).b_iso * tau)) +
313 (*cache).sigmaY;
314}
315
316inline auto hardening_dtau(long double tau) {
317 return (*cache).H +
318 (*cache).Q_inf * (*cache).b_iso * exp(-(*cache).b_iso * tau);
319}
320
321/**
322 *
323
324\f[
325\begin{split}
326f&=\sqrt{s_{ij}s_{ij}}\\
327A_{ij}&=\frac{\partial f}{\partial \sigma_{ij}}=
328\frac{1}{f} s_{kl} \frac{\partial s_{kl}}{\partial \sigma_{ij}}\\
329\frac{\partial A_{ij}}{\partial \sigma_{kl}}&= \frac{\partial^2 f}{\partial
330\sigma_{ij}\partial\sigma_{mn}}= \frac{1}{f} \left( \frac{\partial
331s_{kl}}{\partial \sigma_{mn}}\frac{\partial s_{kl}}{\partial \sigma_{ij}}
332-A_{mn}A_{ij}
333\right)\\
334\frac{\partial f}{\partial \varepsilon_{ij}}&=A_{mn}D_{mnij}
335\\
336\frac{\partial A_{ij}}{\partial \varepsilon_{kl}}&=
337\frac{\partial A_{ij}}{\partial \sigma_{mn}} \frac{\partial
338\sigma_{mn}}{\partial \varepsilon_{kl}}= \frac{\partial A_{ij}}{\partial
339\sigma_{mn}} D_{mnkl}
340\end{split}
341\f]
342
343 */
344inline auto plastic_surface(Tensor2_symmetric<double, 3> &&t_stress_deviator) {
345 return std::sqrt(1.5 * t_stress_deviator(I, J) * t_stress_deviator(I, J)) +
346 std::numeric_limits<double>::epsilon();
347};
348
349inline auto plastic_flow(long double f,
350 Tensor2_symmetric<double, 3> &&t_dev_stress,
351 Ddg<double, 3, 3> &t_diff_deviator) {
352 Tensor2_symmetric<double, 3> t_diff_f;
353 t_diff_f(k, l) =
354 ((1.5) * (t_dev_stress(I, J)) * t_diff_deviator(I, J, k, l)) / f;
355
356 return t_diff_f;
357};
358
359// FIXME: get rid of second projection for efficiency
360// possibly we can also remove D : from the flow rule for the same reason
361
362template <typename T>
363inline auto diff_plastic_flow_dstress(long double f,
364 Tensor2_symmetric<T, 3> &t_flow,
365 Ddg<double, 3, 3> &t_diff_deviator) {
366 Ddg<double, 3, 3> t_diff_flow;
367 t_diff_flow(i, j, k, l) =
368 ((1.5) * (t_diff_deviator(M, N, i, j) * t_diff_deviator(M, N, k, l) -
369 2. / 3. * t_flow(i, j) * t_flow(k, l))) /
370 f;
371
372 return t_diff_flow;
373};
374
375template <typename T>
376inline auto
378 Ddg<double, 3, 3> &&t_diff_plastic_flow_dstress) {
379 Ddg<double, 3, 3> t_diff_flow;
380 t_diff_flow(i, j, k, l) =
381 t_diff_plastic_flow_dstress(i, j, m, n) * t_D(m, n, k, l);
382
383 return t_diff_flow;
384};
385
386template <typename T>
388 Ddg<T, 3, 3> &t_D, Ddg<double, 3, 3> &&t_diff_plastic_flow_dstress) {
389 Ddg<double, 3, 3> t_diff_flow_kin_hard;
390 t_diff_flow_kin_hard(i, j, k, l) =
391 t_diff_plastic_flow_dstress(i, j, m, n) * t_D(m, n, k, l) +
392 1.5 * (*cache).C1_k * t_diff_plastic_flow_dstress(i, j, k, l);
393
394 return t_diff_flow_kin_hard;
395};
396
397inline double constrain_abs(long double x) { return std::abs(x); };
398
399inline double constraint_sign(long double x) {
400 if (x > 0)
401 return 1;
402 else if (x < 0)
403 return -1;
404 else
405 return 0;
406};
407
408inline double w(long double tau_dot, long double f, long double hardening) {
409 auto &cn = (*cache).cn_pl;
410 auto &sigmaY = (*cache).sigmaY;
411 return (f - hardening) / sigmaY + cn * tau_dot;
412};
413
414inline double constraint(long double tau_dot, long double f,
415 long double hardening) {
416 auto &cn = (*cache).cn_pl;
417 auto &visH = (*cache).visH;
418 auto &sigmaY = (*cache).sigmaY;
419 return visH * tau_dot +
420 (sigmaY / 2) * ((cn * tau_dot - (f - hardening) / sigmaY) -
421 constrain_abs(w(tau_dot, f, hardening)));
422};
423
424inline double diff_constrain_dtau_dot(long double tau_dot, long double f,
425 long double hardening) {
426 auto &cn = (*cache).cn_pl;
427 auto &visH = (*cache).visH;
428 auto &sigmaY = (*cache).sigmaY;
429 return visH +
430 (sigmaY / 2) * (cn - cn * constraint_sign(w(tau_dot, f, hardening)));
431};
432
433inline auto diff_constrain_df(long double tau_dot, long double f,
434 long double hardening) {
435 return (-1 - constraint_sign(w(tau_dot, f, hardening))) / 2;
436};
437
438inline auto diff_constrain_dsigma_y(long double tau_dot, long double f,
439 long double hardening) {
440 return (1 + constraint_sign(w(tau_dot, f, hardening))) / 2;
441}
442
443template <typename T>
445 Tensor2_symmetric<T, 3> &t_plastic_flow) {
446 Tensor2_symmetric<double, 3> t_diff_constrain_dstress;
447 t_diff_constrain_dstress(i, j) = diff_constrain_df * t_plastic_flow(i, j);
448 return t_diff_constrain_dstress;
449};
450
451template <typename T1, typename T2>
452inline auto
453diff_constrain_dstrain(Ddg<T1, 3, 3> &t_D,
454 Tensor2_symmetric<T2, 3> &&t_diff_constrain_dstress) {
455 Tensor2_symmetric<double, 3> t_diff_constrain_dstrain;
456 t_diff_constrain_dstrain(k, l) =
457 t_diff_constrain_dstress(i, j) * t_D(i, j, k, l);
458 return t_diff_constrain_dstrain;
459};
460
461template <typename T1, typename T2>
463 Ddg<T1, 3, 3> &t_D, Tensor2_symmetric<T2, 3> &&t_diff_constrain_dstress) {
464 Tensor2_symmetric<double, 3> t_diff_constrain_kin_hard_dstrain;
465 t_diff_constrain_kin_hard_dstrain(k, l) =
466 t_diff_constrain_dstress(i, j) * t_D(i, j, k, l) +
467 1.5 * (*cache).C1_k * t_diff_constrain_dstress(k, l);
468 return t_diff_constrain_kin_hard_dstrain;
469};
470
472
474 boost::shared_ptr<CommonData> common_data_ptr)
476 commonDataPtr(common_data_ptr) {
477 // Operator is only executed for vertices
478 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
479 }
480
481 //! [Calculate stress]
482 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
484 // const size_t nb_gauss_pts = commonDataPtr->mStrainPtr->size2();
485 auto t_tau = getFTensor0FromVec(*(commonDataPtr->plasticTauPtr));
486 auto t_tau_dot = getFTensor0FromVec(*(commonDataPtr->plasticTauDotPtr));
487 auto &data_vec = commonDataPtr->stateArrayPlast;
488 for (auto &f : *(commonDataPtr->plasticSurfacePtr)) {
489
490 if (w(t_tau_dot, f, hardening(t_tau)) > 0)
491 data_vec[0] += 1;
492 data_vec[1] += 1;
493
494 ++t_tau_dot;
495 ++t_tau;
496 }
497
499 }
500
501private:
502 boost::shared_ptr<CommonData> commonDataPtr;
503};
504
506
509
510 /**
511 * @brief Return interface to this class when one ask for for tetrahedron,
512 * otherisw return interface class for generic class.
513 *
514 * @param iface interface class
515 * @return MoFEMErrorCode
516 */
517 MoFEMErrorCode query_interface(boost::typeindex::type_index type_index,
518 UnknownInterface **iface) const {
520 *iface = const_cast<BubbleTauPolynomialBase *>(this);
522 }
523
524 /**
525 * @brief Calculate base functions at intergeneration points
526 *
527 * @param pts
528 * @param ctx_ptr
529 * @return MoFEMErrorCode
530 */
531 MoFEMErrorCode getValue(MatrixDouble &pts,
532 boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
534
535 cTx = ctx_ptr->getInterface<EntPolynomialBaseCtx>();
536
537 int nb_gauss_pts = pts.size2();
538 if (!nb_gauss_pts) {
540 }
541
542 if (pts.size1() < 3) {
543 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
544 "Wrong dimension of pts, should be at least 3 rows with "
545 "coordinates");
546 }
547
548 switch (cTx->sPace) {
549 case L2:
550 // CHKERR getValueL2ForTauBubble(pts);
551 break;
552 default:
553 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
554 }
555
557 }
558
559private:
560 EntPolynomialBaseCtx *cTx;
561
562 // MatrixDouble shapeFun;
563
564 // MoFEMErrorCode getValueL2ForTauBubble(MatrixDouble &pts) {
565 // MoFEMFunctionBegin;
566
567 // const FieldApproximationBase base = cTx->bAse;
568 // // This should be used only in case USER_BASE is selected
569 // if (cTx->bAse != USER_BASE) {
570 // SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
571 // "Wrong base, should be USER_BASE");
572 // }
573
574 // EntitiesFieldData &data = cTx->dAta;
575 // int nb_gauss_pts = pts.size2();
576
577 // data.dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 4, false);
578 // data.dataOnEntities[MBTET][0].getN(base).resize(
579 // nb_gauss_pts,
580 // NBVOLUMETET_L2(data.dataOnEntities[MBTET][0].getDataOrder()), false);
581 // data.dataOnEntities[MBTET][0].getDiffN(base).resize(
582 // nb_gauss_pts,
583 // 3 * NBVOLUMETET_L2(data.dataOnEntities[MBTET][0].getDataOrder()),
584 // false);
585
586 // CHKERR L2_Ainsworth_ShapeFunctions_MBTET(
587 // data.dataOnEntities[MBTET][0].getDataOrder(),
588 // &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
589 // &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
590 // &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
591 // &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
592 // nb_gauss_pts, cTx->basePolynomialsType0);
593
594 // auto tetN = data.dataOnEntities[MBTET][0].getN(base);
595 // data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 2);
596 // for (int gg = 0; gg != nb_gauss_pts; ++gg) {
597 // data.dataOnEntities[MBTET][0].getN(base)(gg, 0) = 1;
598 // auto prod_n0123 = tetN(gg, 0) * tetN(gg, 1) * tetN(gg, 2) * tetN(gg, 3);
599 // data.dataOnEntities[MBTET][0].getN(base)(gg, 1) = prod_n0123;
600 // }
601
602 // MoFEMFunctionReturn(0);
603 // }
604};
605
606//! [Postprocessing]
607}; // namespace OpPlasticTools
static Index< 'M', 3 > M
static Index< 'L', 3 > L
static Index< 'J', 3 > J
static Index< 'K', 3 > K
constexpr double third
Definition: single.cpp:4
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ L2
field with C-1 continuity
Definition: definitions.h:88
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ 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()
Definition: definitions.h:416
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
double trace(FTensor::Tensor2< T, 2, 2 > &t_stress)
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
constexpr auto t_kd
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
auto plastic_flow(long double f, Tensor2_symmetric< double, 3 > &&t_dev_stress, Ddg< double, 3, 3 > &t_diff_deviator)
double constraint(long double tau_dot, long double f, long double hardening)
auto deviator(Tensor2_symmetric< T1, 3 > &t_stress)
auto diff_plastic_flow_dstrain(Ddg< T, 3, 3 > &t_D, Ddg< double, 3, 3 > &&t_diff_plastic_flow_dstress)
auto diff_plastic_flow_kin_hard_dstrain(Ddg< T, 3, 3 > &t_D, Ddg< double, 3, 3 > &&t_diff_plastic_flow_dstress)
auto diff_constrain_dstrain(Ddg< T1, 3, 3 > &t_D, Tensor2_symmetric< T2, 3 > &&t_diff_constrain_dstress)
double w(long double tau_dot, long double f, long double hardening)
auto plastic_surface(Tensor2_symmetric< double, 3 > &&t_stress_deviator)
double constrain_abs(long double x)
double diff_constrain_dtau_dot(long double tau_dot, long double f, long double hardening)
auto diff_constrain_dsigma_y(long double tau_dot, long double f, long double hardening)
auto diff_constrain_kin_hard_dstrain(Ddg< T1, 3, 3 > &t_D, Tensor2_symmetric< T2, 3 > &&t_diff_constrain_dstress)
auto get_effective_stress(Tensor2_symmetric< T, 3 > &t_stress, Tensor2_symmetric< T, 3 > &t_plastic_strain)
double constraint_sign(long double x)
auto hardening_dtau(long double tau)
auto diff_deviator(Ddg< double, 3, 3 > &&t_diff_stress)
auto diff_plastic_flow_dstress(long double f, Tensor2_symmetric< T, 3 > &t_flow, Ddg< double, 3, 3 > &t_diff_deviator)
auto diff_constrain_df(long double tau_dot, long double f, long double hardening)
auto diff_tensor()
[Operators definitions]
auto hardening(long double tau)
auto diff_constrain_dstress(long double &&diff_constrain_df, Tensor2_symmetric< T, 3 > &t_plastic_flow)
constexpr double sqrt2by3
constexpr IntegrationType I
double visH
Viscous hardening.
Definition: plastic.cpp:172
double sigmaY
Yield stress.
Definition: plastic.cpp:170
constexpr auto field_name
const int N
Definition: speed_test.cpp:3
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Data on single entity (This is passed as argument to DataOperator::doWork)
@ OPROW
operator doWork function is executed on FE rows
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Calculate base functions at intergeneration points.
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Return interface to this class when one ask for for tetrahedron, otherisw return interface class for ...
boost::shared_ptr< VectorDouble > plasticTauDotPtr
MoFEMErrorCode calculateDiffDeviator()
boost::shared_ptr< MatrixDouble > plasticGradStrainPtr
boost::shared_ptr< MatrixDouble > guidingVelocityPtr
map< int, VectorDouble > plasticTauSideMap
Ddg< double, 3, 3 > diffDeviator
boost::shared_ptr< VectorDouble > plasticTauPtr
map< int, MatrixDouble > velocityVecSideMap
map< int, MatrixDouble > plasticStrainSideMap
boost::shared_ptr< VectorDouble > plasticSurfacePtr
boost::shared_ptr< MatrixDouble > plasticFlowPtr
boost::shared_ptr< MatrixDouble > plasticStrainDotPtr
boost::shared_ptr< MatrixDouble > plasticStrainPtr
boost::shared_ptr< MatrixDouble > plasticGradTauPtr
boost::shared_ptr< VectorDouble > plasticTauJumpPtr
boost::shared_ptr< MatrixDouble > plasticStrainJumpPtr
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
boost::shared_ptr< CommonData > commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
boost::shared_ptr< CommonData > commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Calculate stress]
OpGetGaussPtsPlasticState(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
boost::shared_ptr< CommonData > commonDataPtr
boost::shared_ptr< MatrixDouble > matDPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Calculate stress]
boost::shared_ptr< CommonData > commonDataPtr
std::vector< EntityHandle > & mapGaussPts
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Postprocessing]
moab::Interface & postProcMesh
boost::shared_ptr< CommonData > commonDataPtr