v0.10.0
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 PlasticOps.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 \\
26 f(\sigma, \tau) \leq 0,\; \dot{\tau} \geq 0,\;\dot{\tau}f(\sigma, \tau)=0\\
27 u_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
39 u_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 
50 namespace OpPlasticTools {
51 
52 //! [Common data]
53 struct CommonData : public OpElasticTools::CommonData {
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  Ddg<double, 3, 2> diffDeviator;
62 
64 };
65 //! [Common data]
66 
73 
80 
81 struct OpCalculatePlasticSurface : public DomainEleOp {
82  OpCalculatePlasticSurface(const std::string field_name,
83  boost::shared_ptr<CommonData> common_data_ptr);
84  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
85 
86 private:
87  boost::shared_ptr<CommonData> commonDataPtr;
88 };
89 
90 struct OpPlasticStress : public DomainEleOp {
91  OpPlasticStress(const std::string field_name,
92  boost::shared_ptr<CommonData> common_data_ptr);
93  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
94 
95 private:
96  boost::shared_ptr<CommonData> commonDataPtr;
97 };
98 
99 struct OpCalculatePlasticFlowRhs : public DomainEleOp {
100  OpCalculatePlasticFlowRhs(const std::string field_name,
101  boost::shared_ptr<CommonData> common_data_ptr);
102  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
103 
104 private:
105  boost::shared_ptr<CommonData> commonDataPtr;
106 };
107 
108 struct OpCalculateContrainsRhs : public DomainEleOp {
109  OpCalculateContrainsRhs(const std::string field_name,
110  boost::shared_ptr<CommonData> common_data_ptr);
111  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
112 
113 private:
114  boost::shared_ptr<CommonData> commonDataPtr;
115 };
116 
117 struct OpCalculatePlasticInternalForceLhs_dEP : public DomainEleOp {
119  const std::string row_field_name, const std::string col_field_name,
120  boost::shared_ptr<CommonData> common_data_ptr);
121  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
122  EntityType col_type, EntData &row_data,
123  EntData &col_data);
124 
125 private:
126  boost::shared_ptr<CommonData> commonDataPtr;
128 };
129 
130 struct OpCalculatePlasticFlowLhs_dU : public DomainEleOp {
131  OpCalculatePlasticFlowLhs_dU(const std::string row_field_name,
132  const std::string col_field_name,
133  boost::shared_ptr<CommonData> common_data_ptr);
134  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
135  EntityType col_type, EntData &row_data,
136  EntData &col_data);
137 
138 private:
139  boost::shared_ptr<CommonData> commonDataPtr;
141 };
142 
143 struct OpCalculatePlasticFlowLhs_dEP : public DomainEleOp {
144  OpCalculatePlasticFlowLhs_dEP(const std::string row_field_name,
145  const std::string col_field_name,
146  boost::shared_ptr<CommonData> common_data_ptr);
147  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
148  EntityType col_type, EntData &row_data,
149  EntData &col_data);
150 
151 private:
152  boost::shared_ptr<CommonData> commonDataPtr;
154 };
155 
156 struct OpCalculatePlasticFlowLhs_dTAU : public DomainEleOp {
157  OpCalculatePlasticFlowLhs_dTAU(const std::string row_field_name,
158  const std::string col_field_name,
159  boost::shared_ptr<CommonData> common_data_ptr);
160  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
161  EntityType col_type, EntData &row_data,
162  EntData &col_data);
163 
164 private:
165  boost::shared_ptr<CommonData> commonDataPtr;
167 };
168 
169 struct OpCalculateContrainsLhs_dU : public DomainEleOp {
170  OpCalculateContrainsLhs_dU(const std::string row_field_name,
171  const std::string col_field_name,
172  boost::shared_ptr<CommonData> common_data_ptr);
173  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
174  EntityType col_type, EntData &row_data,
175  EntData &col_data);
176 
177 private:
178  boost::shared_ptr<CommonData> commonDataPtr;
180 };
181 
182 struct OpCalculateContrainsLhs_dEP : public DomainEleOp {
183  OpCalculateContrainsLhs_dEP(const std::string row_field_name,
184  const std::string col_field_name,
185  boost::shared_ptr<CommonData> common_data_ptr);
186  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
187  EntityType col_type, EntData &row_data,
188  EntData &col_data);
189 
190 private:
191  boost::shared_ptr<CommonData> commonDataPtr;
193 };
194 
195 struct OpCalculateContrainsLhs_dTAU : public DomainEleOp {
196  OpCalculateContrainsLhs_dTAU(const std::string row_field_name,
197  const std::string col_field_name,
198  boost::shared_ptr<CommonData> common_data_ptr);
199  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
200  EntityType col_type, EntData &row_data,
201  EntData &col_data);
202 
203 private:
204  boost::shared_ptr<CommonData> commonDataPtr;
206 };
207 
208 struct OpPostProcPlastic : public DomainEleOp {
209  OpPostProcPlastic(const std::string field_name,
210  moab::Interface &post_proc_mesh,
211  std::vector<EntityHandle> &map_gauss_pts,
212  boost::shared_ptr<CommonData> common_data_ptr);
213  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
214 
215 private:
217  std::vector<EntityHandle> &mapGaussPts;
218  boost::shared_ptr<CommonData> commonDataPtr;
219 };
220 //! [Operators definitions]
221 
222 //! [Lambda functions]
223 inline auto diff_tensor() {
224  Ddg<double, 2, 2> t_diff;
225  t_diff(i, j, k, l) = 0;
226  for (size_t ii = 0; ii != 2; ++ii)
227  for (size_t jj = ii; jj != 2; ++jj)
228  t_diff(ii, jj, ii, jj) = 1;
229  return t_diff;
230 };
231 
232 inline auto diff_symmetrize() {
233  Tensor4<double, 2, 2, 2, 2> t_diff;
234  t_diff(i, j, k, l) = 0;
235  t_diff(0, 0, 0, 0) = 1;
236  t_diff(1, 1, 1, 1) = 1;
237  t_diff(1, 0, 1, 0) = 0.25;
238  t_diff(0, 1, 1, 0) = 0.25;
239  t_diff(0, 1, 0, 1) = 0.25;
240  t_diff(1, 0, 0, 1) = 0.25;
241  return t_diff;
242 };
243 
244 template <typename T1, typename T2>
245 inline auto deviator(Tensor2_symmetric<T1, 2> &t_stress, T2 &stress33) {
246  Tensor2_symmetric<double, 3> t_dev;
247  t_dev(I, J) = 0;
248  constexpr double third = boost::math::constants::third<double>();
249  const double trace = (t_stress(i, i) + stress33) * third;
250 
251  for (int ii = 0; ii != 2; ++ii)
252  for (int jj = ii; jj != 2; ++jj)
253  t_dev(ii, jj) = t_stress(ii, jj);
254 
255  t_dev(2, 2) = stress33; // for plane strain case sigma33 != 0
256  t_dev(0, 0) -= trace;
257  t_dev(1, 1) -= trace;
258  t_dev(2, 2) -= trace;
259  return t_dev;
260 };
261 
262 // inline auto deviator3(Tensor2_symmetric<double, 3> &t_stress) {
263 // Tensor2_symmetric<double, 3> t_dev;
264 // t_dev(I, J) = 0;
265 // const double lambda_p =
266 // poisson / ((1. + poisson) * (1. - 2. * poisson));
267 
268 // double stress33 = (t_stress(0, 0) + t_stress(1, 1)) * (1 - poisson);
269 // stress33 /= ((1. / lambda_p) + 2. * poisson);
270 // // stress33 = 0.;
271 // constexpr double third = boost::math::constants::third<double>();
272 // const double trace = (t_stress(i, i) + stress33) * third;
273 // for (int ii = 0; ii != 2; ++ii)
274 // for (int jj = ii; jj != 2; ++jj)
275 // t_dev(ii, jj) = t_stress(ii, jj);
276 
277 // t_dev(2, 2) = stress33;
278 // t_dev(0, 0) -= trace;
279 // t_dev(1, 1) -= trace;
280 // t_dev(2, 2) -= trace;
281 // return t_dev;
282 // };
283 
284 // inline Tensor2_symmetric<double, 3> getFiniteDiff(
285 // const int &k, const int &l,
286 // const Tensor2_symmetric<double, 3> &B) {
287 // MoFEMFunctionBeginHot;
288 
289 // Tensor2_symmetric<double, 3> out;
290 // Index<'I', 3> I;
291 // Index<'J', 3> J;
292 // Tensor2_symmetric<double, 3> stress = B;
293 
294 // Tensor2_symmetric<double, 3> strainPl = stress;
295 // Tensor2_symmetric<double, 3> strainMin = stress;
296 // const double h = 1e-8;
297 // // if (k != l) {
298 // // strainPl(k, l) = stress(k, l) + 0.5 * h;
299 // // strainMin(k, l) = stress(k, l) - 0.5 * h;
300 // // } else {
301 // strainPl(k, l) += h;
302 // strainMin(k, l) -= h;
303 // // }
304 
305 // auto StressPlusH = deviator3(strainPl);
306 // auto StressMinusH = deviator3(strainMin);
307 
308 // for (int i = 0; i != 3; ++i) {
309 // for (int j = 0; j != 3; ++j) {
310 // out(i, j) = (StressPlusH(i, j) - StressMinusH(i, j)) / (2 * h);
311 // }
312 // }
313 // return out;
314 // }
315 
316 // inline auto getDTensor4Function() {
317 // const Tensor2_symmetric<double, 3> stress(0., 0., 0., 0., 0., 0.);
318 // // Tensor4<double, 3, 3, 2, 2> my_D;
319 // Ddg<double, 3, 2> my_D;
320 // for (int k = 0; k != 2; ++k) {
321 // for (int l = 0; l != 2; ++l) {
322 // auto d_stress = getFiniteDiff(k, l, stress);
323 // for (int i = 0; i != 3; ++i) {
324 // for (int j = 0; j != 3; ++j) {
325 // my_D(i, j, k, l) = d_stress(i, j);
326 // }
327 // }
328 // }
329 // }
330 // return my_D;
331 // };
332 
333 inline auto diff_deviator(Ddg<double, 2, 2> &&t_diff_stress) {
334  Ddg<double, 3, 2> t_diff_deviator;
335  t_diff_deviator(I, J, k, l) = 0;
336  for (int ii = 0; ii != 2; ++ii)
337  for (int jj = ii; jj != 2; ++jj)
338  for (int kk = 0; kk != 2; ++kk)
339  for (int ll = kk; ll != 2; ++ll)
340  t_diff_deviator(ii, jj, kk, ll) = t_diff_stress(ii, jj, kk, ll);
341 
342  constexpr double third = boost::math::constants::third<double>();
343  if (!is_plane_strain) {
344 
345  t_diff_deviator(0, 0, 0, 0) -= third;
346  t_diff_deviator(0, 0, 1, 1) -= third;
347 
348  t_diff_deviator(1, 1, 0, 0) -= third;
349  t_diff_deviator(1, 1, 1, 1) -= third;
350 
351  t_diff_deviator(2, 2, 0, 0) -= third;
352  t_diff_deviator(2, 2, 1, 1) -= third;
353  }
354 
355  if (is_plane_strain) {
356  const double &nu = poisson;
357 
358  t_diff_deviator(0, 0, 0, 0) -= third * (1 + nu);
359  t_diff_deviator(0, 0, 1, 1) -= third * (1 + nu);
360 
361  t_diff_deviator(1, 1, 0, 0) -= third * (1 + nu);
362  t_diff_deviator(1, 1, 1, 1) -= third * (1 + nu);
363 
364  t_diff_deviator(2, 2, 0, 0) = nu - third * (1 + nu);
365  t_diff_deviator(2, 2, 1, 1) = nu - third * (1 + nu);
366 
367  // auto test_my_dev = getDTensor4Function();
368  // t_diff_deviator(I, J, k, l) = test_my_dev(I, J, k, l);
369  }
370 
371  return t_diff_deviator;
372 };
373 
374 template <typename T>
375 inline auto get_total_stress(Tensor2_symmetric<T, 2> &t_stress,
376  Tensor2_symmetric<T, 2> &t_plastic_strain) {
377 
378  Tensor2_symmetric<double, 2> stress_tmp;
379 
380  stress_tmp(i, j) = t_stress(i, j) - 1.5 * C1_k * t_plastic_strain(i, j);
381 
382  return stress_tmp;
383 };
384 
385 constexpr double sqrt2by3 = 0.816496580927726;
386 inline auto hardening(double tau) {
387  const double c = sqrt(2. / 3);
388  return H * tau * c + Q_inf * (1 - exp(-b_iso * c * tau)) + sigmaY;
389 }
390 
391 inline auto hardening_dtau(double tau) {
392  const double c = sqrt(2. / 3);
393  return H * c + Q_inf * b_iso * c * exp(-b_iso * c * tau);
394 }
395 
396 /**
397  *
398 
399 \f[
400 \begin{split}
401 f&=\sqrt{s_{ij}s_{ij}}\\
402 A_{ij}&=\frac{\partial f}{\partial \sigma_{ij}}=
403 \frac{1}{f} s_{kl} \frac{\partial s_{kl}}{\partial \sigma_{ij}}\\
404 \frac{\partial A_{ij}}{\partial \sigma_{kl}}&= \frac{\partial^2 f}{\partial
405 \sigma_{ij}\partial\sigma_{mn}}= \frac{1}{f} \left( \frac{\partial
406 s_{kl}}{\partial \sigma_{mn}}\frac{\partial s_{kl}}{\partial \sigma_{ij}}
407 -A_{mn}A_{ij}
408 \right)\\
409 \frac{\partial f}{\partial \varepsilon_{ij}}&=A_{mn}D_{mnij}
410 \\
411 \frac{\partial A_{ij}}{\partial \varepsilon_{kl}}&=
412 \frac{\partial A_{ij}}{\partial \sigma_{mn}} \frac{\partial
413 \sigma_{mn}}{\partial \varepsilon_{kl}}= \frac{\partial A_{ij}}{\partial
414 \sigma_{mn}} D_{mnkl}
415 \end{split}
416 \f]
417 
418  */
419 inline double
420 plastic_surface(Tensor2_symmetric<double, 3> &&t_stress_deviator) {
421  return std::sqrt(1.5 * t_stress_deviator(I, J) * t_stress_deviator(I, J));
422 };
423 
424 inline auto plastic_flow(double f, Tensor2_symmetric<double, 3> &&t_dev_stress,
425  Ddg<double, 3, 2> &t_diff_deviator) {
426  Tensor2_symmetric<double, 2> t_diff_f;
427  f += std::numeric_limits<double>::epsilon();
428  t_diff_f(k, l) =
429  ((1.5) * (t_dev_stress(I, J)) * t_diff_deviator(I, J, k, l)) / f;
430 
431  return t_diff_f;
432 };
433 
434 template <typename T>
435 inline auto diff_plastic_flow_dstress(double f, Tensor2_symmetric<T, 2> &t_flow,
436  Ddg<double, 3, 2> &t_diff_deviator) {
437  Ddg<double, 2, 2> t_diff_flow;
438  f += std::numeric_limits<double>::epsilon();
439  t_diff_flow(i, j, k, l) =
440  ((1.5) * (t_diff_deviator(M, N, i, j) * t_diff_deviator(M, N, k, l) -
441  2. / 3. * t_flow(i, j) * t_flow(k, l))) /
442  f;
443  return t_diff_flow;
444 };
445 
446 template <typename T>
447 inline auto
448 diff_plastic_flow_dstrain(Ddg<T, 2, 2> &t_D,
449  Ddg<double, 2, 2> &&t_diff_plastic_flow_dstress) {
450  Ddg<double, 2, 2> t_diff_flow;
451  t_diff_flow(i, j, k, l) =
452  t_diff_plastic_flow_dstress(i, j, m, n) * t_D(m, n, k, l);
453 
454  return t_diff_flow;
455 };
456 
457 template <typename T>
459  Ddg<T, 2, 2> &t_D, Ddg<double, 2, 2> &&t_diff_plastic_flow_dstress) {
460  Ddg<double, 2, 2> t_diff_flow_kin_hard;
461  t_diff_flow_kin_hard(i, j, k, l) =
462  t_diff_plastic_flow_dstress(i, j, m, n) * t_D(m, n, k, l) +
463  1.5 * C1_k * t_diff_plastic_flow_dstress(i, j, k, l);
464 
465  return t_diff_flow_kin_hard;
466 };
467 
468 // inline double getFiniteDiff2(const int &k, const int &l,
469 // const Tensor2_symmetric<double, 2> &B) {
470 // MoFEMFunctionBeginHot;
471 
472 // double out;
473 // Tensor2_symmetric<double, 2> stress = B;
474 // Tensor2_symmetric<double, 2> strainPl = stress;
475 // Tensor2_symmetric<double, 2> strainMin = stress;
476 // const double h = 1e-8;
477 
478 // strainPl(k, l) += h;
479 // strainMin(k, l) -= h;
480 // double stress33 = 0;
481 
482 // auto dev_plus = plastic_surface(deviator(strainPl, stress33));
483 // auto dev_minus = plastic_surface(deviator(strainMin, stress33));
484 
485 // out = (dev_plus - dev_minus) / (2 * h);
486 // return out;
487 // }
488 
489 // template <typename T>
490 // inline auto get2Tensor2Function(Tensor2_symmetric<T, 2> &stress_in) {
491 // Tensor2_symmetric<double, 2> stress;
492 // // Tensor4<double, 3, 3, 2, 2> my_D;
493 // stress(i, j) = stress_in(i, j);
494 // Tensor2_symmetric<double, 2> fdm_flow;
495 // for (int k = 0; k != 2; ++k) {
496 // for (int l = 0; l != 2; ++l) {
497 // double d_stress = getFiniteDiff2(k, l, stress);
498 // fdm_flow(k, l) = d_stress;
499 // }
500 // }
501 
502 // return fdm_flow;
503 // };
504 
505 /**
506 
507 \f[
508 \dot{\tau} - \frac{1}{2}\left\{\dot{\tau} + (f(\pmb\sigma) - \sigma_y) +
509 \| \dot{\tau} + (f(\pmb\sigma) - \sigma_y) \|\right\} = 0
510 \f]
511 
512  */
513 inline double contrains(double tau, double f) {
514  if ((f + cn_pl * tau) >= 0)
515  return -f;
516  else
517  return cn_pl * tau;
518 };
519 
520 inline double sign(double x) {
521  if (x == 0)
522  return 0;
523  else if (x > 0)
524  return 1;
525  else
526  return -1;
527 };
528 
529 inline double diff_constrain_dtau(double tau, double f) {
530  return (cn_pl - cn_pl * sign(f + cn_pl * tau)) / 2.;
531 };
532 
533 inline auto diff_constrain_df(double tau, double f) {
534  return (-1 - sign(f + cn_pl * tau)) / 2.;
535 };
536 
537 template <typename T>
539  Tensor2_symmetric<T, 2> &t_plastic_flow) {
540  Tensor2_symmetric<double, 2> t_diff_constrain_dstress;
541  t_diff_constrain_dstress(i, j) = diff_constrain_df * t_plastic_flow(i, j);
542  return t_diff_constrain_dstress;
543 };
544 
545 template <typename T>
546 inline auto
547 diff_constrain_dstrain(Ddg<T, 2, 2> &t_D,
548  Tensor2_symmetric<T, 2> &&t_diff_constrain_dstress) {
549  Tensor2_symmetric<double, 2> t_diff_constrain_dstrain;
550  t_diff_constrain_dstrain(k, l) =
551  t_diff_constrain_dstress(i, j) * t_D(i, j, k, l);
552  return t_diff_constrain_dstrain;
553 };
554 
555 template <typename T>
557  Ddg<T, 2, 2> &t_D, Tensor2_symmetric<T, 2> &&t_diff_constrain_dstress) {
558  Tensor2_symmetric<double, 2> t_diff_constrain_kin_hard_dstrain;
559  t_diff_constrain_kin_hard_dstrain(k, l) =
560  t_diff_constrain_dstress(i, j) * t_D(i, j, k, l) +
561  1.5 * C1_k * t_diff_constrain_dstress(k, l);
562  return t_diff_constrain_kin_hard_dstrain;
563 };
564 
567  auto diff_dev = diff_deviator(diff_tensor());
568  this->diffDeviator(I, J, k, l) = diff_dev(I, J, k, l);
570 }
571 
573  const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
574  : DomainEleOp(field_name, DomainEleOp::OPROW),
575  commonDataPtr(common_data_ptr) {
576  // Operator is only executed for vertices
577  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
578 }
579 
581  EntData &data) {
583 
584  const size_t nb_gauss_pts = commonDataPtr->mStressPtr->size2();
585  auto t_stress = getFTensor2SymmetricFromMat<2>(*(commonDataPtr->mStressPtr));
586  auto t_strain = getFTensor2SymmetricFromMat<2>(*(commonDataPtr->mStrainPtr));
587  auto &diff_dev = commonDataPtr->diffDeviator;
588  commonDataPtr->plasticSurfacePtr->resize(nb_gauss_pts, false);
589  commonDataPtr->plasticFlowPtr->resize(3, nb_gauss_pts, false);
590  auto t_flow =
591  getFTensor2SymmetricFromMat<2>(*(commonDataPtr->plasticFlowPtr));
592  auto t_plastic_strain =
593  getFTensor2SymmetricFromMat<2>(*(commonDataPtr->plasticStrainPtr));
594  const double lambda = young_modulus * poisson /
595  ((1. + poisson) * (1. - 2. * poisson));
596 
597  for (auto &f : *(commonDataPtr->plasticSurfacePtr)) {
598 
599  auto stress_tot = get_total_stress(t_stress, t_plastic_strain);
600  const double stress33 =
601  is_plane_strain * (stress_tot(0, 0) + stress_tot(1, 1)) * poisson;
602  Tensor2_symmetric<double, 2> el_strain;
603  el_strain(k, l) = t_strain(k, l) - t_plastic_strain(k, l);
604 
605  const double stress33_test =
606  (el_strain(0, 0) + el_strain(1, 1)) * lambda;
607 
608  f = plastic_surface(deviator(stress_tot, stress33));
609 
610  auto t_flow_tmp = plastic_flow(f, deviator(stress_tot, stress33), diff_dev);
611  t_flow(i, j) = t_flow_tmp(i, j);
612 
613  ++t_plastic_strain;
614  ++t_flow;
615  ++t_stress;
616  ++t_strain;
617  }
618 
620 }
621 
622 OpPlasticStress::OpPlasticStress(const std::string field_name,
623  boost::shared_ptr<CommonData> common_data_ptr)
624  : DomainEleOp(field_name, DomainEleOp::OPROW),
625  commonDataPtr(common_data_ptr) {
626  // Operator is only executed for vertices
627  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
628 }
629 
630 //! [Calculate stress]
631 MoFEMErrorCode OpPlasticStress::doWork(int side, EntityType type,
632  EntData &data) {
634  const size_t nb_gauss_pts = commonDataPtr->mStrainPtr->size2();
635  commonDataPtr->mStressPtr->resize(3, nb_gauss_pts);
636 
637  auto &t_D = commonDataPtr->tD;
638  auto t_strain = getFTensor2SymmetricFromMat<2>(*(commonDataPtr->mStrainPtr));
639  auto t_plastic_strain =
640  getFTensor2SymmetricFromMat<2>(*(commonDataPtr->plasticStrainPtr));
641 
642  auto t_stress = getFTensor2SymmetricFromMat<2>(*(commonDataPtr->mStressPtr));
643  plastic_strains = *(commonDataPtr->plasticStrainPtr); //FIXME:
644  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
645  t_stress(i, j) =
646  t_D(i, j, k, l) * (t_strain(k, l) - t_plastic_strain(k, l));
647 
648  ++t_strain;
649  ++t_plastic_strain;
650  ++t_stress;
651  }
652 
654 }
655 //! [Calculate stress]
656 
658  const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
659  : DomainEleOp(field_name, DomainEleOp::OPROW),
660  commonDataPtr(common_data_ptr) {}
661 
663  EntData &data) {
665  const size_t nb_dofs = data.getIndices().size();
666  if (nb_dofs) {
667 
668  auto get_dt = [&]() {
669  double dt;
670  CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
671  return dt;
672  };
673  const auto dt = get_dt();
674 
675  auto t_flow =
676  getFTensor2SymmetricFromMat<2>(*(commonDataPtr->plasticFlowPtr));
677  auto t_plastic_strain_dot =
678  getFTensor2SymmetricFromMat<2>(*(commonDataPtr->plasticStrainDotPtr));
679  auto t_tau_dot = getFTensor0FromVec(*(commonDataPtr->plasticTauDotPtr));
680  auto &t_D = commonDataPtr->tD;
681 
682  const size_t nb_integration_pts = data.getN().size1();
683  const size_t nb_base_functions = data.getN().size2();
684  auto t_w = getFTensor0IntegrationWeight();
685  auto t_base = data.getFTensor0N();
686 
687  std::array<double, MAX_DOFS_ON_ENTITY> nf;
688  std::fill(&nf[0], &nf[nb_dofs], 0);
689 
690  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
691  double alpha = dt * getMeasure() * t_w;
692 
693  Tensor2_symmetric<PackPtr<double *, 3>, 2> t_nf{&nf[0], &nf[1], &nf[2]};
694 
695  Tensor2_symmetric<double, 2> t_flow_stress_diff;
696  t_flow_stress_diff(i, j) = t_D(i, j, k, l) * (t_plastic_strain_dot(k, l) -
697  t_tau_dot * t_flow(k, l));
698 
699  size_t bb = 0;
700  for (; bb != nb_dofs / 3; ++bb) {
701  t_nf(i, j) += alpha * t_base * t_flow_stress_diff(i, j);
702  ++t_base;
703  ++t_nf;
704  }
705  for (; bb < nb_base_functions; ++bb)
706  ++t_base;
707 
708  ++t_flow;
709  ++t_plastic_strain_dot;
710  ++t_tau_dot;
711  ++t_w;
712  }
713 
714  CHKERR VecSetValues(getTSf(), data, nf.data(), ADD_VALUES);
715  }
716 
718 }
719 
721  const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
722  : DomainEleOp(field_name, DomainEleOp::OPROW),
723  commonDataPtr(common_data_ptr) {}
724 
726  EntData &data) {
728 
729  const size_t nb_dofs = data.getIndices().size();
730  if (nb_dofs) {
731 
732  auto get_dt = [&]() {
733  double dt;
734  CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
735  return dt;
736  };
737  const auto dt = get_dt();
738 
739  auto t_tau = getFTensor0FromVec(*(commonDataPtr->plasticTauPtr));
740  auto t_tau_dot = getFTensor0FromVec(*(commonDataPtr->plasticTauDotPtr));
741  auto t_f = getFTensor0FromVec(*(commonDataPtr->plasticSurfacePtr));
742  auto t_w = getFTensor0IntegrationWeight();
743 
744  std::array<double, MAX_DOFS_ON_ENTITY> nf;
745  std::fill(&nf[0], &nf[nb_dofs], 0);
746 
747  // SNES snes;
748  // CHKERR TSGetSNES(getFEMethod()->ts, &snes);
749  // CHKERR SNESGetIterationNumber(snes, &iter);
750 
751  auto t_base = data.getFTensor0N();
752  const size_t nb_integration_pts = data.getN().size1();
753  const size_t nb_base_functions = data.getN().size2();
754  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
755  const double alpha = dt * getMeasure() * t_w;
756  const double beta = alpha * contrains(t_tau_dot, t_f - hardening(t_tau));
757 
758  size_t bb = 0;
759  for (; bb != nb_dofs; ++bb) {
760  nf[bb] += beta * t_base;
761  ++t_base;
762  }
763  for (; bb < nb_base_functions; ++bb)
764  ++t_base;
765 
766  ++t_tau;
767  ++t_tau_dot;
768  ++t_f;
769  ++t_w;
770  }
771 
772  CHKERR VecSetValues(getTSf(), data, nf.data(), ADD_VALUES);
773  }
774 
776 }
777 
779  const std::string row_field_name, const std::string col_field_name,
780  boost::shared_ptr<CommonData> common_data_ptr)
781  : DomainEleOp(row_field_name, col_field_name, DomainEleOp::OPROWCOL),
782  commonDataPtr(common_data_ptr) {
783  sYmm = false;
784 }
785 
787  int row_side, int col_side, EntityType row_type, EntityType col_type,
788  EntData &row_data, EntData &col_data) {
790 
791  const size_t nb_row_dofs = row_data.getIndices().size();
792  const size_t nb_col_dofs = col_data.getIndices().size();
793  if (nb_row_dofs && nb_col_dofs) {
794 
795  locMat.resize(nb_row_dofs, nb_col_dofs, false);
796  locMat.clear();
797 
798  const size_t nb_integration_pts = row_data.getN().size1();
799  const size_t nb_row_base_functions = row_data.getN().size2();
800  auto t_w = getFTensor0IntegrationWeight();
801  auto t_row_diff_base = row_data.getFTensor1DiffN<2>();
802  auto &t_D = commonDataPtr->tD;
803  // auto F = getFTensor2FromMat<2,
804  // 2>(*(commonDataPtr->mDeformationGradient));
805 
806  MatrixDouble &dE_dF = *(commonDataPtr->dE_dF_mat);
807  Tensor4<PackPtr<double *, 1>, 2, 2, 2, 2> D2(MAT_TO_TENSOR4_2D(dE_dF));
808  auto t_diff_symmetrize = diff_symmetrize();
809 
810  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
811  double alpha = getMeasure() * t_w;
812 
813  size_t rr = 0;
814  for (; rr != nb_row_dofs / 2; ++rr) {
815 
816  Christof<PackPtr<double *, 3>, 2, 2> t_mat{
817 
818  &locMat(2 * rr + 0, 0), &locMat(2 * rr + 0, 1),
819  &locMat(2 * rr + 0, 2),
820 
821  &locMat(2 * rr + 1, 0), &locMat(2 * rr + 1, 1),
822  &locMat(2 * rr + 1, 2)
823 
824  };
825 
826  auto t_col_base = col_data.getFTensor0N(gg, 0);
827  for (size_t cc = 0; cc != nb_col_dofs / 3; ++cc) {
828 
829  Christof<double, 2, 2> t_tmp;
830  t_tmp(l, i, k) =
831  (t_D(i, j, k, l) * ((alpha * t_col_base) * t_row_diff_base(j)));
832 
833  if (is_large_strain) {
834  t_tmp(l, i, k) = 0.;
835  Tensor4<double, 2, 2, 2, 2> D2_symm;
836  D2_symm(i, j, k, l) = 0.;
837 
838  for (int i = 0; i != 2; ++i)
839  for (int j = 0; j != 2; ++j)
840  for (int k = 0; k != 2; ++k)
841  for (int l = 0; l != 2; ++l)
842  for (int m = 0; m != 2; ++m)
843  for (int n = 0; n != 2; ++n)
844  // D2_symm(i, j, k, l) +=
845  // D2(i, j, m, n) * t_diff_symmetrize(m, n, k, l);
846  D2_symm(k, l, m, n) += t_D(i, j, m, n) * D2(i, j, k, l);
847 
848  for (int i = 0; i != 2; ++i)
849  for (int j = 0; j != 2; ++j)
850  for (int k = 0; k != 2; ++k)
851  for (int l = 0; l != 2; ++l)
852  for (int m = 0; m != 2; ++m)
853  for (int n = 0; n != 2; ++n)
854  // t_tmp(k, m, n) +=
855  // (t_D(i, j, m, n) * D2_symm(i, j, k, l) *
856  // ((alpha * t_col_base) * t_row_diff_base(l)));
857  t_tmp(i, k, l) +=
858  (D2_symm(i, j, m, n) *
859  t_diff_symmetrize(m, n, k, l) *
860  ((alpha * t_col_base) * t_row_diff_base(j)));
861  }
862 
863  for (int ii = 0; ii != 2; ++ii)
864  for (int kk = 0; kk != 2; ++kk)
865  for (int ll = 0; ll != 2; ++ll)
866  t_mat(ii, kk, ll) -= t_tmp(ii, kk, ll);
867 
868  ++t_mat;
869  ++t_col_base;
870  }
871 
872  ++t_row_diff_base;
873  }
874 
875  if (is_large_strain)
876  ++D2;
877 
878  for (; rr < nb_row_base_functions; ++rr)
879  ++t_row_diff_base;
880  // ++F;
881  ++t_w;
882  }
883 
884  MatSetValues(getSNESB(), row_data, col_data, &*locMat.data().begin(),
885  ADD_VALUES);
886  }
887 
889 }
890 
892  const std::string row_field_name, const std::string col_field_name,
893  boost::shared_ptr<CommonData> common_data_ptr)
894  : DomainEleOp(row_field_name, col_field_name, DomainEleOp::OPROWCOL),
895  commonDataPtr(common_data_ptr) {
896  sYmm = false;
897 }
898 
899 MoFEMErrorCode OpCalculatePlasticFlowLhs_dU::doWork(int row_side, int col_side,
900  EntityType row_type,
901  EntityType col_type,
902  EntData &row_data,
903  EntData &col_data) {
905 
906  const size_t nb_row_dofs = row_data.getIndices().size();
907  const size_t nb_col_dofs = col_data.getIndices().size();
908  if (nb_row_dofs && nb_col_dofs) {
909 
910  auto get_dt = [&]() {
911  double dt;
912  CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
913  return dt;
914  };
915  const auto dt = get_dt();
916 
917  locMat.resize(nb_row_dofs, nb_col_dofs, false);
918  locMat.clear();
919 
920  const size_t nb_integration_pts = row_data.getN().size1();
921  const size_t nb_row_base_functions = row_data.getN().size2();
922  auto t_w = getFTensor0IntegrationWeight();
923  auto t_row_base = row_data.getFTensor0N();
924  auto t_f = getFTensor0FromVec(*(commonDataPtr->plasticSurfacePtr));
925  auto t_tau_dot = getFTensor0FromVec(*(commonDataPtr->plasticTauDotPtr));
926  auto &diff_dev = commonDataPtr->diffDeviator;
927  auto t_flow =
928  getFTensor2SymmetricFromMat<2>(*(commonDataPtr->plasticFlowPtr));
929  auto &t_D = commonDataPtr->tD;
930 
931  MatrixDouble &dE_dF = *(commonDataPtr->dE_dF_mat);
932  Tensor4<PackPtr<double *, 1>, 2, 2, 2, 2> D2(MAT_TO_TENSOR4_2D(dE_dF));
933 
934  auto t_diff_symmetrize = diff_symmetrize();
935 
936  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
937 
938  double alpha = dt * getMeasure() * t_w;
939  auto t_diff_plastic_flow_dstrain = diff_plastic_flow_dstrain(
940  t_D, diff_plastic_flow_dstress(t_f, t_flow, diff_dev));
941  Ddg<double, 2, 2> t_flow_stress_dstrain;
942  t_flow_stress_dstrain(i, j, k, l) =
943  t_D(i, j, m, n) * t_diff_plastic_flow_dstrain(m, n, k, l);
944  Tensor4<double, 2, 2, 2, 2> t_diff_plastic_flow_stress_dgrad;
945 
946  if (!is_large_strain)
947  t_diff_plastic_flow_stress_dgrad(i, j, k, l) =
948  t_flow_stress_dstrain(i, j, m, n) * t_diff_symmetrize(m, n, k, l);
949  else {
950 
951  Tensor4<double, 2, 2, 2, 2> t_diff_plastic_flow_stress_dgrads;
952  t_diff_plastic_flow_stress_dgrads(i, j, k, l) =
953  t_flow_stress_dstrain(i, j, m, n) * t_diff_symmetrize(m, n, k, l);
954  t_diff_plastic_flow_stress_dgrad(i, j, k, l) = 0.;
955  for (int i : {0, 1})
956  for (int j : {0, 1})
957  for (int k : {0, 1})
958  for (int l : {0, 1})
959  for (int m : {0, 1})
960  for (int n : {0, 1})
961  t_diff_plastic_flow_stress_dgrad(i, j, k, l) +=
962  t_diff_plastic_flow_stress_dgrads(i, j, m, n) *
963  D2(m, n, k, l);
964 
965  // Tensor4<double, 2, 2, 2, 2> tmp_t123;
966  // Tensor4<double, 2, 2, 2, 2> tmp_t4 =
967  // t_diff_plastic_flow_stress_dgrad;
968  // t_diff_plastic_flow_stress_dgrad(i, j, k, l) = 0.;
969  // for (int i : {0, 1})
970  // for (int j : {0, 1})
971  // for (int k : {0, 1})
972  // for (int l : {0, 1})
973  // for (int m : {0, 1})
974  // for (int n : {0, 1}) {
975 
976  // t_diff_plastic_flow_stress_dgrad(i, j, k, l) +=
977  // tmp_t4(i, j, m, n) * t_diff_symmetrize(m, n, k, l);
978 
979  // }
980  }
981 
982  size_t rr = 0;
983  for (; rr != nb_row_dofs / 3; ++rr) {
984 
985  Dg<PackPtr<double *, 2>, 2, 2> t_mat{&locMat(3 * rr + 0, 0),
986  &locMat(3 * rr + 0, 1),
987 
988  &locMat(3 * rr + 1, 0),
989  &locMat(3 * rr + 1, 1),
990 
991  &locMat(3 * rr + 2, 0),
992  &locMat(3 * rr + 2, 1)
993 
994  };
995 
996  const double c0 = alpha * t_row_base * t_tau_dot;
997 
998  auto t_col_diff_base = col_data.getFTensor1DiffN<2>(gg, 0);
999  for (size_t cc = 0; cc != nb_col_dofs / 2; ++cc) {
1000 
1001  Tensor3<double, 2, 2, 2> t_tmp;
1002  t_tmp(i, j, l) = c0 * (t_diff_plastic_flow_stress_dgrad(i, j, l, k) *
1003  t_col_diff_base(k));
1004 
1005  for (int ii = 0; ii != 2; ++ii)
1006  for (int jj = ii; jj < 2; ++jj)
1007  for (int ll = 0; ll != 2; ++ll)
1008  t_mat(ii, jj, ll) -= t_tmp(ii, jj, ll);
1009 
1010  ++t_mat;
1011  ++t_col_diff_base;
1012  }
1013 
1014  ++t_row_base;
1015  }
1016  for (; rr < nb_row_base_functions; ++rr)
1017  ++t_row_base;
1018 
1019  if (is_large_strain)
1020  ++D2;
1021 
1022  ++t_w;
1023  ++t_f;
1024  ++t_flow;
1025  ++t_tau_dot;
1026  }
1027 
1028  CHKERR MatSetValues(getSNESB(), row_data, col_data, &*locMat.data().begin(),
1029  ADD_VALUES);
1030  }
1031 
1033 }
1034 
1036  const std::string row_field_name, const std::string col_field_name,
1037  boost::shared_ptr<CommonData> common_data_ptr)
1038  : DomainEleOp(row_field_name, col_field_name, DomainEleOp::OPROWCOL),
1039  commonDataPtr(common_data_ptr) {
1040  sYmm = false;
1041 }
1042 
1043 MoFEMErrorCode OpCalculatePlasticFlowLhs_dEP::doWork(int row_side, int col_side,
1044  EntityType row_type,
1045  EntityType col_type,
1046  EntData &row_data,
1047  EntData &col_data) {
1049 
1050  const size_t nb_row_dofs = row_data.getIndices().size();
1051  const size_t nb_col_dofs = col_data.getIndices().size();
1052  if (nb_row_dofs && nb_col_dofs) {
1053 
1054  auto get_dt = [&]() {
1055  double dt;
1056  CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
1057  return dt;
1058  };
1059  const auto dt = get_dt();
1060 
1061  locMat.resize(nb_row_dofs, nb_col_dofs, false);
1062  locMat.clear();
1063 
1064  const size_t nb_integration_pts = row_data.getN().size1();
1065  const size_t nb_row_base_functions = row_data.getN().size2();
1066  auto t_w = getFTensor0IntegrationWeight();
1067  auto t_row_base = row_data.getFTensor0N();
1068  auto &diff_dev = commonDataPtr->diffDeviator;
1069  auto t_f = getFTensor0FromVec(*(commonDataPtr->plasticSurfacePtr));
1070  auto t_tau_dot = getFTensor0FromVec(*(commonDataPtr->plasticTauDotPtr));
1071  auto t_flow =
1072  getFTensor2SymmetricFromMat<2>(*(commonDataPtr->plasticFlowPtr));
1073 
1074  auto &t_D = commonDataPtr->tD;
1075  auto t_diff_plastic_strain = diff_tensor();
1076 
1077  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1078  double alpha = dt * getMeasure() * t_w;
1079  double beta = alpha * getTSa();
1080 
1081  size_t rr = 0;
1082  for (; rr != nb_row_dofs / 3; ++rr) {
1083 
1084  const double c0 = alpha * t_row_base * t_tau_dot;
1085  const double c1 = beta * t_row_base;
1086  auto t_diff_plastic_flow_dstrain = diff_plastic_flow_kin_hard_dstrain(
1087  t_D, diff_plastic_flow_dstress(t_f, t_flow, diff_dev));
1088 
1089  Ddg<PackPtr<double *, 3>, 2, 2> t_mat{
1090 
1091  &locMat(3 * rr + 0, 0), &locMat(3 * rr + 0, 1),
1092  &locMat(3 * rr + 0, 2),
1093 
1094  &locMat(3 * rr + 1, 0), &locMat(3 * rr + 1, 1),
1095  &locMat(3 * rr + 1, 2),
1096 
1097  &locMat(3 * rr + 2, 0), &locMat(3 * rr + 2, 1),
1098  &locMat(3 * rr + 2, 2)
1099 
1100  };
1101 
1102  auto t_col_base = col_data.getFTensor0N(gg, 0);
1103  for (size_t cc = 0; cc != nb_col_dofs / 3; ++cc) {
1104 
1105  t_mat(i, j, k, l) +=
1106  t_col_base * (t_D(i, j, m, n) *
1107  (c1 * t_diff_plastic_strain(m, n, k, l) +
1108  c0 * t_diff_plastic_flow_dstrain(m, n, k, l)));
1109 
1110  ++t_mat;
1111  ++t_col_base;
1112  }
1113 
1114  ++t_row_base;
1115  }
1116  for (; rr < nb_row_base_functions; ++rr)
1117  ++t_row_base;
1118 
1119  ++t_w;
1120  ++t_f;
1121  ++t_flow;
1122  ++t_tau_dot;
1123  }
1124 
1125  CHKERR MatSetValues(getTSB(), row_data, col_data, &*locMat.data().begin(),
1126  ADD_VALUES);
1127  }
1128 
1130 }
1131 
1133  const std::string row_field_name, const std::string col_field_name,
1134  boost::shared_ptr<CommonData> common_data_ptr)
1135  : DomainEleOp(row_field_name, col_field_name, DomainEleOp::OPROWCOL),
1136  commonDataPtr(common_data_ptr) {
1137  sYmm = false;
1138 }
1139 
1141 OpCalculatePlasticFlowLhs_dTAU::doWork(int row_side, int col_side,
1142  EntityType row_type, EntityType col_type,
1143  EntData &row_data, EntData &col_data) {
1145 
1146  const size_t nb_row_dofs = row_data.getIndices().size();
1147  const size_t nb_col_dofs = col_data.getIndices().size();
1148  if (nb_row_dofs && nb_col_dofs) {
1149 
1150  auto get_dt = [&]() {
1151  double dt;
1152  CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
1153  return dt;
1154  };
1155  const auto dt = get_dt();
1156 
1157  locMat.resize(nb_row_dofs, nb_col_dofs, false);
1158  locMat.clear();
1159 
1160  const size_t nb_integration_pts = row_data.getN().size1();
1161  auto t_w = getFTensor0IntegrationWeight();
1162  auto t_flow =
1163  getFTensor2SymmetricFromMat<2>(*(commonDataPtr->plasticFlowPtr));
1164 
1165  auto t_row_base = row_data.getFTensor0N();
1166 
1167  auto &t_D = commonDataPtr->tD;
1168 
1169  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1170  double alpha = dt * t_w * getMeasure() * getTSa();
1171 
1172  Tensor2_symmetric<double, 2> t_flow_stress;
1173  t_flow_stress(i, j) = t_D(i, j, m, n) * t_flow(m, n);
1174 
1175  for (size_t rr = 0; rr != nb_row_dofs / 3; ++rr) {
1176 
1177  Tensor2_symmetric<PackPtr<double *, 1>, 2> t_mat{
1178  &locMat(3 * rr + 0, 0), &locMat(3 * rr + 1, 0),
1179  &locMat(3 * rr + 2, 0)};
1180 
1181  auto t_col_base = col_data.getFTensor0N(gg, 0);
1182  for (size_t cc = 0; cc != nb_col_dofs; cc++) {
1183  t_mat(i, j) -= alpha * t_row_base * t_col_base * t_flow_stress(i, j);
1184  ++t_mat;
1185  ++t_col_base;
1186  }
1187 
1188  ++t_row_base;
1189  }
1190 
1191  ++t_w;
1192  ++t_flow;
1193  }
1194 
1195  CHKERR MatSetValues(getSNESB(), row_data, col_data, &*locMat.data().begin(),
1196  ADD_VALUES);
1197  }
1198 
1200 }
1201 
1203  const std::string row_field_name, const std::string col_field_name,
1204  boost::shared_ptr<CommonData> common_data_ptr)
1205  : DomainEleOp(row_field_name, col_field_name, DomainEleOp::OPROWCOL),
1206  commonDataPtr(common_data_ptr) {
1207  sYmm = false;
1208 }
1209 
1210 MoFEMErrorCode OpCalculateContrainsLhs_dU::doWork(int row_side, int col_side,
1211  EntityType row_type,
1212  EntityType col_type,
1213  EntData &row_data,
1214  EntData &col_data) {
1216 
1217  const size_t nb_row_dofs = row_data.getIndices().size();
1218  const size_t nb_col_dofs = col_data.getIndices().size();
1219  if (nb_row_dofs && nb_col_dofs) {
1220 
1221  auto get_dt = [&]() {
1222  double dt;
1223  CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
1224  return dt;
1225  };
1226  const auto dt = get_dt();
1227 
1228  locMat.resize(nb_row_dofs, nb_col_dofs, false);
1229  locMat.clear();
1230  const size_t nb_integration_pts = row_data.getN().size1();
1231  const size_t nb_row_base_functions = row_data.getN().size2();
1232  auto t_w = getFTensor0IntegrationWeight();
1233  auto t_row_base = row_data.getFTensor0N();
1234  auto t_f = getFTensor0FromVec(*(commonDataPtr->plasticSurfacePtr));
1235  auto t_tau = getFTensor0FromVec(*(commonDataPtr->plasticTauPtr));
1236  auto t_tau_dot = getFTensor0FromVec(*(commonDataPtr->plasticTauDotPtr));
1237  auto t_flow =
1238  getFTensor2SymmetricFromMat<2>(*(commonDataPtr->plasticFlowPtr));
1239  auto t_stress =
1240  getFTensor2SymmetricFromMat<2>(*(commonDataPtr->mStressPtr));
1241  auto &t_D = commonDataPtr->tD;
1242  auto t_diff_symmetrize = diff_symmetrize();
1243 
1244  // Tensor2<double, 2, 2> t_test_grad_u(2., 4., 1., -1.);
1245  // Tensor2<double, 2, 2> test1;
1246  // Tensor2<double, 2, 2> test2;
1247  // Tensor2<double, 2, 2> test3;
1248  // Tensor2<double, 2, 2> test4;
1249  // Tensor2<double, 2, 2> test5;
1250  // test1(i, j) = Is(i, j, k, l) * t_test_grad_u(k, l);
1251  // test3(i, j) = 0.5 * (t_test_grad_u(i, j) + t_test_grad_u(j, i));
1252  // test2(i, j) = t_diff_symmetrize(i, j, k, l) * t_test_grad_u(k, l);
1253  // test4(k, l) = t_diff_symmetrize(i, j, k, l) * t_test_grad_u(i, j);
1254  // test4(k, l) = Is(i, j, k, l) * t_test_grad_u(i, j);
1255  // bool _is_IS_ddg = is_ddg(Is);
1256 
1257  MatrixDouble &dE_dF = *(commonDataPtr->dE_dF_mat);
1258  Tensor4<PackPtr<double *, 1>, 2, 2, 2, 2> D2(MAT_TO_TENSOR4_2D(dE_dF));
1259 
1260  // if (getTStime() > 0.16)
1261  // string wait;
1262 
1263  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1264  double alpha = dt * getMeasure() * t_w;
1265 
1266  auto t_diff_constrain_dstrain = diff_constrain_dstrain(
1267  t_D,
1269  diff_constrain_df(t_tau_dot, t_f - hardening(t_tau)), t_flow));
1270 
1271  // if (t_diff_constrain_dstrain(0, 0) != 0)
1272  // string wait;
1273 
1274  Tensor2<double, 2, 2> t_diff_constrain_dgrad;
1275  if (!is_large_strain) {
1276 
1277  t_diff_constrain_dgrad(k, l) =
1278  t_diff_constrain_dstrain(i, j) * t_diff_symmetrize(i, j, k, l);
1279  } else {
1280 
1281  Tensor2<double, 2, 2> t_diff_constrain_dgrads;
1282  t_diff_constrain_dgrads(k, l) =
1283  t_diff_constrain_dstrain(i, j) * t_diff_symmetrize(i, j, k, l);
1284 
1285  t_diff_constrain_dgrad(k, l) =
1286  t_diff_constrain_dgrads(i, j) * D2(i, j, k, l);
1287  // Tensor2<double, 2, 2> tmp_tt = t_diff_constrain_dgrad;
1288  // t_diff_constrain_dgrad(k, l) =
1289  // tmp_tt(i, j) * t_diff_symmetrize(i, j, k, l);
1290  }
1291 
1292  Tensor1<PackPtr<double *, 2>, 2> t_mat{&locMat(0, 0), &locMat(0, 1)};
1293 
1294  size_t rr = 0;
1295  for (; rr != nb_row_dofs; ++rr) {
1296 
1297  auto t_col_diff_base = col_data.getFTensor1DiffN<2>(gg, 0);
1298  for (size_t cc = 0; cc != nb_col_dofs / 2; cc++) {
1299 
1300  t_mat(i) += alpha * t_row_base * t_diff_constrain_dgrad(i, j) *
1301  t_col_diff_base(j);
1302 
1303  ++t_mat;
1304  ++t_col_diff_base;
1305  }
1306 
1307  ++t_row_base;
1308  }
1309  for (; rr != nb_row_base_functions; ++rr)
1310  ++t_row_base;
1311 
1312  if (is_large_strain)
1313  ++D2;
1314 
1315  ++t_f;
1316  ++t_tau;
1317  ++t_tau_dot;
1318  ++t_flow;
1319  ++t_stress;
1320  ++t_w;
1321  }
1322 
1323  CHKERR MatSetValues(getSNESB(), row_data, col_data, &*locMat.data().begin(),
1324  ADD_VALUES);
1325  }
1326 
1328 }
1329 
1331  const std::string row_field_name, const std::string col_field_name,
1332  boost::shared_ptr<CommonData> common_data_ptr)
1333  : DomainEleOp(row_field_name, col_field_name, DomainEleOp::OPROWCOL),
1334  commonDataPtr(common_data_ptr) {
1335  sYmm = false;
1336 }
1337 
1338 MoFEMErrorCode OpCalculateContrainsLhs_dEP::doWork(int row_side, int col_side,
1339  EntityType row_type,
1340  EntityType col_type,
1341  EntData &row_data,
1342  EntData &col_data) {
1344 
1345  const size_t nb_row_dofs = row_data.getIndices().size();
1346  const size_t nb_col_dofs = col_data.getIndices().size();
1347  if (nb_row_dofs && nb_col_dofs) {
1348 
1349  auto get_dt = [&]() {
1350  double dt;
1351  CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
1352  return dt;
1353  };
1354  const auto dt = get_dt();
1355 
1356  locMat.resize(nb_row_dofs, nb_col_dofs, false);
1357  locMat.clear();
1358 
1359  const size_t nb_integration_pts = row_data.getN().size1();
1360  const size_t nb_row_base_functions = row_data.getN().size2();
1361  auto t_w = getFTensor0IntegrationWeight();
1362  auto t_row_base = row_data.getFTensor0N();
1363  auto t_f = getFTensor0FromVec(*(commonDataPtr->plasticSurfacePtr));
1364  auto t_tau = getFTensor0FromVec(*(commonDataPtr->plasticTauPtr));
1365  auto t_tau_dot = getFTensor0FromVec(*(commonDataPtr->plasticTauDotPtr));
1366  auto t_flow =
1367  getFTensor2SymmetricFromMat<2>(*(commonDataPtr->plasticFlowPtr));
1368  auto &t_D = commonDataPtr->tD;
1369  auto t_diff_symmetrize = diff_symmetrize();
1370 
1371  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1372  double alpha = dt * getMeasure() * t_w;
1373 
1374  auto mat_ptr = locMat.data().begin();
1375  auto t_diff_constrain_dstrain = diff_constrain_kin_hard_dstrain(
1376  t_D,
1378  diff_constrain_df(t_tau_dot, t_f - hardening(t_tau)), t_flow));
1379 
1380  Tensor2_symmetric<PackPtr<double *, 3>, 2> t_mat{
1381  &locMat(0, 0), &locMat(0, 1), &locMat(0, 2)};
1382 
1383  size_t rr = 0;
1384  for (; rr != nb_row_dofs; ++rr) {
1385 
1386  auto t_col_base = col_data.getFTensor0N(gg, 0);
1387  for (size_t cc = 0; cc != nb_col_dofs / 3; cc++) {
1388 
1389  t_mat(i, j) -=
1390  alpha * t_row_base * t_col_base * t_diff_constrain_dstrain(i, j);
1391 
1392  ++t_mat;
1393  ++t_col_base;
1394  }
1395 
1396  ++t_row_base;
1397  }
1398  for (; rr != nb_row_base_functions; ++rr)
1399  ++t_row_base;
1400 
1401  ++t_f;
1402  ++t_tau;
1403  ++t_tau_dot;
1404  ++t_flow;
1405  ++t_w;
1406  }
1407 
1408  CHKERR MatSetValues(getSNESB(), row_data, col_data, &*locMat.data().begin(),
1409  ADD_VALUES);
1410  }
1411 
1413 }
1414 
1416  const std::string row_field_name, const std::string col_field_name,
1417  boost::shared_ptr<CommonData> common_data_ptr)
1418  : DomainEleOp(row_field_name, col_field_name, DomainEleOp::OPROWCOL),
1419  commonDataPtr(common_data_ptr) {
1420  sYmm = false;
1421 }
1422 
1423 MoFEMErrorCode OpCalculateContrainsLhs_dTAU::doWork(int row_side, int col_side,
1424  EntityType row_type,
1425  EntityType col_type,
1426  EntData &row_data,
1427  EntData &col_data) {
1429 
1430  const size_t nb_row_dofs = row_data.getIndices().size();
1431  const size_t nb_col_dofs = col_data.getIndices().size();
1432  if (nb_row_dofs && nb_col_dofs) {
1433 
1434  auto get_dt = [&]() {
1435  double dt;
1436  CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
1437  return dt;
1438  };
1439  const auto dt = get_dt();
1440 
1441  locMat.resize(nb_row_dofs, nb_col_dofs, false);
1442  locMat.clear();
1443 
1444  const size_t nb_integration_pts = row_data.getN().size1();
1445  const size_t nb_row_base_functions = row_data.getN().size2();
1446  auto t_w = getFTensor0IntegrationWeight();
1447  auto t_f = getFTensor0FromVec(*(commonDataPtr->plasticSurfacePtr));
1448  auto t_tau = getFTensor0FromVec(*(commonDataPtr->plasticTauPtr));
1449  auto t_tau_dot = getFTensor0FromVec(*(commonDataPtr->plasticTauDotPtr));
1450 
1451  auto t_row_base = row_data.getFTensor0N();
1452  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1453  const double alpha = dt * getMeasure() * t_w;
1454  const double c0 = alpha * getTSa() *
1455  diff_constrain_dtau(t_tau_dot, t_f - hardening(t_tau));
1456  const double c1 = alpha *
1457  diff_constrain_df(t_tau_dot, t_f - hardening(t_tau)) *
1458  hardening_dtau(t_tau);
1459 
1460  auto mat_ptr = locMat.data().begin();
1461 
1462  size_t rr = 0;
1463  for (; rr != nb_row_dofs; ++rr) {
1464 
1465  auto t_col_base = col_data.getFTensor0N(gg, 0);
1466  for (size_t cc = 0; cc != nb_col_dofs; ++cc) {
1467  *mat_ptr += (c0 - c1) * t_row_base * t_col_base;
1468  ++mat_ptr;
1469  ++t_col_base;
1470  }
1471  ++t_row_base;
1472  }
1473  for (; rr < nb_row_base_functions; ++rr)
1474  ++t_row_base;
1475 
1476  ++t_w;
1477  ++t_f;
1478  ++t_tau;
1479  ++t_tau_dot;
1480  }
1481 
1482  CHKERR MatSetValues(getSNESB(), row_data, col_data, &*locMat.data().begin(),
1483  ADD_VALUES);
1484  }
1485 
1487 }
1488 
1490  const std::string field_name, moab::Interface &post_proc_mesh,
1491  std::vector<EntityHandle> &map_gauss_pts,
1492  boost::shared_ptr<CommonData> common_data_ptr)
1493  : DomainEleOp(field_name, DomainEleOp::OPROW), postProcMesh(post_proc_mesh),
1494  mapGaussPts(map_gauss_pts), commonDataPtr(common_data_ptr) {
1495  // Operator is only executed for vertices
1496  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1497 }
1498 
1499 //! [Postprocessing]
1500 MoFEMErrorCode OpPostProcPlastic::doWork(int side, EntityType type,
1501  EntData &data) {
1503 
1504  std::array<double, 9> def;
1505  std::fill(def.begin(), def.end(), 0);
1506 
1507  auto get_tag = [&](const std::string name, size_t size) {
1508  Tag th;
1509  CHKERR postProcMesh.tag_get_handle(name.c_str(), size, MB_TYPE_DOUBLE, th,
1510  MB_TAG_CREAT | MB_TAG_SPARSE,
1511  def.data());
1512  return th;
1513  };
1514 
1515  MatrixDouble3by3 mat(3, 3);
1516 
1517  auto set_matrix_2d = [&](auto &t) -> MatrixDouble3by3 & {
1518  mat.clear();
1519  for (size_t r = 0; r != 2; ++r)
1520  for (size_t c = 0; c != 2; ++c)
1521  mat(r, c) = t(r, c);
1522  return mat;
1523  };
1524 
1525  auto set_matrix_3d = [&](auto &t) -> MatrixDouble3by3 & {
1526  mat.clear();
1527  for (size_t r = 0; r != 3; ++r)
1528  for (size_t c = 0; c != 3; ++c)
1529  mat(r, c) = t(r, c);
1530  return mat;
1531  };
1532 
1533  auto set_matrix_2d_symm = [&](auto &t) -> MatrixDouble3by3 & {
1534  mat.clear();
1535  for (size_t r = 0; r != 2; ++r)
1536  for (size_t c = 0; c <= r; ++c)
1537  mat(r, c) = t(r, c);
1538  return mat;
1539  };
1540 
1541  auto set_scalar = [&](auto t) -> MatrixDouble3by3 & {
1542  mat.clear();
1543  mat(0, 0) = t;
1544  return mat;
1545  };
1546 
1547  auto set_tag = [&](auto th, auto gg, MatrixDouble3by3 &mat) {
1548  return postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
1549  &*mat.data().begin());
1550  };
1551 
1552  auto th_plastic_surface = get_tag("PLASTIC_SURFACE", 1);
1553  auto th_tau = get_tag("PLASTIC_MULTIPLIER", 1);
1554 
1555  auto th_log_plastic_strain = get_tag("PLASTIC_STRAIN_LOG", 9);
1556  auto th_eqv_ep = get_tag("EQUIV_EP", 1);
1557  auto th_det_g = get_tag("DET_G", 1);
1558  auto th_trace_ep = get_tag("TRACE_EP", 1);
1559  auto th_hardening = get_tag("HARDENING", 1);
1560 
1561  auto th_plastic_flow = get_tag("PLASTIC_FLOW", 9);
1562  auto th_plastic_flow_trace = get_tag("PLASTIC_FLOW_TRACE", 1);
1563  string strain_name = "PLASTIC_STRAIN";
1564  if (is_large_strain)
1565  strain_name = "PLASTIC_DEFORMATION";
1566 
1567  auto th_plastic_strain = get_tag(strain_name, 9);
1568 
1569  auto t_flow =
1570  getFTensor2SymmetricFromMat<2>(*(commonDataPtr->plasticFlowPtr));
1571  auto t_plastic_strain =
1572  getFTensor2SymmetricFromMat<2>(*(commonDataPtr->plasticStrainPtr));
1573  auto t_stress = getFTensor2SymmetricFromMat<2>(*(commonDataPtr->mStressPtr));
1574 
1575  size_t gg = 0;
1576  for (; gg != commonDataPtr->plasticSurfacePtr->size(); ++gg) {
1577  const double f = (*(commonDataPtr->plasticSurfacePtr))[gg];
1578  const double tau = (*(commonDataPtr->plasticTauPtr))[gg];
1579 
1580 
1581  CHKERR set_tag(th_plastic_surface, gg, set_scalar(f - hardening(tau)));
1582  CHKERR set_tag(th_tau, gg, set_scalar(tau));
1583  CHKERR set_tag(th_plastic_flow, gg, set_matrix_2d_symm(t_flow));
1584  if (is_large_strain) {
1585  const double flow_trace = t_flow(i, i);
1586  CHKERR set_tag(th_plastic_flow_trace, gg, set_scalar(flow_trace));
1587 
1588  CHKERR set_tag(th_log_plastic_strain, gg,
1589  set_matrix_2d_symm(t_plastic_strain));
1590  t_plastic_strain(i, j) *= 2.;
1591  auto exp_plast_strain = exp_map(t_plastic_strain);
1592  double det_g;
1593  CHKERR determinantTensor2by2(exp_plast_strain, det_g);
1594  CHKERR set_tag(th_det_g, gg, set_scalar(det_g));
1595 
1596  const double trace_ep = t_plastic_strain(i, i);
1597  CHKERR set_tag(th_trace_ep, gg, set_scalar(trace_ep));
1598  CHKERR set_tag(th_plastic_strain, gg,
1599  set_matrix_2d_symm(exp_plast_strain));
1600  const double eqv_ep = std::sqrt(t_plastic_strain(i, j) * t_plastic_strain(i, j));
1601  CHKERR set_tag(th_eqv_ep, gg, set_scalar(eqv_ep));
1602  const double hard = hardening(tau) / young_modulus_inv;
1603  CHKERR set_tag(th_hardening, gg, set_scalar(hard));
1604  } else
1605  CHKERR set_tag(th_plastic_strain, gg,
1606  set_matrix_2d_symm(t_plastic_strain));
1607 
1608  ++t_flow;
1609  ++t_stress;
1610  ++t_plastic_strain;
1611  }
1612 
1614 }
1615 
1617 
1618  OpGetGaussPtsPlasticState(const std::string field_name,
1619  boost::shared_ptr<CommonData> common_data_ptr)
1620  : DomainEleOp(field_name, DomainEleOp::OPROW),
1621  commonDataPtr(common_data_ptr) {
1622  // Operator is only executed for vertices
1623  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1624  }
1625 
1626  //! [Calculate stress]
1627  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
1629  // const size_t nb_gauss_pts = commonDataPtr->mStrainPtr->size2();
1630  auto t_tau = getFTensor0FromVec(*(commonDataPtr->plasticTauPtr));
1631  vector<double> data_vec{0, 0}; // active and all
1632 
1633  for (auto &f : *(commonDataPtr->plasticSurfacePtr)) {
1634 
1635  // if ((f + cn_pl * t_tau) >= 0)
1636  if (cn_pl * t_tau <= f)
1637  data_vec[0] += 1;
1638  data_vec[1] += 1;
1639 
1640  ++t_tau;
1641  }
1642 
1643  constexpr std::array<int, 2> indices = {0, 1};
1644  CHKERR VecSetValues(commonDataPtr->stateVecPlast, 2, indices.data(),
1645  &data_vec[0], ADD_VALUES);
1646 
1648  }
1649 
1650 private:
1651  boost::shared_ptr<CommonData> commonDataPtr;
1652 };
1653 
1654 //! [Postprocessing]
1655 }; // namespace OpPlasticTools
OpPostProcPlastic(const std::string field_name, moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, boost::shared_ptr< CommonData > common_data_ptr)
OpCalculatePlasticFlowLhs_dEP(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: PlasticOps.hpp:750
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Definition: PlasticOps.hpp:758
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
Definition: PlasticOps.hpp:912
constexpr double sqrt2by3
static FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:143
boost::shared_ptr< VectorDouble > plasticTauDotPtr
Definition: PlasticOps.hpp:57
OpCalculateContrainsRhs(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: PlasticOps.hpp:528
OpCalculateContrainsLhs_dU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: PlasticOps.hpp:904
boost::shared_ptr< MatrixDouble > plasticFlowPtr
Definition: PlasticOps.hpp:55
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: PlasticOps.hpp:411
FTensor::Index< 'm', 2 > m
Definition: PlasticOps.hpp:67
double diff_constrain_dtau(double tau, double f)
Definition: PlasticOps.hpp:375
boost::shared_ptr< MatrixDouble > plasticStrainPtr
Definition: PlasticOps.hpp:58
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:134
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:173
boost::shared_ptr< VectorDouble > plasticTauPtr
Definition: PlasticOps.hpp:56
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Postprocessing]
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:100
auto hardening(double tau)
Definition: PlasticOps.hpp:281
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Definition: PlasticOps.hpp:583
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
FTensor::Index< 'J', 3 > J
Definition: PlasticOps.hpp:71
double plastic_surface(Tensor2_symmetric< double, 3 > &&t_stress_deviator)
FTensor::Tensor1< double *, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
OpPlasticStress(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: PlasticOps.hpp:440
constexpr double H
OpCalculatePlasticFlowLhs_dU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: PlasticOps.hpp:653
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
Definition: PlasticOps.hpp:661
OpCalculatePlasticFlowLhs_dTAU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: PlasticOps.hpp:841
FTensor::Index< 'I', 3 > I
Definition: PlasticOps.hpp:70
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Calculate stress]
FTensor::Index< 'N', 3 > N
Definition: PlasticOps.hpp:73
std::vector< EntityHandle > & mapGaussPts
Definition: PlasticOps.hpp:212
auto diff_constrain_dstress(double &&diff_constrain_df, FTensor::Tensor2_symmetric< T, 2 > &t_plastic_flow)
Definition: PlasticOps.hpp:385
MoFEMErrorCode calculateDiffDeviator()
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Definition: PlasticOps.hpp:995
MoFEMErrorCode determinantTensor2by2(T1 &t, T2 &det)
Calculate determinant 2 by 2.
Definition: Templates.hpp:584
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
boost::shared_ptr< CommonData > commonDataPtr
OpCalculateContrainsLhs_dEP(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: PlasticOps.hpp:987
const VectorInt & getIndices() const
Get global indices of dofs on entity.
double contrains(double tau, double f)
Definition: PlasticOps.hpp:359
double sign(double x)
Definition: PlasticOps.hpp:366
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: PlasticOps.hpp:533
Definition: single.cpp:3
boost::shared_ptr< MatrixDouble > plasticStrainDotPtr
Definition: PlasticOps.hpp:59
OpCalculatePlasticSurface(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
[Lambda functions]
Definition: PlasticOps.hpp:403
FTensor::Index< 'i', 2 > i
[Common data]
Definition: PlasticOps.hpp:63
auto diff_tensor()
[Operators definitions]
Definition: PlasticOps.hpp:218
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
Definition: PlasticOps.hpp:850
auto plastic_flow(double f, FTensor::Tensor2_symmetric< double, 3 > &&t_dev_stress, FTensor::Ddg< double, 3, 2 > &&t_diff_deviator)
Definition: PlasticOps.hpp:313
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:186
auto getFTensor0IntegrationWeight()
Get integration weights.
double trace(FTensor::Tensor2_symmetric< T, 2 > &t_stress)
Definition: PlasticOps.hpp:240
auto hardening_dtau()
Definition: PlasticOps.hpp:283
auto diff_deviator(FTensor::Ddg< double, 2, 2 > &&t_diff_stress)
Definition: PlasticOps.hpp:258
auto diff_constrain_dstrain(FTensor::Ddg< T, 2, 2 > &t_D, FTensor::Tensor2_symmetric< T, 2 > &&t_diff_constrain_dstress)
Definition: PlasticOps.hpp:393
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:102
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
constexpr double young_modulus
#define CHKERR
Inline error check.
Definition: definitions.h:604
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:82
const double r
rate factor
auto diff_plastic_flow_dstrain(FTensor::Ddg< T, 2, 2 > &t_D, FTensor::Ddg< double, 2, 2 > &&t_diff_plastic_flow_dstress)
Definition: PlasticOps.hpp:341
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:91
auto exp_map(const Tensor2_symmetric< T, 3 > &X)
Index< 'L', 3 > L
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:147
boost::shared_ptr< VectorDouble > plasticSurfacePtr
Definition: PlasticOps.hpp:54
OpCalculatePlasticFlowRhs(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
[Calculate stress]
Definition: PlasticOps.hpp:471
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: PlasticOps.hpp:476
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Calculate stress]
Definition: PlasticOps.hpp:449
OpGetGaussPtsPlasticState(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
auto diff_plastic_flow_kin_hard_dstrain(Ddg< T, 2, 2 > &t_D, Ddg< double, 2, 2 > &&t_diff_plastic_flow_dstress)
auto diff_constrain_kin_hard_dstrain(Ddg< T, 2, 2 > &t_D, Tensor2_symmetric< T, 2 > &&t_diff_constrain_dstress)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1943
Data on single entity (This is passed as argument to DataOperator::doWork)
OpCalculateContrainsLhs_dTAU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
auto get_total_stress(Tensor2_symmetric< T, 2 > &t_stress, Tensor2_symmetric< T, 2 > &t_plastic_strain)
auto diff_symmetrize()
Definition: PlasticOps.hpp:227
constexpr double third
FTensor::Index< 'n', 2 > n
Definition: PlasticOps.hpp:68
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Ddg< double, 3, 2 > diffDeviator
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:199
constexpr double sigmaY
OpCalculatePlasticInternalForceLhs_dEP(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: PlasticOps.hpp:575
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:109
FTensor::Index< 'k', 2 > k
Definition: PlasticOps.hpp:65
Index< 'K', 3 > K
auto diff_plastic_flow_dstress(double f, FTensor::Tensor2_symmetric< T, 2 > &t_flow, FTensor::Ddg< double, 3, 2 > &&t_diff_deviator)
Definition: PlasticOps.hpp:327
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:213
FTensor::Index< 'l', 2 > l
Definition: PlasticOps.hpp:66
moab::Interface & postProcMesh
Definition: PlasticOps.hpp:211
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:160
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
auto deviator(FTensor::Tensor2_symmetric< T, 2 > &t_stress, double trace)
Definition: PlasticOps.hpp:246
FTensor::Index< 'j', 2 > j
Definition: PlasticOps.hpp:64
FTensor::Index< 'M', 3 > M
Definition: PlasticOps.hpp:72
auto diff_constrain_df(double tau, double f)
Definition: PlasticOps.hpp:379