v0.9.1
PlasticOps.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]
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 //! [Common data]
62 
69 
74 
75 //! [Operators definitions]
77  OpCalculatePlasticSurface(const std::string field_name,
78  boost::shared_ptr<CommonData> common_data_ptr);
79  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
80 
81 private:
82  boost::shared_ptr<CommonData> commonDataPtr;
83 };
84 
85 struct OpPlasticStress : public DomainEleOp {
86  OpPlasticStress(const std::string field_name,
87  boost::shared_ptr<CommonData> common_data_ptr);
88  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
89 
90 private:
91  boost::shared_ptr<CommonData> commonDataPtr;
92 };
93 
95  OpCalculatePlasticFlowRhs(const std::string field_name,
96  boost::shared_ptr<CommonData> common_data_ptr);
97  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
98 
99 private:
100  boost::shared_ptr<CommonData> commonDataPtr;
101 };
102 
104  OpCalculateContrainsRhs(const std::string field_name,
105  boost::shared_ptr<CommonData> common_data_ptr);
106  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
107 
108 private:
109  boost::shared_ptr<CommonData> commonDataPtr;
110 };
111 
114  const std::string row_field_name, const std::string col_field_name,
115  boost::shared_ptr<CommonData> common_data_ptr);
116  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
117  EntityType col_type, EntData &row_data,
118  EntData &col_data);
119 
120 private:
121  boost::shared_ptr<CommonData> commonDataPtr;
123 };
124 
126  OpCalculatePlasticFlowLhs_dU(const std::string row_field_name,
127  const std::string col_field_name,
128  boost::shared_ptr<CommonData> common_data_ptr);
129  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
130  EntityType col_type, EntData &row_data,
131  EntData &col_data);
132 
133 private:
134  boost::shared_ptr<CommonData> commonDataPtr;
136 };
137 
139  OpCalculatePlasticFlowLhs_dEP(const std::string row_field_name,
140  const std::string col_field_name,
141  boost::shared_ptr<CommonData> common_data_ptr);
142  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
143  EntityType col_type, EntData &row_data,
144  EntData &col_data);
145 
146 private:
147  boost::shared_ptr<CommonData> commonDataPtr;
149 };
150 
152  OpCalculatePlasticFlowLhs_dTAU(const std::string row_field_name,
153  const std::string col_field_name,
154  boost::shared_ptr<CommonData> common_data_ptr);
155  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
156  EntityType col_type, EntData &row_data,
157  EntData &col_data);
158 
159 private:
160  boost::shared_ptr<CommonData> commonDataPtr;
162 };
163 
165  OpCalculateContrainsLhs_dU(const std::string row_field_name,
166  const std::string col_field_name,
167  boost::shared_ptr<CommonData> common_data_ptr);
168  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
169  EntityType col_type, EntData &row_data,
170  EntData &col_data);
171 
172 private:
173  boost::shared_ptr<CommonData> commonDataPtr;
175 };
176 
178  OpCalculateContrainsLhs_dEP(const std::string row_field_name,
179  const std::string col_field_name,
180  boost::shared_ptr<CommonData> common_data_ptr);
181  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
182  EntityType col_type, EntData &row_data,
183  EntData &col_data);
184 
185 private:
186  boost::shared_ptr<CommonData> commonDataPtr;
188 };
189 
191  OpCalculateContrainsLhs_dTAU(const std::string row_field_name,
192  const std::string col_field_name,
193  boost::shared_ptr<CommonData> common_data_ptr);
194  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
195  EntityType col_type, EntData &row_data,
196  EntData &col_data);
197 
198 private:
199  boost::shared_ptr<CommonData> commonDataPtr;
201 };
202 
204  OpPostProcPlastic(const std::string field_name,
205  moab::Interface &post_proc_mesh,
206  std::vector<EntityHandle> &map_gauss_pts,
207  boost::shared_ptr<CommonData> common_data_ptr);
208  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
209 
210 private:
212  std::vector<EntityHandle> &mapGaussPts;
213  boost::shared_ptr<CommonData> commonDataPtr;
214 };
215 //! [Operators definitions]
216 
217 //! [Lambda functions]
218 inline auto diff_tensor() {
220  t_diff(i, j, k, l) = 0;
221  for (size_t ii = 0; ii != 2; ++ii)
222  for (size_t jj = ii; jj != 2; ++jj)
223  t_diff(ii, jj, ii, jj) = 1;
224  return t_diff;
225 };
226 
227 inline auto diff_symmetrize() {
229  t_diff(i, j, k, l) = 0;
230  t_diff(0, 0, 0, 0) = 1;
231  t_diff(1, 1, 1, 1) = 1;
232  t_diff(1, 0, 1, 0) = 0.25;
233  t_diff(0, 1, 1, 0) = 0.25;
234  t_diff(0, 1, 0, 1) = 0.25;
235  t_diff(1, 0, 0, 1) = 0.25;
236  return t_diff;
237 };
238 
239 template <typename T>
240 inline double trace(FTensor::Tensor2_symmetric<T, 2> &t_stress) {
241  constexpr double third = boost::math::constants::third<double>();
242  return (t_stress(0, 0) + t_stress(1, 1)) * third;
243 };
244 
245 template <typename T>
246 inline auto deviator(FTensor::Tensor2_symmetric<T, 2> &t_stress, double trace) {
248  t_dev(I, J) = 0;
249  for (int ii = 0; ii != 2; ++ii)
250  for (int jj = ii; jj != 2; ++jj)
251  t_dev(ii, jj) = t_stress(ii, jj);
252  t_dev(0, 0) -= trace;
253  t_dev(1, 1) -= trace;
254  t_dev(2, 2) -= trace;
255  return t_dev;
256 };
257 
258 inline auto diff_deviator(FTensor::Ddg<double, 2, 2> &&t_diff_stress) {
259  FTensor::Ddg<double, 3, 2> t_diff_deviator;
260  t_diff_deviator(I, J, k, l) = 0;
261  for (int ii = 0; ii != 2; ++ii)
262  for (int jj = ii; jj != 2; ++jj)
263  for (int kk = 0; kk != 2; ++kk)
264  for (int ll = kk; ll != 2; ++ll)
265  t_diff_deviator(ii, jj, kk, ll) = t_diff_stress(ii, jj, kk, ll);
266 
267  constexpr double third = boost::math::constants::third<double>();
268 
269  t_diff_deviator(0, 0, 0, 0) -= third;
270  t_diff_deviator(0, 0, 1, 1) -= third;
271 
272  t_diff_deviator(1, 1, 0, 0) -= third;
273  t_diff_deviator(1, 1, 1, 1) -= third;
274 
275  t_diff_deviator(2, 2, 0, 0) -= third;
276  t_diff_deviator(2, 2, 1, 1) -= third;
277 
278  return t_diff_deviator;
279 };
280 
281 inline auto hardening(double tau) { return H * tau + sigmaY; }
282 
283 inline auto hardening_dtau() { return H; }
284 
285 /**
286  *
287 
288 \f[
289 \begin{split}
290 f&=\sqrt{s_{ij}s_{ij}}\\
291 A_{ij}&=\frac{\partial f}{\partial \sigma_{ij}}=
292 \frac{1}{f} s_{kl} \frac{\partial s_{kl}}{\partial \sigma_{ij}}\\
293 \frac{\partial A_{ij}}{\partial \sigma_{kl}}&= \frac{\partial^2 f}{\partial
294 \sigma_{ij}\partial\sigma_{mn}}= \frac{1}{f} \left( \frac{\partial
295 s_{kl}}{\partial \sigma_{mn}}\frac{\partial s_{kl}}{\partial \sigma_{ij}}
296 -A_{mn}A_{ij}
297 \right)\\
298 \frac{\partial f}{\partial \varepsilon_{ij}}&=A_{mn}D_{mnij}
299 \\
300 \frac{\partial A_{ij}}{\partial \varepsilon_{kl}}&=
301 \frac{\partial A_{ij}}{\partial \sigma_{mn}} \frac{\partial
302 \sigma_{mn}}{\partial \varepsilon_{kl}}= \frac{\partial A_{ij}}{\partial
303 \sigma_{mn}} D_{mnkl}
304 \end{split}
305 \f]
306 
307  */
308 inline double
310  return std::sqrt(1.5 * t_stress_deviator(I, J) * t_stress_deviator(I, J));
311 };
312 
313 inline auto plastic_flow(double f,
315  FTensor::Ddg<double, 3, 2> &&t_diff_deviator) {
317  if (std::abs(f) > std::numeric_limits<double>::epsilon())
318  t_diff_f(k, l) =
319  (1.5 / f) * (t_dev_stress(I, J) * t_diff_deviator(I, J, k, l));
320  else
321  t_diff_f(k, l) = 0;
322  return t_diff_f;
323 };
324 
325 template <typename T>
326 inline auto
328  FTensor::Ddg<double, 3, 2> &&t_diff_deviator) {
329  FTensor::Ddg<double, 2, 2> t_diff_flow;
330  if (std::abs(f) > std::numeric_limits<double>::epsilon())
331  t_diff_flow(i, j, k, l) =
332  (1.5 / f) * (t_diff_deviator(M, N, i, j) * t_diff_deviator(M, N, k, l) -
333  2. / 3. * t_flow(i, j) * t_flow(k, l));
334  else
335  t_diff_flow(i, j, k, l) = 0;
336 
337  return t_diff_flow;
338 };
339 
340 template <typename T>
343  FTensor::Ddg<double, 2, 2> &&t_diff_plastic_flow_dstress) {
344  FTensor::Ddg<double, 2, 2> t_diff_flow;
345  t_diff_flow(i, j, k, l) =
346  t_diff_plastic_flow_dstress(i, j, m, n) * t_D(m, n, k, l);
347 
348  return t_diff_flow;
349 };
350 
351 /**
352 
353 \f[
354 \dot{\tau} - \frac{1}{2}\left\{\dot{\tau} + (f(\pmb\sigma) - \sigma_y) +
355 \| \dot{\tau} + (f(\pmb\sigma) - \sigma_y) \|\right\} = 0
356 \f]
357 
358  */
359 inline double contrains(double tau, double f) {
360  if ((f + cn * tau) >= 0)
361  return -f;
362  else
363  return cn * tau;
364 };
365 
366 inline double sign(double x) {
367  if (x == 0)
368  return 0;
369  else if (x > 0)
370  return 1;
371  else
372  return -1;
373 };
374 
375 inline double diff_constrain_dtau(double tau, double f) {
376  return (cn - cn * sign(f + cn * tau)) / 2.;
377 };
378 
379 inline auto diff_constrain_df(double tau, double f) {
380  return (-1 - sign(f + cn * tau)) / 2.;
381 };
382 
383 template <typename T>
384 inline auto
386  FTensor::Tensor2_symmetric<T, 2> &t_plastic_flow) {
387  FTensor::Tensor2_symmetric<double, 2> t_diff_constrain_dstress;
388  t_diff_constrain_dstress(i, j) = diff_constrain_df * t_plastic_flow(i, j);
389  return t_diff_constrain_dstress;
390 };
391 
392 template <typename T>
395  FTensor::Tensor2_symmetric<T, 2> &&t_diff_constrain_dstress) {
396  FTensor::Tensor2_symmetric<double, 2> t_diff_constrain_dstrain;
397  t_diff_constrain_dstrain(k, l) =
398  t_diff_constrain_dstress(i, j) * t_D(i, j, k, l);
399  return t_diff_constrain_dstrain;
400 };
401 //! [Lambda functions]
402 
404  const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
405  : DomainEleOp(field_name, DomainEleOp::OPROW),
406  commonDataPtr(common_data_ptr) {
407  // Opetor is only executed for vertices
408  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
409 }
410 
412  EntData &data) {
414 
415  const size_t nb_gauss_pts = commonDataPtr->mStressPtr->size2();
416  auto t_stress = getFTensor2SymmetricFromMat<2>(*(commonDataPtr->mStressPtr));
417  auto t_strain = getFTensor2SymmetricFromMat<2>(*(commonDataPtr->mStrainPtr));
418 
419  commonDataPtr->plasticSurfacePtr->resize(nb_gauss_pts, false);
420  commonDataPtr->plasticFlowPtr->resize(3, nb_gauss_pts, false);
421  auto t_flow =
422  getFTensor2SymmetricFromMat<2>(*(commonDataPtr->plasticFlowPtr));
423 
424  for (auto &f : *(commonDataPtr->plasticSurfacePtr)) {
425  f = platsic_surface(deviator(t_stress, trace(t_stress)));
426  auto t_flow_tmp = plastic_flow(f,
427 
428  deviator(t_stress, trace(t_stress)),
429 
431  t_flow(i, j) = t_flow_tmp(i, j);
432  ++t_flow;
433  ++t_stress;
434  ++t_strain;
435  }
436 
438 }
439 
440 OpPlasticStress::OpPlasticStress(const std::string field_name,
441  boost::shared_ptr<CommonData> common_data_ptr)
442  : DomainEleOp(field_name, DomainEleOp::OPROW),
443  commonDataPtr(common_data_ptr) {
444  // Opetor is only executed for vertices
445  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
446 }
447 
448 //! [Calculate stress]
450  EntData &data) {
452  const size_t nb_gauss_pts = commonDataPtr->mStrainPtr->size2();
453  commonDataPtr->mStressPtr->resize(3, nb_gauss_pts);
454  auto &t_D = commonDataPtr->tD;
455  auto t_strain = getFTensor2SymmetricFromMat<2>(*(commonDataPtr->mStrainPtr));
456  auto t_plastic_strain =
457  getFTensor2SymmetricFromMat<2>(*(commonDataPtr->plasticStrainPtr));
458  auto t_stress = getFTensor2SymmetricFromMat<2>(*(commonDataPtr->mStressPtr));
459 
460  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
461  t_stress(i, j) = t_D(i, j, k, l) * (t_strain(k, l) - t_plastic_strain(k, l));
462  ++t_strain;
463  ++t_plastic_strain;
464  ++t_stress;
465  }
466 
468 }
469 //! [Calculate stress]
470 
472  const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
473  : DomainEleOp(field_name, DomainEleOp::OPROW),
474  commonDataPtr(common_data_ptr) {}
475 
477  EntData &data) {
479  const size_t nb_dofs = data.getIndices().size();
480  if (nb_dofs) {
481 
482  auto t_flow =
483  getFTensor2SymmetricFromMat<2>(*(commonDataPtr->plasticFlowPtr));
484  auto t_plastic_strain_dot =
485  getFTensor2SymmetricFromMat<2>(*(commonDataPtr->plasticStrainDotPtr));
486  auto t_tau_dot = getFTensor0FromVec(*(commonDataPtr->plasticTauDotPtr));
487  auto &t_D = commonDataPtr->tD;
488 
489  const size_t nb_integration_pts = data.getN().size1();
490  const size_t nb_base_functions = data.getN().size2();
491  auto t_w = getFTensor0IntegrationWeight();
492  auto t_base = data.getFTensor0N();
493 
494  std::array<double, MAX_DOFS_ON_ENTITY> nf;
495  std::fill(&nf[0], &nf[nb_dofs], 0);
496 
497  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
498  double alpha = getMeasure() * t_w;
499 
501  &nf[0], &nf[1], &nf[2]};
502 
503  FTensor::Tensor2_symmetric<double, 2> t_flow_stress_diff;
504  t_flow_stress_diff(i, j) = t_D(i, j, k, l) * (t_plastic_strain_dot(k, l) -
505  t_tau_dot * t_flow(k, l));
506 
507  size_t bb = 0;
508  for (; bb != nb_dofs / 3; ++bb) {
509  t_nf(i, j) += alpha * t_base * t_flow_stress_diff(i, j);
510  ++t_base;
511  ++t_nf;
512  }
513  for (; bb < nb_base_functions; ++bb)
514  ++t_base;
515 
516  ++t_flow;
517  ++t_plastic_strain_dot;
518  ++t_tau_dot;
519  ++t_w;
520  }
521 
522  CHKERR VecSetValues(getTSf(), data, nf.data(), ADD_VALUES);
523  }
524 
526 }
527 
529  const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
530  : DomainEleOp(field_name, DomainEleOp::OPROW),
531  commonDataPtr(common_data_ptr) {}
532 
534  EntData &data) {
536 
537  const size_t nb_dofs = data.getIndices().size();
538  if (nb_dofs) {
539 
540  auto t_tau = getFTensor0FromVec(*(commonDataPtr->plasticTauPtr));
541  auto t_tau_dot = getFTensor0FromVec(*(commonDataPtr->plasticTauDotPtr));
542  auto t_f = getFTensor0FromVec(*(commonDataPtr->plasticSurfacePtr));
543  auto t_w = getFTensor0IntegrationWeight();
544 
545  std::array<double, MAX_DOFS_ON_ENTITY> nf;
546  std::fill(&nf[0], &nf[nb_dofs], 0);
547 
548  auto t_base = data.getFTensor0N();
549  const size_t nb_integration_pts = data.getN().size1();
550  const size_t nb_base_functions = data.getN().size2();
551  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
552  const double alpha = getMeasure() * t_w;
553  const double beta = alpha * contrains(t_tau_dot, t_f - hardening(t_tau));
554 
555  size_t bb = 0;
556  for (; bb != nb_dofs; ++bb) {
557  nf[bb] += beta * t_base;
558  ++t_base;
559  }
560  for (; bb < nb_base_functions; ++bb)
561  ++t_base;
562 
563  ++t_tau;
564  ++t_tau_dot;
565  ++t_f;
566  ++t_w;
567  }
568 
569  CHKERR VecSetValues(getTSf(), data, nf.data(), ADD_VALUES);
570  }
571 
573 }
574 
576  const std::string row_field_name, const std::string col_field_name,
577  boost::shared_ptr<CommonData> common_data_ptr)
578  : DomainEleOp(row_field_name, col_field_name, DomainEleOp::OPROWCOL),
579  commonDataPtr(common_data_ptr) {
580  sYmm = false;
581 }
582 
584  int row_side, int col_side, EntityType row_type, EntityType col_type,
585  EntData &row_data, EntData &col_data) {
587 
588  const size_t nb_row_dofs = row_data.getIndices().size();
589  const size_t nb_col_dofs = col_data.getIndices().size();
590  if (nb_row_dofs && nb_col_dofs) {
591 
592  locMat.resize(nb_row_dofs, nb_col_dofs, false);
593  locMat.clear();
594 
595  const size_t nb_integration_pts = row_data.getN().size1();
596  const size_t nb_row_base_functions = row_data.getN().size2();
597  auto t_w = getFTensor0IntegrationWeight();
598  auto t_row_diff_base = row_data.getFTensor1DiffN<2>();
599  auto &t_D = commonDataPtr->tD;
600 
601  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
602  double alpha = getMeasure() * t_w;
603 
604  size_t rr = 0;
605  for (; rr != nb_row_dofs / 2; ++rr) {
606 
608 
609  &locMat(2 * rr + 0, 0), &locMat(2 * rr + 0, 1),
610  &locMat(2 * rr + 0, 2),
611 
612  &locMat(2 * rr + 1, 0), &locMat(2 * rr + 1, 1),
613  &locMat(2 * rr + 1, 2)
614 
615  };
616 
617  auto t_col_base = col_data.getFTensor0N(gg, 0);
618  for (size_t cc = 0; cc != nb_col_dofs / 3; ++cc) {
619 
620  // I mix up the indices here so that it behaves like a
621  // Dg. That way I don't have to have a separate wrapper
622  // class Christof_Expr, which simplifies things.
623  // You cyclicly has to shift index, i, j, k -> l, i, j
625  t_tmp(l, i, k) =
626  (t_D(i, j, k, l) * ((alpha * t_col_base) * t_row_diff_base(j)));
627 
628  for (int ii = 0; ii != 2; ++ii)
629  for (int kk = 0; kk != 2; ++kk)
630  for (int ll = 0; ll != 2; ++ll)
631  t_mat(ii, kk, ll) -= t_tmp(ii, kk, ll);
632 
633  ++t_mat;
634  ++t_col_base;
635  }
636 
637  ++t_row_diff_base;
638  }
639 
640  for (; rr < nb_row_base_functions; ++rr)
641  ++t_row_diff_base;
642 
643  ++t_w;
644  }
645 
646  MatSetValues(getSNESB(), row_data, col_data, &*locMat.data().begin(),
647  ADD_VALUES);
648  }
649 
651 }
652 
654  const std::string row_field_name, const std::string col_field_name,
655  boost::shared_ptr<CommonData> common_data_ptr)
656  : DomainEleOp(row_field_name, col_field_name, DomainEleOp::OPROWCOL),
657  commonDataPtr(common_data_ptr) {
658  sYmm = false;
659 }
660 
662  EntityType row_type,
663  EntityType col_type,
664  EntData &row_data,
665  EntData &col_data) {
667 
668  const size_t nb_row_dofs = row_data.getIndices().size();
669  const size_t nb_col_dofs = col_data.getIndices().size();
670  if (nb_row_dofs && nb_col_dofs) {
671 
672  locMat.resize(nb_row_dofs, nb_col_dofs, false);
673  locMat.clear();
674 
675  const size_t nb_integration_pts = row_data.getN().size1();
676  const size_t nb_row_base_functions = row_data.getN().size2();
677  auto t_w = getFTensor0IntegrationWeight();
678  auto t_row_base = row_data.getFTensor0N();
679  auto t_f = getFTensor0FromVec(*(commonDataPtr->plasticSurfacePtr));
680  auto t_tau_dot = getFTensor0FromVec(*(commonDataPtr->plasticTauDotPtr));
681  auto t_flow =
682  getFTensor2SymmetricFromMat<2>(*(commonDataPtr->plasticFlowPtr));
683  auto &t_D = commonDataPtr->tD;
684  auto t_diff_symmetrize = diff_symmetrize();
685 
686  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
687 
688  double alpha = getMeasure() * t_w;
689  auto t_diff_plastic_flow_dstrain = diff_plastic_flow_dstrain(
690  t_D,
692  FTensor::Ddg<double, 2, 2> t_flow_stress_dstrain;
693  t_flow_stress_dstrain(i, j, k, l) =
694  t_D(i, j, m, n) * t_diff_plastic_flow_dstrain(m, n, k, l);
695  FTensor::Tensor4<double, 2, 2, 2, 2> t_diff_plastic_flow_stress_dgrad;
696  t_diff_plastic_flow_stress_dgrad(i, j, k, l) =
697  t_flow_stress_dstrain(i, j, m, n) * t_diff_symmetrize(m, n, k, l);
698 
699  size_t rr = 0;
700  for (; rr != nb_row_dofs / 3; ++rr) {
701 
703  &locMat(3 * rr + 0, 0),
704  &locMat(3 * rr + 0, 1),
705 
706  &locMat(3 * rr + 1, 0),
707  &locMat(3 * rr + 1, 1),
708 
709  &locMat(3 * rr + 2, 0),
710  &locMat(3 * rr + 2, 1)
711 
712  };
713 
714  const double c0 = alpha * t_row_base * t_tau_dot;
715 
716  auto t_col_diff_base = col_data.getFTensor1DiffN<2>(gg, 0);
717  for (size_t cc = 0; cc != nb_col_dofs / 2; ++cc) {
718 
720  t_tmp(i, j, l) = c0 * (t_diff_plastic_flow_stress_dgrad(i, j, l, k) *
721  t_col_diff_base(k));
722 
723  for (int ii = 0; ii != 2; ++ii)
724  for (int jj = ii; jj < 2; ++jj)
725  for (int ll = 0; ll != 2; ++ll)
726  t_mat(ii, jj, ll) -= t_tmp(ii, jj, ll);
727 
728  ++t_mat;
729  ++t_col_diff_base;
730  }
731 
732  ++t_row_base;
733  }
734  for (; rr < nb_row_base_functions; ++rr)
735  ++t_row_base;
736 
737  ++t_w;
738  ++t_f;
739  ++t_flow;
740  ++t_tau_dot;
741  }
742 
743  CHKERR MatSetValues(getSNESB(), row_data, col_data, &*locMat.data().begin(),
744  ADD_VALUES);
745  }
746 
748 }
749 
751  const std::string row_field_name, const std::string col_field_name,
752  boost::shared_ptr<CommonData> common_data_ptr)
753  : DomainEleOp(row_field_name, col_field_name, DomainEleOp::OPROWCOL),
754  commonDataPtr(common_data_ptr) {
755  sYmm = false;
756 }
757 
759  EntityType row_type,
760  EntityType col_type,
761  EntData &row_data,
762  EntData &col_data) {
764 
765  const size_t nb_row_dofs = row_data.getIndices().size();
766  const size_t nb_col_dofs = col_data.getIndices().size();
767  if (nb_row_dofs && nb_col_dofs) {
768 
769  locMat.resize(nb_row_dofs, nb_col_dofs, false);
770  locMat.clear();
771 
772  const size_t nb_integration_pts = row_data.getN().size1();
773  const size_t nb_row_base_functions = row_data.getN().size2();
774  auto t_w = getFTensor0IntegrationWeight();
775  auto t_row_base = row_data.getFTensor0N();
776  auto t_f = getFTensor0FromVec(*(commonDataPtr->plasticSurfacePtr));
777  auto t_tau_dot = getFTensor0FromVec(*(commonDataPtr->plasticTauDotPtr));
778  auto t_flow =
779  getFTensor2SymmetricFromMat<2>(*(commonDataPtr->plasticFlowPtr));
780 
781  auto &t_D = commonDataPtr->tD;
782  auto t_diff_plastic_strain = diff_tensor();
783 
784  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
785  double alpha = getMeasure() * t_w;
786  double beta = alpha * getTSa();
787 
788  size_t rr = 0;
789  for (; rr != nb_row_dofs / 3; ++rr) {
790 
791  const double c0 = alpha * t_row_base * t_tau_dot;
792  const double c1 = beta * t_row_base;
793 
794  auto t_diff_plastic_flow_dstrain = diff_plastic_flow_dstrain(
795  t_D, diff_plastic_flow_dstress(t_f, t_flow,
797 
799 
800  &locMat(3 * rr + 0, 0), &locMat(3 * rr + 0, 1),
801  &locMat(3 * rr + 0, 2),
802 
803  &locMat(3 * rr + 1, 0), &locMat(3 * rr + 1, 1),
804  &locMat(3 * rr + 1, 2),
805 
806  &locMat(3 * rr + 2, 0), &locMat(3 * rr + 2, 1),
807  &locMat(3 * rr + 2, 2)
808 
809  };
810 
811  auto t_col_base = col_data.getFTensor0N(gg, 0);
812  for (size_t cc = 0; cc != nb_col_dofs / 3; ++cc) {
813 
814  t_mat(i, j, k, l) +=
815  t_col_base * (t_D(i, j, m, n) *
816  (c1 * t_diff_plastic_strain(m, n, k, l) +
817  c0 * t_diff_plastic_flow_dstrain(m, n, k, l)));
818 
819  ++t_mat;
820  ++t_col_base;
821  }
822 
823  ++t_row_base;
824  }
825  for (; rr < nb_row_base_functions; ++rr)
826  ++t_row_base;
827 
828  ++t_w;
829  ++t_f;
830  ++t_flow;
831  ++t_tau_dot;
832  }
833 
834  CHKERR MatSetValues(getTSB(), row_data, col_data, &*locMat.data().begin(),
835  ADD_VALUES);
836  }
837 
839 }
840 
842  const std::string row_field_name, const std::string col_field_name,
843  boost::shared_ptr<CommonData> common_data_ptr)
844  : DomainEleOp(row_field_name, col_field_name, DomainEleOp::OPROWCOL),
845  commonDataPtr(common_data_ptr) {
846  sYmm = false;
847 }
848 
850 OpCalculatePlasticFlowLhs_dTAU::doWork(int row_side, int col_side,
851  EntityType row_type, EntityType col_type,
852  EntData &row_data, EntData &col_data) {
854 
855  const size_t nb_row_dofs = row_data.getIndices().size();
856  const size_t nb_col_dofs = col_data.getIndices().size();
857  if (nb_row_dofs && nb_col_dofs) {
858 
859  locMat.resize(nb_row_dofs, nb_col_dofs, false);
860  locMat.clear();
861 
862  const size_t nb_integration_pts = row_data.getN().size1();
863  auto t_w = getFTensor0IntegrationWeight();
864  auto t_flow =
865  getFTensor2SymmetricFromMat<2>(*(commonDataPtr->plasticFlowPtr));
866 
867  auto t_row_base = row_data.getFTensor0N();
868 
869  auto &t_D = commonDataPtr->tD;
870 
871  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
872  double alpha = getMeasure() * t_w * getTSa();
873 
875  t_flow_stress(i, j) = t_D(i, j, m, n) * t_flow(m, n);
876 
877  for (size_t rr = 0; rr != nb_row_dofs / 3; ++rr) {
878 
880  &locMat(3 * rr + 0, 0), &locMat(3 * rr + 1, 0),
881  &locMat(3 * rr + 2, 0)};
882 
883  auto t_col_base = col_data.getFTensor0N(gg, 0);
884  for (size_t cc = 0; cc != nb_col_dofs; cc++) {
885  t_mat(i, j) -= alpha * t_row_base * t_col_base * t_flow_stress(i, j);
886  ++t_mat;
887  ++t_col_base;
888  }
889 
890  ++t_row_base;
891  }
892 
893  ++t_w;
894  ++t_flow;
895  }
896 
897  CHKERR MatSetValues(getSNESB(), row_data, col_data, &*locMat.data().begin(),
898  ADD_VALUES);
899  }
900 
902 }
903 
905  const std::string row_field_name, const std::string col_field_name,
906  boost::shared_ptr<CommonData> common_data_ptr)
907  : DomainEleOp(row_field_name, col_field_name, DomainEleOp::OPROWCOL),
908  commonDataPtr(common_data_ptr) {
909  sYmm = false;
910 }
911 
913  EntityType row_type,
914  EntityType col_type,
915  EntData &row_data,
916  EntData &col_data) {
918 
919  const size_t nb_row_dofs = row_data.getIndices().size();
920  const size_t nb_col_dofs = col_data.getIndices().size();
921  if (nb_row_dofs && nb_col_dofs) {
922 
923  locMat.resize(nb_row_dofs, nb_col_dofs, false);
924  locMat.clear();
925 
926  const size_t nb_integration_pts = row_data.getN().size1();
927  const size_t nb_row_base_functions = row_data.getN().size2();
928  auto t_w = getFTensor0IntegrationWeight();
929  auto t_row_base = row_data.getFTensor0N();
930  auto t_f = getFTensor0FromVec(*(commonDataPtr->plasticSurfacePtr));
931  auto t_tau = getFTensor0FromVec(*(commonDataPtr->plasticTauPtr));
932  auto t_tau_dot = getFTensor0FromVec(*(commonDataPtr->plasticTauDotPtr));
933  auto t_flow =
934  getFTensor2SymmetricFromMat<2>(*(commonDataPtr->plasticFlowPtr));
935  auto t_stress =
936  getFTensor2SymmetricFromMat<2>(*(commonDataPtr->mStressPtr));
937  auto &t_D = commonDataPtr->tD;
938  auto t_diff_symmetrize = diff_symmetrize();
939 
940  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
941  double alpha = getMeasure() * t_w;
942 
943  auto t_diff_constrain_dstrain = diff_constrain_dstrain(
944  t_D,
946  diff_constrain_df(t_tau_dot, t_f - hardening(t_tau)), t_flow));
947  FTensor::Tensor2<double, 2, 2> t_diff_constrain_dgrad;
948  t_diff_constrain_dgrad(k, l) =
949  t_diff_constrain_dstrain(i, j) * t_diff_symmetrize(i, j, k, l);
950 
952  &locMat(0, 1)};
953 
954  size_t rr = 0;
955  for (; rr != nb_row_dofs; ++rr) {
956 
957  auto t_col_diff_base = col_data.getFTensor1DiffN<2>(gg, 0);
958  for (size_t cc = 0; cc != nb_col_dofs / 2; cc++) {
959 
960  t_mat(i) += alpha * t_row_base * t_diff_constrain_dgrad(i, j) *
961  t_col_diff_base(j);
962 
963  ++t_mat;
964  ++t_col_diff_base;
965  }
966 
967  ++t_row_base;
968  }
969  for (; rr != nb_row_base_functions; ++rr)
970  ++t_row_base;
971 
972  ++t_f;
973  ++t_tau;
974  ++t_tau_dot;
975  ++t_flow;
976  ++t_stress;
977  ++t_w;
978  }
979 
980  CHKERR MatSetValues(getSNESB(), row_data, col_data, &*locMat.data().begin(),
981  ADD_VALUES);
982  }
983 
985 }
986 
988  const std::string row_field_name, const std::string col_field_name,
989  boost::shared_ptr<CommonData> common_data_ptr)
990  : DomainEleOp(row_field_name, col_field_name, DomainEleOp::OPROWCOL),
991  commonDataPtr(common_data_ptr) {
992  sYmm = false;
993 }
994 
996  EntityType row_type,
997  EntityType col_type,
998  EntData &row_data,
999  EntData &col_data) {
1001 
1002  const size_t nb_row_dofs = row_data.getIndices().size();
1003  const size_t nb_col_dofs = col_data.getIndices().size();
1004  if (nb_row_dofs && nb_col_dofs) {
1005 
1006  locMat.resize(nb_row_dofs, nb_col_dofs, false);
1007  locMat.clear();
1008 
1009  const size_t nb_integration_pts = row_data.getN().size1();
1010  const size_t nb_row_base_functions = row_data.getN().size2();
1011  auto t_w = getFTensor0IntegrationWeight();
1012  auto t_row_base = row_data.getFTensor0N();
1013  auto t_f = getFTensor0FromVec(*(commonDataPtr->plasticSurfacePtr));
1014  auto t_tau = getFTensor0FromVec(*(commonDataPtr->plasticTauPtr));
1015  auto t_tau_dot = getFTensor0FromVec(*(commonDataPtr->plasticTauDotPtr));
1016  auto t_flow =
1017  getFTensor2SymmetricFromMat<2>(*(commonDataPtr->plasticFlowPtr));
1018  auto &t_D = commonDataPtr->tD;
1019 
1020  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1021  double alpha = getMeasure() * t_w;
1022 
1023  auto mat_ptr = locMat.data().begin();
1024  auto t_diff_constrain_dstrain = diff_constrain_dstrain(
1025  t_D,
1027  diff_constrain_df(t_tau_dot, t_f - hardening(t_tau)), t_flow));
1028 
1030  &locMat(0, 0), &locMat(0, 1), &locMat(0, 2)};
1031 
1032  size_t rr = 0;
1033  for (; rr != nb_row_dofs; ++rr) {
1034 
1035  auto t_col_base = col_data.getFTensor0N(gg, 0);
1036  for (size_t cc = 0; cc != nb_col_dofs / 3; cc++) {
1037 
1038  t_mat(i, j) -=
1039  alpha * t_row_base * t_col_base * t_diff_constrain_dstrain(i, j);
1040 
1041  ++t_mat;
1042  ++t_col_base;
1043  }
1044 
1045  ++t_row_base;
1046  }
1047  for (; rr != nb_row_base_functions; ++rr)
1048  ++t_row_base;
1049 
1050  ++t_f;
1051  ++t_tau;
1052  ++t_tau_dot;
1053  ++t_flow;
1054  ++t_w;
1055  }
1056 
1057  CHKERR MatSetValues(getSNESB(), row_data, col_data, &*locMat.data().begin(),
1058  ADD_VALUES);
1059  }
1060 
1062 }
1063 
1065  const std::string row_field_name, const std::string col_field_name,
1066  boost::shared_ptr<CommonData> common_data_ptr)
1067  : DomainEleOp(row_field_name, col_field_name, DomainEleOp::OPROWCOL),
1068  commonDataPtr(common_data_ptr) {
1069  sYmm = false;
1070 }
1071 
1073  EntityType row_type,
1074  EntityType col_type,
1075  EntData &row_data,
1076  EntData &col_data) {
1078 
1079  const size_t nb_row_dofs = row_data.getIndices().size();
1080  const size_t nb_col_dofs = col_data.getIndices().size();
1081  if (nb_row_dofs && nb_col_dofs) {
1082 
1083  locMat.resize(nb_row_dofs, nb_col_dofs, false);
1084  locMat.clear();
1085 
1086  const size_t nb_integration_pts = row_data.getN().size1();
1087  const size_t nb_row_base_functions = row_data.getN().size2();
1088  auto t_w = getFTensor0IntegrationWeight();
1089  auto t_f = getFTensor0FromVec(*(commonDataPtr->plasticSurfacePtr));
1090  auto t_tau = getFTensor0FromVec(*(commonDataPtr->plasticTauPtr));
1091  auto t_tau_dot = getFTensor0FromVec(*(commonDataPtr->plasticTauDotPtr));
1092 
1093  auto t_row_base = row_data.getFTensor0N();
1094  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1095  const double alpha = getMeasure() * t_w;
1096  const double c0 = alpha * getTSa() *
1097  diff_constrain_dtau(t_tau_dot, t_f - hardening(t_tau));
1098  const double c1 = alpha *
1099  diff_constrain_df(t_tau_dot, t_f - hardening(t_tau)) *
1100  hardening_dtau();
1101 
1102  auto mat_ptr = locMat.data().begin();
1103 
1104  size_t rr = 0;
1105  for (; rr != nb_row_dofs; ++rr) {
1106 
1107  auto t_col_base = col_data.getFTensor0N(gg, 0);
1108  for (size_t cc = 0; cc != nb_col_dofs; ++cc) {
1109  *mat_ptr += (c0 - c1) * t_row_base * t_col_base;
1110  ++mat_ptr;
1111  ++t_col_base;
1112  }
1113  ++t_row_base;
1114  }
1115  for (; rr < nb_row_base_functions; ++rr)
1116  ++t_row_base;
1117 
1118  ++t_w;
1119  ++t_f;
1120  ++t_tau;
1121  ++t_tau_dot;
1122  }
1123 
1124  CHKERR MatSetValues(getSNESB(), row_data, col_data, &*locMat.data().begin(),
1125  ADD_VALUES);
1126  }
1127 
1129 }
1130 
1132  const std::string field_name, moab::Interface &post_proc_mesh,
1133  std::vector<EntityHandle> &map_gauss_pts,
1134  boost::shared_ptr<CommonData> common_data_ptr)
1135  : DomainEleOp(field_name, DomainEleOp::OPROW), postProcMesh(post_proc_mesh),
1136  mapGaussPts(map_gauss_pts), commonDataPtr(common_data_ptr) {
1137  // Opetor is only executed for vertices
1138  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1139 }
1140 
1141 //! [Postprocessing]
1143  EntData &data) {
1145 
1146  std::array<double, 9> def;
1147  std::fill(def.begin(), def.end(), 0);
1148 
1149  auto get_tag = [&](const std::string name, size_t size) {
1150  Tag th;
1151  CHKERR postProcMesh.tag_get_handle(name.c_str(), size, MB_TYPE_DOUBLE, th,
1152  MB_TAG_CREAT | MB_TAG_SPARSE,
1153  def.data());
1154  return th;
1155  };
1156 
1157  MatrixDouble3by3 mat(3, 3);
1158 
1159  auto set_matrix_2d = [&](auto &t) -> MatrixDouble3by3 & {
1160  mat.clear();
1161  for (size_t r = 0; r != 2; ++r)
1162  for (size_t c = 0; c != 2; ++c)
1163  mat(r, c) = t(r, c);
1164  return mat;
1165  };
1166 
1167  auto set_matrix_3d = [&](auto &t) -> MatrixDouble3by3 & {
1168  mat.clear();
1169  for (size_t r = 0; r != 3; ++r)
1170  for (size_t c = 0; c != 3; ++c)
1171  mat(r, c) = t(r, c);
1172  return mat;
1173  };
1174 
1175  auto set_matrix_2d_symm = [&](auto &t) -> MatrixDouble3by3 & {
1176  mat.clear();
1177  for (size_t r = 0; r != 2; ++r)
1178  for (size_t c = 0; c != 2; ++c)
1179  mat(r, c) = t(r, c);
1180  return mat;
1181  };
1182 
1183  auto set_scalar = [&](auto t) -> MatrixDouble3by3 & {
1184  mat.clear();
1185  mat(0, 0) = t;
1186  return mat;
1187  };
1188 
1189  auto set_tag = [&](auto th, auto gg, MatrixDouble3by3 &mat) {
1190  return postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
1191  &*mat.data().begin());
1192  };
1193 
1194  auto th_plastic_surface = get_tag("PLASTIC_SURFACE", 1);
1195  auto th_tau = get_tag("PLASTIC_MULTIPLIER", 1);
1196  auto th_plastic_flow = get_tag("PLASTIC_FLOW", 9);
1197  auto th_plastic_strain = get_tag("PLASTIC_STRAIN", 9);
1198 
1199  auto t_flow =
1200  getFTensor2SymmetricFromMat<2>(*(commonDataPtr->plasticFlowPtr));
1201  auto t_plastic_strain =
1202  getFTensor2SymmetricFromMat<2>(*(commonDataPtr->plasticStrainPtr));
1203  size_t gg = 0;
1204  for (int gg = 0; gg != commonDataPtr->plasticSurfacePtr->size(); ++gg) {
1205  const double f = (*(commonDataPtr->plasticSurfacePtr))[gg];
1206  const double tau = (*(commonDataPtr->plasticTauPtr))[gg];
1207  CHKERR set_tag(th_plastic_surface, gg, set_scalar(f - hardening(tau)));
1208  CHKERR set_tag(th_tau, gg, set_scalar(tau));
1209  CHKERR set_tag(th_plastic_flow, gg, set_matrix_2d(t_flow));
1210  CHKERR set_tag(th_plastic_strain, gg, set_matrix_2d(t_plastic_strain));
1211  ++t_flow;
1212  ++t_plastic_strain;
1213  }
1214 
1216 }
1217 
1218 struct Monitor : public FEMethod {
1219 
1220  Monitor(SmartPetscObj<DM> &dm,
1221  boost::shared_ptr<PostProcFaceOnRefinedMesh> &post_proc_fe,
1222  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> ux_scatter,
1223  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uy_scatter
1224 
1225  )
1226  : dM(dm), postProcFe(post_proc_fe), uXScatter(ux_scatter),
1227  uYScatter(uy_scatter){};
1228 
1229  MoFEMErrorCode preProcess() { return 0; }
1230  MoFEMErrorCode operator()() { return 0; }
1231 
1234 
1235  auto make_vtk = [&]() {
1238  CHKERR postProcFe->writeFile(
1239  "out_plastic_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
1241  };
1242 
1243  auto print_max_min = [&](auto &tuple, const std::string msg) {
1245  CHKERR VecScatterBegin(std::get<1>(tuple), ts_u, std::get<0>(tuple),
1246  INSERT_VALUES, SCATTER_FORWARD);
1247  CHKERR VecScatterEnd(std::get<1>(tuple), ts_u, std::get<0>(tuple),
1248  INSERT_VALUES, SCATTER_FORWARD);
1249  double max, min;
1250  CHKERR VecMax(std::get<0>(tuple), PETSC_NULL, &max);
1251  CHKERR VecMin(std::get<0>(tuple), PETSC_NULL, &min);
1252  PetscPrintf(PETSC_COMM_WORLD, "%s time %3.4e min %3.4e max %3.4e\n",
1253  msg.c_str(), ts_t, min, max);
1255  };
1256 
1257  CHKERR make_vtk();
1258  CHKERR print_max_min(uXScatter, "Ux");
1259  CHKERR print_max_min(uYScatter, "Uy");
1260 
1262  }
1263 
1264 private:
1265  SmartPetscObj<DM> dM;
1266  boost::shared_ptr<PostProcFaceOnRefinedMesh> postProcFe;
1267  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
1268  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
1269 };
1270 
1271 //! [Postprocessing]
1272 }; // 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
boost::shared_ptr< PostProcFaceOnRefinedMesh > postProcFe
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Definition: PlasticOps.hpp:758
FTensor::Index< 'k', 2 > k
Definition: PlasticOps.hpp:65
MoFEMErrorCode postProcess()
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
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:186
static FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:141
FTensor::Index< 'j', 2 > j
Definition: PlasticOps.hpp:64
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
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: PlasticOps.hpp:411
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:134
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:173
double diff_constrain_dtau(double tau, double f)
Definition: PlasticOps.hpp:375
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
FTensor::Index< 'l', 2 > l
Definition: PlasticOps.hpp:66
boost::shared_ptr< MatrixDouble > plasticStrainDotPtr
Definition: PlasticOps.hpp:59
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:147
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:82
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Postprocessing]
auto hardening(double tau)
Definition: PlasticOps.hpp:281
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:75
FTensor::Index< 'n', 2 > n
Definition: PlasticOps.hpp:68
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Definition: PlasticOps.hpp:583
boost::shared_ptr< VectorDouble > plasticSurfacePtr
Definition: PlasticOps.hpp:54
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
FTensor::Index< 'I', 3 > I
Definition: PlasticOps.hpp:70
FTensor::Tensor1< double *, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
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
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:100
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:199
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:213
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:507
boost::shared_ptr< VectorDouble > plasticTauDotPtr
Definition: PlasticOps.hpp:57
FTensor::Index< 'J', 3 > J
Definition: PlasticOps.hpp:71
auto diff_constrain_dstress(double &&diff_constrain_df, FTensor::Tensor2_symmetric< T, 2 > &t_plastic_flow)
Definition: PlasticOps.hpp:385
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 preProcess()
double platsic_surface(FTensor::Tensor2_symmetric< double, 3 > &&t_stress_deviator)
Definition: PlasticOps.hpp:309
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
moab::Interface & postProcMesh
Definition: PlasticOps.hpp:211
bool sYmm
If true assume that matrix is symmetric structure.
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
OpCalculatePlasticSurface(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
[Lambda functions]
Definition: PlasticOps.hpp:403
Definition: PlasticOps.hpp:203
auto diff_tensor()
[Operators definitions]
Definition: PlasticOps.hpp:218
std::vector< EntityHandle > & mapGaussPts
Definition: PlasticOps.hpp:212
FTensor::Index< 'm', 2 > m
Definition: PlasticOps.hpp:67
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
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
boost::shared_ptr< MatrixDouble > plasticFlowPtr
Definition: PlasticOps.hpp:55
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
auto getFTensor0IntegrationWeight()
Get integration weights.
MoFEMErrorCode MatSetValues(Mat M, const DataForcesAndSourcesCore::EntData &row_data, const DataForcesAndSourcesCore::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
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
boost::shared_ptr< VectorDouble > plasticTauPtr
Definition: PlasticOps.hpp:56
auto diff_constrain_dstrain(FTensor::Ddg< T, 2, 2 > &t_D, FTensor::Tensor2_symmetric< T, 2 > &&t_diff_constrain_dstress)
Definition: PlasticOps.hpp:393
constexpr double cn
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:101
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:160
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.
MoFEMErrorCode VecSetValues(Vec V, const DataForcesAndSourcesCore::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
#define CHKERR
Inline error check.
Definition: definitions.h:602
FTensor::Index< 'M', 3 > M
Definition: PlasticOps.hpp:72
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
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
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
FTensor::Index< 'i', 2 > i
[Common data]
Definition: PlasticOps.hpp:63
SmartPetscObj< DM > dM
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1788
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:413
auto diff_symmetrize()
Definition: PlasticOps.hpp:227
constexpr double third
Monitor(SmartPetscObj< DM > &dm, boost::shared_ptr< PostProcFaceOnRefinedMesh > &post_proc_fe, std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter >> ux_scatter, std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter >> uy_scatter)
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOps.hpp:91
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
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< MatrixDouble > plasticStrainPtr
Definition: PlasticOps.hpp:58
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< 'N', 3 > N
Definition: PlasticOps.hpp:73
MoFEMErrorCode operator()()
auto diff_constrain_df(double tau, double f)
Definition: PlasticOps.hpp:379