v0.14.0
PlasticOpsGeneric.hpp
Go to the documentation of this file.
1 
2 
3 /** \file PlasticOpsGeneric.hpp
4  * \example PlasticOpsGeneric.hpp
5  */
6 
7 namespace PlasticOps {
8 
9 //! [Lambda functions]
10 template <int DIM> inline auto diff_tensor(FTensor::Number<DIM>) {
17  t_diff(i, j, k, l) = (t_kd(i, k) ^ t_kd(j, l)) / 4.;
18  return t_diff;
19 };
20 
21 template <int DIM> inline auto symm_L_tensor(FTensor::Number<DIM>) {
24  constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
27  t_L(i, j, L) = 0;
28  if constexpr (DIM == 2) {
29  t_L(0, 0, 0) = 1;
30  t_L(1, 0, 1) = 1;
31  t_L(1, 1, 2) = 1;
32  } else {
33  t_L(0, 0, 0) = 1;
34  t_L(1, 0, 1) = 1;
35  t_L(2, 0, 2) = 1;
36  t_L(1, 1, 3) = 1;
37  t_L(2, 1, 4) = 1;
38  t_L(2, 2, 5) = 1;
39  }
40  return t_L;
41 }
42 
43 template <int DIM> inline auto diff_symmetrize(FTensor::Number<DIM>) {
44 
49 
51 
52  t_diff(i, j, k, l) = 0;
53  t_diff(0, 0, 0, 0) = 1;
54  t_diff(1, 1, 1, 1) = 1;
55 
56  t_diff(1, 0, 1, 0) = 0.5;
57  t_diff(1, 0, 0, 1) = 0.5;
58 
59  t_diff(0, 1, 0, 1) = 0.5;
60  t_diff(0, 1, 1, 0) = 0.5;
61 
62  if constexpr (DIM == 3) {
63  t_diff(2, 2, 2, 2) = 1;
64 
65  t_diff(2, 0, 2, 0) = 0.5;
66  t_diff(2, 0, 0, 2) = 0.5;
67  t_diff(0, 2, 0, 2) = 0.5;
68  t_diff(0, 2, 2, 0) = 0.5;
69 
70  t_diff(2, 1, 2, 1) = 0.5;
71  t_diff(2, 1, 1, 2) = 0.5;
72  t_diff(1, 2, 1, 2) = 0.5;
73  t_diff(1, 2, 2, 1) = 0.5;
74  }
75 
76  return t_diff;
77 };
78 
79 template <typename T>
80 inline double trace(FTensor::Tensor2_symmetric<T, 2> &t_stress) {
81  constexpr double third = boost::math::constants::third<double>();
82  return (t_stress(0, 0) + t_stress(1, 1)) * third;
83 };
84 
85 template <typename T>
86 inline double trace(FTensor::Tensor2_symmetric<T, 3> &t_stress) {
87  constexpr double third = boost::math::constants::third<double>();
88  return (t_stress(0, 0) + t_stress(1, 1) + t_stress(2, 2)) * third;
89 };
90 
91 template <typename T, int DIM>
92 inline auto deviator(
93 
94  FTensor::Tensor2_symmetric<T, DIM> &t_stress, double trace,
96 
97 ) {
99  t_dev(I, J) = 0;
100  for (int ii = 0; ii != DIM; ++ii)
101  for (int jj = ii; jj != DIM; ++jj)
102  t_dev(ii, jj) = t_stress(ii, jj);
103  t_dev(0, 0) -= trace;
104  t_dev(1, 1) -= trace;
105  t_dev(2, 2) -= trace;
106  for (int ii = 0; ii != DIM; ++ii)
107  for (int jj = ii; jj != DIM; ++jj)
108  t_dev(ii, jj) -= t_alpha(ii, jj);
109  return t_dev;
110 };
111 
112 template <typename T>
113 inline auto deviator(FTensor::Tensor2_symmetric<T, 2> &t_stress, double trace,
115  return deviator(t_stress, trace, t_alpha, FTensor::Number<2>());
116 };
117 
118 template <typename T>
119 inline auto deviator(FTensor::Tensor2_symmetric<T, 3> &t_stress, double trace,
121  return deviator(t_stress, trace, t_alpha, FTensor::Number<3>());
122 };
123 
124 template <int DIM>
125 inline auto diff_deviator(FTensor::Ddg<double, DIM, DIM> &&t_diff_stress,
129  FTensor::Ddg<double, 3, DIM> t_diff_deviator;
130  t_diff_deviator(I, J, k, l) = 0;
131  for (int ii = 0; ii != DIM; ++ii)
132  for (int jj = ii; jj != DIM; ++jj)
133  for (int kk = 0; kk != DIM; ++kk)
134  for (int ll = kk; ll != DIM; ++ll)
135  t_diff_deviator(ii, jj, kk, ll) = t_diff_stress(ii, jj, kk, ll);
136 
137  constexpr double third = boost::math::constants::third<double>();
138 
139  t_diff_deviator(0, 0, 0, 0) -= third;
140  t_diff_deviator(0, 0, 1, 1) -= third;
141 
142  t_diff_deviator(1, 1, 0, 0) -= third;
143  t_diff_deviator(1, 1, 1, 1) -= third;
144 
145  t_diff_deviator(2, 2, 0, 0) -= third;
146  t_diff_deviator(2, 2, 1, 1) -= third;
147 
148  if constexpr (DIM == 3) {
149  t_diff_deviator(0, 0, 2, 2) -= third;
150  t_diff_deviator(1, 1, 2, 2) -= third;
151  t_diff_deviator(2, 2, 2, 2) -= third;
152  }
153 
154  return t_diff_deviator;
155 };
156 
157 inline auto diff_deviator(FTensor::Ddg<double, 2, 2> &&t_diff_stress) {
158  return diff_deviator(std::move(t_diff_stress), FTensor::Number<2>());
159 }
160 
161 inline auto diff_deviator(FTensor::Ddg<double, 3, 3> &&t_diff_stress) {
162  return diff_deviator(std::move(t_diff_stress), FTensor::Number<3>());
163 }
164 
165 /**
166  *
167 
168 \f[
169 \begin{split}
170 f&=\sqrt{s_{ij}s_{ij}}\\
171 A_{ij}&=\frac{\partial f}{\partial \sigma_{ij}}=
172 \frac{1}{f} s_{kl} \frac{\partial s_{kl}}{\partial \sigma_{ij}}\\
173 \frac{\partial A_{ij}}{\partial \sigma_{kl}}&= \frac{\partial^2 f}{\partial
174 \sigma_{ij}\partial\sigma_{mn}}= \frac{1}{f} \left( \frac{\partial
175 s_{kl}}{\partial \sigma_{mn}}\frac{\partial s_{kl}}{\partial \sigma_{ij}}
176 -A_{mn}A_{ij}
177 \right)\\
178 \frac{\partial f}{\partial \varepsilon_{ij}}&=A_{mn}D_{mnij}
179 \\
180 \frac{\partial A_{ij}}{\partial \varepsilon_{kl}}&=
181 \frac{\partial A_{ij}}{\partial \sigma_{mn}} \frac{\partial
182 \sigma_{mn}}{\partial \varepsilon_{kl}}= \frac{\partial A_{ij}}{\partial
183 \sigma_{mn}} D_{mnkl}
184 \end{split}
185 \f]
186 
187  */
188 inline double
190  return std::sqrt(1.5 * t_stress_deviator(I, J) * t_stress_deviator(I, J)) +
191  std::numeric_limits<double>::epsilon();
192 };
193 
194 template <int DIM>
195 inline auto plastic_flow(long double f,
197  FTensor::Ddg<double, 3, DIM> &&t_diff_deviator) {
201  t_diff_f(k, l) =
202  (1.5 * (t_dev_stress(I, J) * t_diff_deviator(I, J, k, l))) / f;
203  return t_diff_f;
204 };
205 
206 template <typename T, int DIM>
207 inline auto
210  FTensor::Ddg<double, 3, DIM> &&t_diff_deviator) {
215  FTensor::Ddg<double, DIM, DIM> t_diff_flow;
216  t_diff_flow(i, j, k, l) =
217  (1.5 * (t_diff_deviator(M, N, i, j) * t_diff_deviator(M, N, k, l) -
218  (2. / 3.) * t_flow(i, j) * t_flow(k, l))) /
219  f;
220  return t_diff_flow;
221 };
222 
223 template <typename T, int DIM>
226  FTensor::Ddg<double, DIM, DIM> &&t_diff_plastic_flow_dstress) {
233  FTensor::Ddg<double, DIM, DIM> t_diff_flow;
234  t_diff_flow(i, j, k, l) =
235  t_diff_plastic_flow_dstress(i, j, m, n) * t_D(m, n, k, l);
236  return t_diff_flow;
237 };
238 
239 // inline double constrain_diff_sign(double x) {
240 // const auto y = x / zeta;
241 // if (y > std::numeric_limits<float>::max_exponent10 ||
242 // y < std::numeric_limits<float>::min_exponent10) {
243 // return 0;
244 // } else {
245 // const auto e = std::exp(y);
246 // const auto ep1 = e + 1;
247 // return (2 / zeta) * (e / (ep1 * ep1));
248 // }
249 // };
250 
251 inline double constrian_sign(double x) {
252  const auto y = x / zeta;
253  if (y > std::numeric_limits<float>::max_exponent10 ||
254  y < std::numeric_limits<float>::min_exponent10) {
255  if (x > 0)
256  return 1.;
257  else
258  return -1.;
259  } else {
260  const auto e = std::exp(y);
261  return (e - 1) / (1 + e);
262  }
263 };
264 
265 inline double constrain_abs(double x) {
266  const auto y = -x / zeta;
267  if (y > std::numeric_limits<float>::max_exponent10 ||
268  y < std::numeric_limits<float>::min_exponent10) {
269  return std::abs(x);
270  } else {
271  const double e = std::exp(y);
272  return x + 2 * zeta * std::log1p(e);
273  }
274 };
275 
276 inline double w(double eqiv, double dot_tau, double f, double sigma_y,
277  double sigma_Y) {
278  return (f - sigma_y) / sigma_Y + cn1 * (dot_tau);
279 };
280 
281 /**
282 
283 \f[
284 \dot{\tau} - \frac{1}{2}\left\{\dot{\tau} + (f(\pmb\sigma) - \sigma_y) +
285 \| \dot{\tau} + (f(\pmb\sigma) - \sigma_y) \|\right\} = 0 \\
286 c_n \sigma_y \dot{\tau} - \frac{1}{2}\left\{c_n\sigma_y \dot{\tau} +
287 (f(\pmb\sigma) - \sigma_y) +
288 \| c_n \sigma_y \dot{\tau} + (f(\pmb\sigma) - \sigma_y) \|\right\} = 0
289 \f]
290 
291  */
292 inline double constraint(double eqiv, double dot_tau, double f, double sigma_y,
293  double abs_w, double vis_H, double sigma_Y) {
294  return vis_H * dot_tau +
295  (sigma_Y / 2) * ((cn0 * (dot_tau - eqiv) + cn1 * (dot_tau) -
296  (f - sigma_y) / sigma_Y) -
297  abs_w);
298 };
299 
300 inline double diff_constrain_ddot_tau(double sign, double eqiv, double dot_tau,
301  double vis_H, double sigma_Y) {
302  return vis_H + (sigma_Y / 2) * (cn0 + cn1 * (1 - sign));
303 };
304 
305 inline double diff_constrain_deqiv(double sign, double eqiv, double dot_tau,
306  double sigma_Y) {
307  return (sigma_Y / 2) * (-cn0);
308 };
309 
310 inline auto diff_constrain_df(double sign) { return (-1 - sign) / 2; };
311 
312 inline auto diff_constrain_dsigma_y(double sign) { return (1 + sign) / 2; }
313 
314 template <typename T, int DIM>
315 inline auto
317  FTensor::Tensor2_symmetric<T, DIM> &t_plastic_flow) {
320  FTensor::Tensor2_symmetric<double, DIM> t_diff_constrain_dstress;
321  t_diff_constrain_dstress(i, j) = diff_constrain_df * t_plastic_flow(i, j);
322  return t_diff_constrain_dstress;
323 };
324 
325 template <typename T1, typename T2, int DIM>
326 inline auto diff_constrain_dstrain(T1 &t_D, T2 &&t_diff_constrain_dstress) {
331  FTensor::Tensor2_symmetric<double, DIM> t_diff_constrain_dstrain;
332  t_diff_constrain_dstrain(k, l) =
333  t_diff_constrain_dstress(i, j) * t_D(i, j, k, l);
334  return t_diff_constrain_dstrain;
335 };
336 
337 template <typename T, int DIM>
339  FTensor::Tensor2_symmetric<T, DIM> &t_plastic_strain_dot) {
342  constexpr double A = 2. / 3;
343  return std::sqrt(A * t_plastic_strain_dot(i, j) *
344  t_plastic_strain_dot(i, j)) +
345  std::numeric_limits<double>::epsilon();
346 };
347 
348 template <typename T1, typename T2, typename T3, int DIM>
349 inline auto diff_equivalent_strain_dot(const T1 eqiv, T2 &t_plastic_strain_dot,
350  T3 &t_diff_plastic_strain,
356  constexpr double A = 2. / 3;
358  t_diff_eqiv(i, j) = A * (t_plastic_strain_dot(k, l) / eqiv) *
359  t_diff_plastic_strain(k, l, i, j);
360  return t_diff_eqiv;
361 };
362 
363 //! [Lambda functions]
364 
365 //! [Auxiliary functions functions
366 static inline auto get_mat_tensor_sym_dvector(size_t rr, MatrixDouble &mat,
369  &mat(3 * rr + 0, 0), &mat(3 * rr + 0, 1), &mat(3 * rr + 1, 0),
370  &mat(3 * rr + 1, 1), &mat(3 * rr + 2, 0), &mat(3 * rr + 2, 1)};
371 }
372 
373 static inline auto get_mat_tensor_sym_dvector(size_t rr, MatrixDouble &mat,
376  &mat(6 * rr + 0, 0), &mat(6 * rr + 0, 1), &mat(6 * rr + 0, 2),
377  &mat(6 * rr + 1, 0), &mat(6 * rr + 1, 1), &mat(6 * rr + 1, 2),
378  &mat(6 * rr + 2, 0), &mat(6 * rr + 2, 1), &mat(6 * rr + 2, 2),
379  &mat(6 * rr + 3, 0), &mat(6 * rr + 3, 1), &mat(6 * rr + 3, 2),
380  &mat(6 * rr + 4, 0), &mat(6 * rr + 4, 1), &mat(6 * rr + 4, 2),
381  &mat(6 * rr + 5, 0), &mat(6 * rr + 5, 1), &mat(6 * rr + 5, 2)};
382 }
383 
384 //! [Auxiliary functions functions
385 
386 template <int DIM, typename DomainEleOp>
388  : public DomainEleOp {
389  OpCalculatePlasticSurfaceImpl(const std::string field_name,
390  boost::shared_ptr<CommonData> common_data_ptr);
391  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
392 
393 protected:
394  boost::shared_ptr<CommonData> commonDataPtr;
395 };
396 
397 template <int DIM, typename DomainEleOp>
400  boost::shared_ptr<CommonData> common_data_ptr)
402  commonDataPtr(common_data_ptr) {
403  // Operator is only executed for vertices
404  std::fill(&DomainEleOp::doEntities[MBEDGE],
405  &DomainEleOp::doEntities[MBMAXTYPE], false);
406 }
407 
408 template <int DIM, typename DomainEleOp>
410  int side, EntityType type, EntData &data) {
412 
415 
416  const size_t nb_gauss_pts = commonDataPtr->mStressPtr->size2();
417  auto t_stress =
418  getFTensor2SymmetricFromMat<DIM>(*(commonDataPtr->mStressPtr));
419  auto t_plastic_strain =
420  getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticStrain);
421 
422  commonDataPtr->plasticSurface.resize(nb_gauss_pts, false);
423  commonDataPtr->plasticFlow.resize(size_symm, nb_gauss_pts, false);
424  auto t_flow = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticFlow);
425 
426  auto &params = commonDataPtr->blockParams;
427 
428  for (auto &f : commonDataPtr->plasticSurface) {
429 
430  f = platsic_surface(
431 
432  deviator(
433  t_stress, trace(t_stress),
434  kinematic_hardening(t_plastic_strain, params[CommonData::C1_k]))
435 
436  );
437 
438  auto t_flow_tmp =
439  plastic_flow(f,
440 
441  deviator(t_stress, trace(t_stress),
442  kinematic_hardening(t_plastic_strain,
443  params[CommonData::C1_k])),
444 
446  t_flow(i, j) = t_flow_tmp(i, j);
447 
448  ++t_plastic_strain;
449  ++t_flow;
450  ++t_stress;
451  }
452 
454 }
455 
456 template <int DIM, typename DomainEleOp>
458  OpCalculatePlasticityImpl(const std::string field_name,
459  boost::shared_ptr<CommonData> common_data_ptr,
460  boost::shared_ptr<MatrixDouble> m_D_ptr);
461  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
462 
463 protected:
464  boost::shared_ptr<CommonData> commonDataPtr;
465  boost::shared_ptr<MatrixDouble> mDPtr;
466 };
467 
468 template <int DIM, typename DomainEleOp>
470  const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr,
471  boost::shared_ptr<MatrixDouble> m_D_ptr)
473  commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
474  // Opetor is only executed for vertices
475  std::fill(&DomainEleOp::doEntities[MBEDGE],
476  &DomainEleOp::doEntities[MBMAXTYPE], false);
477 }
478 
479 template <int DIM, typename DomainEleOp>
481  int side, EntityType type, EntData &data) {
483 
490 
491  auto &params = commonDataPtr->blockParams; ///< material parameters
492 
493  const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
494  auto t_w = DomainEleOp::getFTensor0IntegrationWeight();
495  auto t_tau = getFTensor0FromVec(commonDataPtr->plasticTau);
496  auto t_tau_dot = getFTensor0FromVec(commonDataPtr->plasticTauDot);
497  auto t_f = getFTensor0FromVec(commonDataPtr->plasticSurface);
498  auto t_flow = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticFlow);
499  auto t_plastic_strain =
500  getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticStrain);
501  auto t_plastic_strain_dot =
502  getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticStrainDot);
503  auto t_stress =
504  getFTensor2SymmetricFromMat<DIM>(*(commonDataPtr->mStressPtr));
505 
506  auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*commonDataPtr->mDPtr);
507  auto t_D_Op = getFTensor4DdgFromMat<DIM, DIM, 0>(*mDPtr);
508 
509  auto t_diff_plastic_strain = diff_tensor(FTensor::Number<DIM>());
510  auto t_diff_deviator = diff_deviator(diff_tensor(FTensor::Number<DIM>()));
511 
512  FTensor::Ddg<double, DIM, DIM> t_flow_dir_dstress;
513  FTensor::Ddg<double, DIM, DIM> t_flow_dir_dstrain;
514  t_flow_dir_dstress(i, j, k, l) =
515  1.5 * (t_diff_deviator(M, N, i, j) * t_diff_deviator(M, N, k, l));
516  t_flow_dir_dstrain(i, j, k, l) =
517  t_flow_dir_dstress(i, j, m, n) * t_D_Op(m, n, k, l);
518 
519 
520  auto t_alpha_dir =
521  kinematic_hardening_dplastic_strain<DIM>(params[CommonData::C1_k]);
522 
523  commonDataPtr->resC.resize(nb_gauss_pts, false);
524  commonDataPtr->resCdTau.resize(nb_gauss_pts, false);
525  commonDataPtr->resCdStrain.resize(size_symm, nb_gauss_pts, false);
526  commonDataPtr->resCdPlasticStrain.resize(size_symm, nb_gauss_pts, false);
527  commonDataPtr->resFlow.resize(size_symm, nb_gauss_pts, false);
528  commonDataPtr->resFlowDtau.resize(size_symm, nb_gauss_pts, false);
529  commonDataPtr->resFlowDstrain.resize(size_symm * size_symm, nb_gauss_pts,
530  false);
531  commonDataPtr->resFlowDstrainDot.resize(size_symm * size_symm, nb_gauss_pts,
532  false);
533 
534  commonDataPtr->resC.clear();
535  commonDataPtr->resCdTau.clear();
536  commonDataPtr->resCdStrain.clear();
537  commonDataPtr->resCdPlasticStrain.clear();
538  commonDataPtr->resFlow.clear();
539  commonDataPtr->resFlowDtau.clear();
540  commonDataPtr->resFlowDstrain.clear();
541  commonDataPtr->resFlowDstrainDot.clear();
542 
543  auto t_res_c = getFTensor0FromVec(commonDataPtr->resC);
544  auto t_res_c_dtau = getFTensor0FromVec(commonDataPtr->resCdTau);
545  auto t_res_c_dstrain =
546  getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resCdStrain);
547  auto t_res_c_plastic_strain =
548  getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resCdPlasticStrain);
549  auto t_res_flow = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resFlow);
550  auto t_res_flow_dtau =
551  getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resFlowDtau);
552  auto t_res_flow_dstrain =
553  getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->resFlowDstrain);
554  auto t_res_flow_dplastic_strain =
555  getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->resFlowDstrainDot);
556 
557  auto next = [&]() {
558  ++t_tau;
559  ++t_tau_dot;
560  ++t_f;
561  ++t_flow;
562  ++t_plastic_strain;
563  ++t_plastic_strain_dot;
564  ++t_stress;
565  ++t_res_c;
566  ++t_res_c_dtau;
567  ++t_res_c_dstrain;
568  ++t_res_c_plastic_strain;
569  ++t_res_flow;
570  ++t_res_flow_dtau;
571  ++t_res_flow_dstrain;
572  ++t_res_flow_dplastic_strain;
573  ++t_w;
574  };
575 
576  auto get_avtive_pts = [&]() {
577  int nb_points_avtive_on_elem = 0;
578  int nb_points_on_elem = 0;
579 
580  auto t_tau = getFTensor0FromVec(commonDataPtr->plasticTau);
581  auto t_tau_dot = getFTensor0FromVec(commonDataPtr->plasticTauDot);
582  auto t_f = getFTensor0FromVec(commonDataPtr->plasticSurface);
583  auto t_plastic_strain_dot =
584  getFTensor2SymmetricFromMat<SPACE_DIM>(commonDataPtr->plasticStrainDot);
585 
586  for (auto &f : commonDataPtr->plasticSurface) {
587  auto eqiv = equivalent_strain_dot(t_plastic_strain_dot);
588  const auto ww = w(
589  eqiv, t_tau_dot, t_f,
590  iso_hardening(t_tau, params[CommonData::H], params[CommonData::QINF],
591  params[CommonData::BISO], params[CommonData::SIGMA_Y]),
592  params[CommonData::SIGMA_Y]);
593  const auto sign_ww = constrian_sign(ww);
594 
595  ++nb_points_on_elem;
596  if (sign_ww > 0) {
597  ++nb_points_avtive_on_elem;
598  }
599 
600  ++t_tau;
601  ++t_tau_dot;
602  ++t_f;
603  ++t_plastic_strain_dot;
604  }
605 
606  int &active_points = PlasticOps::CommonData::activityData[0];
607  int &avtive_full_elems = PlasticOps::CommonData::activityData[1];
608  int &avtive_elems = PlasticOps::CommonData::activityData[2];
609  int &nb_points = PlasticOps::CommonData::activityData[3];
610  int &nb_elements = PlasticOps::CommonData::activityData[4];
611 
612  ++nb_elements;
613  nb_points += nb_points_on_elem;
614  if (nb_points_avtive_on_elem > 0) {
615  ++avtive_elems;
616  active_points += nb_points_avtive_on_elem;
617  if (nb_points_avtive_on_elem == nb_points_on_elem) {
618  ++avtive_full_elems;
619  }
620  }
621 
622  if (nb_points_avtive_on_elem != nb_points_on_elem)
623  return 1;
624  else
625  return 0;
626  };
627 
628  if (DomainEleOp::getTSCtx() == TSMethod::TSContext::CTX_TSSETIJACOBIAN) {
629  get_avtive_pts();
630  }
631 
632  for (auto &f : commonDataPtr->plasticSurface) {
633 
634  auto eqiv = equivalent_strain_dot(t_plastic_strain_dot);
635  auto t_diff_eqiv = diff_equivalent_strain_dot(eqiv, t_plastic_strain_dot,
636  t_diff_plastic_strain,
638 
639  const auto sigma_y =
640  iso_hardening(t_tau, params[CommonData::H], params[CommonData::QINF],
641  params[CommonData::BISO], params[CommonData::SIGMA_Y]);
642  const auto d_sigma_y =
643  iso_hardening_dtau(t_tau, params[CommonData::H],
644  params[CommonData::QINF], params[CommonData::BISO]);
645 
646  auto ww = w(eqiv, t_tau_dot, t_f, sigma_y, params[CommonData::SIGMA_Y]);
647  auto abs_ww = constrain_abs(ww);
648  auto sign_ww = constrian_sign(ww);
649 
650  auto c = constraint(eqiv, t_tau_dot, t_f, sigma_y, abs_ww,
651  params[CommonData::VIS_H], params[CommonData::SIGMA_Y]);
652  auto c_dot_tau = diff_constrain_ddot_tau(sign_ww, eqiv, t_tau_dot,
653  params[CommonData::VIS_H],
654  params[CommonData::SIGMA_Y]);
655  auto c_equiv = diff_constrain_deqiv(sign_ww, eqiv, t_tau_dot,
656  params[CommonData::SIGMA_Y]);
657  auto c_sigma_y = diff_constrain_dsigma_y(sign_ww);
658  auto c_f = diff_constrain_df(sign_ww);
659 
660  auto t_dev_stress = deviator(
661 
662  t_stress, trace(t_stress),
663 
664  kinematic_hardening(t_plastic_strain, params[CommonData::C1_k])
665 
666  );
667 
669  t_flow_dir(k, l) = 1.5 * (t_dev_stress(I, J) * t_diff_deviator(I, J, k, l));
671  t_flow_dstrain(i, j) = t_flow(k, l) * t_D_Op(k, l, i, j);
672 
673  auto get_res_c = [&]() { return c; };
674 
675  auto get_res_c_dstrain = [&](auto &t_diff_res) {
676  t_diff_res(i, j) = c_f * t_flow_dstrain(i, j);
677  };
678 
679  auto get_res_c_dplastic_strain = [&](auto &t_diff_res) {
680  t_diff_res(i, j) = (this->getTSa() * c_equiv) * t_diff_eqiv(i, j);
681  t_diff_res(k, l) -= c_f * t_flow(i, j) * t_alpha_dir(i, j, k, l);
682  };
683 
684  auto get_res_c_dtau = [&]() {
685  return this->getTSa() * c_dot_tau + c_sigma_y * d_sigma_y;
686  };
687 
688  auto get_res_c_plastic_strain = [&](auto &t_diff_res) {
689  t_diff_res(k, l) = -c_f * t_flow(i, j) * t_alpha_dir(i, j, k, l);
690  };
691 
692  auto get_res_flow = [&](auto &t_res_flow) {
693  const auto a = sigma_y;
694  const auto b = t_tau_dot;
695  t_res_flow(k, l) = a * t_plastic_strain_dot(k, l) - b * t_flow_dir(k, l);
696  };
697 
698  auto get_res_flow_dtau = [&](auto &t_res_flow_dtau) {
699  const auto da = d_sigma_y;
700  const auto db = this->getTSa();
701  t_res_flow_dtau(k, l) =
702  da * t_plastic_strain_dot(k, l) - db * t_flow_dir(k, l);
703  };
704 
705  auto get_res_flow_dstrain = [&](auto &t_res_flow_dstrain) {
706  const auto b = t_tau_dot;
707  t_res_flow_dstrain(m, n, k, l) = -t_flow_dir_dstrain(m, n, k, l) * b;
708  };
709 
710  auto get_res_flow_dplastic_strain = [&](auto &t_res_flow_dplastic_strain) {
711  const auto a = sigma_y;
712  t_res_flow_dplastic_strain(m, n, k, l) =
713  (a * this->getTSa()) * t_diff_plastic_strain(m, n, k, l);
714  const auto b = t_tau_dot;
715  t_res_flow_dplastic_strain(m, n, i, j) +=
716  (t_flow_dir_dstrain(m, n, k, l) * t_alpha_dir(k, l, i, j)) * b;
717  };
718 
719  t_res_c = get_res_c();
720  get_res_flow(t_res_flow);
721 
722  if (this->getTSCtx() == TSMethod::TSContext::CTX_TSSETIJACOBIAN) {
723  t_res_c_dtau = get_res_c_dtau();
724  get_res_c_dstrain(t_res_c_dstrain);
725  get_res_c_dplastic_strain(t_res_c_plastic_strain);
726  get_res_flow_dtau(t_res_flow_dtau);
727  get_res_flow_dstrain(t_res_flow_dstrain);
728  get_res_flow_dplastic_strain(t_res_flow_dplastic_strain);
729  }
730 
731  next();
732  }
733 
735 }
736 
737 template <int DIM, typename DomainEleOp>
739 
740  /**
741  * @deprecated do not use this constructor
742  */
743  DEPRECATED OpPlasticStressImpl(const std::string field_name,
744  boost::shared_ptr<CommonData> common_data_ptr,
745  boost::shared_ptr<MatrixDouble> mDPtr);
746  OpPlasticStressImpl(boost::shared_ptr<CommonData> common_data_ptr,
747  boost::shared_ptr<MatrixDouble> mDPtr);
748 
749  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
750 
751 private:
752  boost::shared_ptr<MatrixDouble> mDPtr;
753  boost::shared_ptr<CommonData> commonDataPtr;
754 };
755 
756 template <int DIM, typename DomainEleOp>
758  const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr,
759  boost::shared_ptr<MatrixDouble> m_D_ptr)
761  commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
762  // Operator is only executed for vertices
763  std::fill(&DomainEleOp::doEntities[MBEDGE],
764  &DomainEleOp::doEntities[MBMAXTYPE], false);
765 }
766 
767 template <int DIM, typename DomainEleOp>
769  boost::shared_ptr<CommonData> common_data_ptr,
770  boost::shared_ptr<MatrixDouble> m_D_ptr)
771  : DomainEleOp(NOSPACE, DomainEleOp::OPSPACE),
772  commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {}
773 
774 //! [Calculate stress]
775 template <int DIM, typename DomainEleOp>
778  EntData &data) {
780 
785 
786  const size_t nb_gauss_pts = commonDataPtr->mStrainPtr->size2();
787  commonDataPtr->mStressPtr->resize((DIM * (DIM + 1)) / 2, nb_gauss_pts);
788  auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*mDPtr);
789  auto t_strain =
790  getFTensor2SymmetricFromMat<DIM>(*(commonDataPtr->mStrainPtr));
791  auto t_plastic_strain =
792  getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticStrain);
793  auto t_stress =
794  getFTensor2SymmetricFromMat<DIM>(*(commonDataPtr->mStressPtr));
795 
796  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
797  t_stress(i, j) =
798  t_D(i, j, k, l) * (t_strain(k, l) - t_plastic_strain(k, l));
799  ++t_strain;
800  ++t_plastic_strain;
801  ++t_stress;
802  }
803 
805 }
806 //! [Calculate stress]
807 
808 template <int DIM, typename AssemblyDomainEleOp>
810  : public AssemblyDomainEleOp {
811  OpCalculatePlasticFlowRhsImpl(const std::string field_name,
812  boost::shared_ptr<CommonData> common_data_ptr,
813  boost::shared_ptr<MatrixDouble> m_D_ptr);
815 
816 private:
817  boost::shared_ptr<CommonData> commonDataPtr;
818  boost::shared_ptr<MatrixDouble> mDPtr;
819 };
820 
821 template <int DIM, typename AssemblyDomainEleOp>
824  boost::shared_ptr<CommonData> common_data_ptr,
825  boost::shared_ptr<MatrixDouble> m_D_ptr)
827  commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {}
828 
829 template <int DIM, typename AssemblyDomainEleOp>
834 
839  constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
841 
842  const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
843  const auto nb_base_functions = data.getN().size2();
844 
845  auto t_res_flow = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resFlow);
846 
847  auto t_L = symm_L_tensor(FTensor::Number<DIM>());
848 
849  auto next = [&]() { ++t_res_flow; };
850 
851  auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
852  auto t_base = data.getFTensor0N();
853  auto &nf = AssemblyDomainEleOp::locF;
854  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
855  double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
856  ++t_w;
857 
859  t_rhs(L) = alpha * (t_res_flow(i, j) * t_L(i, j, L));
860  next();
861 
862  auto t_nf = getFTensor1FromArray<size_symm, size_symm>(nf);
863  size_t bb = 0;
864  for (; bb != AssemblyDomainEleOp::nbRows / size_symm; ++bb) {
865  t_nf(L) += t_base * t_rhs(L);
866  ++t_base;
867  ++t_nf;
868  }
869  for (; bb < nb_base_functions; ++bb)
870  ++t_base;
871  }
872 
874 }
875 
876 template <typename AssemblyDomainEleOp>
878  : public AssemblyDomainEleOp {
879  OpCalculateConstraintsRhsImpl(const std::string field_name,
880  boost::shared_ptr<CommonData> common_data_ptr,
881  boost::shared_ptr<MatrixDouble> m_D_ptr);
883 
884 private:
885  boost::shared_ptr<CommonData> commonDataPtr;
886  boost::shared_ptr<MatrixDouble> mDPtr;
887 };
888 
889 template <typename AssemblyDomainEleOp>
892  boost::shared_ptr<CommonData> common_data_ptr,
893  boost::shared_ptr<MatrixDouble> m_D_ptr)
895  commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {}
896 
897 template <typename AssemblyDomainEleOp>
902 
903  const size_t nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
904  const size_t nb_base_functions = data.getN().size2();
905 
906  auto t_res_c = getFTensor0FromVec(commonDataPtr->resC);
907 
908  auto next = [&]() { ++t_res_c; };
909 
910  auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
911  auto &nf = AssemblyDomainEleOp::locF;
912  auto t_base = data.getFTensor0N();
913  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
914  const double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
915  ++t_w;
916  const auto res = alpha * t_res_c;
917  next();
918 
919  size_t bb = 0;
920  for (; bb != AssemblyDomainEleOp::nbRows; ++bb) {
921  nf[bb] += t_base * res;
922  ++t_base;
923  }
924  for (; bb < nb_base_functions; ++bb)
925  ++t_base;
926  }
927 
929 }
930 
931 template <int DIM, typename AssemblyDomainEleOp>
933  : public AssemblyDomainEleOp {
935  const std::string row_field_name, const std::string col_field_name,
936  boost::shared_ptr<CommonData> common_data_ptr,
937  boost::shared_ptr<MatrixDouble> m_D_ptr);
938  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
939  EntitiesFieldData::EntData &col_data);
940 
941 private:
942  boost::shared_ptr<CommonData> commonDataPtr;
943  boost::shared_ptr<MatrixDouble> mDPtr;
944 };
945 
946 template <int DIM, typename AssemblyDomainEleOp>
949  const std::string row_field_name, const std::string col_field_name,
950  boost::shared_ptr<CommonData> common_data_ptr,
951  boost::shared_ptr<MatrixDouble> m_D_ptr)
952  : AssemblyDomainEleOp(row_field_name, col_field_name,
953  AssemblyDomainEleOp::OPROWCOL),
954  commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
955  AssemblyDomainEleOp::sYmm = false;
956 }
957 
958 static inline auto get_mat_tensor_sym_dtensor_sym(size_t rr, MatrixDouble &mat,
961  &mat(3 * rr + 0, 0), &mat(3 * rr + 0, 1), &mat(3 * rr + 0, 2),
962  &mat(3 * rr + 1, 0), &mat(3 * rr + 1, 1), &mat(3 * rr + 1, 2),
963  &mat(3 * rr + 2, 0), &mat(3 * rr + 2, 1), &mat(3 * rr + 2, 2)};
964 }
965 
966 static inline auto get_mat_tensor_sym_dtensor_sym(size_t rr, MatrixDouble &mat,
969  &mat(6 * rr + 0, 0), &mat(6 * rr + 0, 1), &mat(6 * rr + 0, 2),
970  &mat(6 * rr + 0, 3), &mat(6 * rr + 0, 4), &mat(6 * rr + 0, 5),
971  &mat(6 * rr + 1, 0), &mat(6 * rr + 1, 1), &mat(6 * rr + 1, 2),
972  &mat(6 * rr + 1, 3), &mat(6 * rr + 1, 4), &mat(6 * rr + 1, 5),
973  &mat(6 * rr + 2, 0), &mat(6 * rr + 2, 1), &mat(6 * rr + 2, 2),
974  &mat(6 * rr + 2, 3), &mat(6 * rr + 2, 4), &mat(6 * rr + 2, 5),
975  &mat(6 * rr + 3, 0), &mat(6 * rr + 3, 1), &mat(6 * rr + 3, 2),
976  &mat(6 * rr + 3, 3), &mat(6 * rr + 3, 4), &mat(6 * rr + 3, 5),
977  &mat(6 * rr + 4, 0), &mat(6 * rr + 4, 1), &mat(6 * rr + 4, 2),
978  &mat(6 * rr + 4, 3), &mat(6 * rr + 4, 4), &mat(6 * rr + 4, 5),
979  &mat(6 * rr + 5, 0), &mat(6 * rr + 5, 1), &mat(6 * rr + 5, 2),
980  &mat(6 * rr + 5, 3), &mat(6 * rr + 5, 4), &mat(6 * rr + 5, 5)};
981 }
982 
983 template <int DIM, typename AssemblyDomainEleOp>
986  EntitiesFieldData::EntData &row_data,
987  EntitiesFieldData::EntData &col_data) {
989 
994  constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
997 
998  auto &locMat = AssemblyDomainEleOp::locMat;
999 
1000  const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
1001  const auto nb_row_base_functions = row_data.getN().size2();
1002 
1003  auto t_res_flow_dstrain =
1004  getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->resFlowDstrain);
1005  auto t_res_flow_dplastic_strain =
1006  getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->resFlowDstrainDot);
1007  auto t_L = symm_L_tensor(FTensor::Number<DIM>());
1008 
1009  auto next = [&]() {
1010  ++t_res_flow_dstrain;
1011  ++t_res_flow_dplastic_strain;
1012  };
1013 
1014  auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
1015  auto t_row_base = row_data.getFTensor0N();
1016  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1017  double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
1018  ++t_w;
1019 
1021  t_res_mat(O, L) =
1022  alpha * (t_L(i, j, O) * ((t_res_flow_dplastic_strain(i, j, k, l) -
1023  t_res_flow_dstrain(i, j, k, l)) *
1024  t_L(k, l, L)));
1025  next();
1026 
1027  size_t rr = 0;
1028  for (; rr != AssemblyDomainEleOp::nbRows / size_symm; ++rr) {
1029  auto t_mat = get_mat_tensor_sym_dtensor_sym(rr, locMat,
1031  auto t_col_base = col_data.getFTensor0N(gg, 0);
1032  for (size_t cc = 0; cc != AssemblyDomainEleOp::nbCols / size_symm; ++cc) {
1033  t_mat(O, L) += ((t_row_base * t_col_base) * t_res_mat(O, L));
1034  ++t_mat;
1035  ++t_col_base;
1036  }
1037 
1038  ++t_row_base;
1039  }
1040 
1041  for (; rr < nb_row_base_functions; ++rr)
1042  ++t_row_base;
1043  }
1044 
1046 }
1047 
1048 template <int DIM, typename AssemblyDomainEleOp>
1050  : public AssemblyDomainEleOp {
1052  const std::string row_field_name, const std::string col_field_name,
1053  boost::shared_ptr<CommonData> common_data_ptr,
1054  boost::shared_ptr<MatrixDouble> m_D_ptr);
1055  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
1056  EntitiesFieldData::EntData &col_data);
1057 
1058 private:
1059  boost::shared_ptr<CommonData> commonDataPtr;
1060  boost::shared_ptr<MatrixDouble> mDPtr;
1061 };
1062 
1063 template <int DIM, typename AssemblyDomainEleOp>
1065  : public AssemblyDomainEleOp {
1067  const std::string row_field_name, const std::string col_field_name,
1068  boost::shared_ptr<CommonData> common_data_ptr,
1069  boost::shared_ptr<MatrixDouble> m_D_ptr);
1070  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
1071  EntitiesFieldData::EntData &col_data);
1072 
1073 private:
1074  boost::shared_ptr<CommonData> commonDataPtr;
1075  boost::shared_ptr<MatrixDouble> mDPtr;
1076 };
1077 
1078 template <int DIM, typename AssemblyDomainEleOp>
1081  const std::string row_field_name, const std::string col_field_name,
1082  boost::shared_ptr<CommonData> common_data_ptr,
1083  boost::shared_ptr<MatrixDouble> m_D_ptr)
1084  : AssemblyDomainEleOp(row_field_name, col_field_name,
1085  DomainEleOp::OPROWCOL),
1086  commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
1087  AssemblyDomainEleOp::sYmm = false;
1088 }
1089 
1090 static inline auto get_mat_tensor_sym_dscalar(size_t rr, MatrixDouble &mat,
1093  &mat(3 * rr + 0, 0), &mat(3 * rr + 1, 0), &mat(3 * rr + 2, 0)};
1094 }
1095 
1096 static inline auto get_mat_tensor_sym_dscalar(size_t rr, MatrixDouble &mat,
1099  &mat(6 * rr + 0, 0), &mat(6 * rr + 1, 0), &mat(6 * rr + 2, 0),
1100  &mat(6 * rr + 3, 0), &mat(6 * rr + 4, 0), &mat(6 * rr + 5, 0)};
1101 }
1102 
1103 template <int DIM, typename AssemblyDomainEleOp>
1106  EntitiesFieldData::EntData &row_data,
1107  EntitiesFieldData::EntData &col_data) {
1109 
1112  constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
1114 
1115  const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
1116  const size_t nb_row_base_functions = row_data.getN().size2();
1117  auto &locMat = AssemblyDomainEleOp::locMat;
1118 
1119  const auto type = AssemblyDomainEleOp::getFEType();
1120  const auto nb_nodes = moab::CN::VerticesPerEntity(type);
1121 
1122  auto t_res_flow_dtau =
1123  getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resFlowDtau);
1124 
1125  auto t_L = symm_L_tensor(FTensor::Number<DIM>());
1126 
1127  auto next = [&]() { ++t_res_flow_dtau; };
1128 
1129  auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
1130  auto t_row_base = row_data.getFTensor0N();
1131  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1132  double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
1133  ++t_w;
1135  t_res_vec(L) = alpha * (t_res_flow_dtau(i, j) * t_L(i, j, L));
1136  next();
1137 
1138  size_t rr = 0;
1139  for (; rr != AssemblyDomainEleOp::nbRows / size_symm; ++rr) {
1140  auto t_mat =
1142  auto t_col_base = col_data.getFTensor0N(gg, 0);
1143  for (size_t cc = 0; cc != AssemblyDomainEleOp::nbCols; cc++) {
1144  t_mat(L) += t_row_base * t_col_base * t_res_vec(L);
1145  ++t_mat;
1146  ++t_col_base;
1147  }
1148  ++t_row_base;
1149  }
1150  for (; rr != nb_row_base_functions; ++rr)
1151  ++t_row_base;
1152  }
1153 
1155 }
1156 
1157 template <int DIM, typename AssemblyDomainEleOp>
1159  : public AssemblyDomainEleOp {
1161  const std::string row_field_name, const std::string col_field_name,
1162  boost::shared_ptr<CommonData> common_data_ptr,
1163  boost::shared_ptr<MatrixDouble> mat_D_ptr);
1164  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
1165  EntitiesFieldData::EntData &col_data);
1166 
1167 private:
1168  boost::shared_ptr<CommonData> commonDataPtr;
1169  boost::shared_ptr<MatrixDouble> mDPtr;
1170 };
1171 
1172 template <int DIM, typename AssemblyDomainEleOp>
1175  const std::string row_field_name, const std::string col_field_name,
1176  boost::shared_ptr<CommonData> common_data_ptr,
1177  boost::shared_ptr<MatrixDouble> m_D_ptr)
1178  : AssemblyDomainEleOp(row_field_name, col_field_name,
1179  DomainEleOp::OPROWCOL),
1180  commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
1181  AssemblyDomainEleOp::sYmm = false;
1182 }
1183 
1186  &mat(0, 0), &mat(0, 1), &mat(0, 2)};
1187 }
1188 
1191  &mat(0, 0), &mat(0, 1), &mat(0, 2), &mat(0, 3), &mat(0, 4), &mat(0, 5)};
1192 }
1193 
1194 template <int DIM, typename AssemblyDomainEleOp>
1197  EntitiesFieldData::EntData &row_data,
1198  EntitiesFieldData::EntData &col_data) {
1200 
1203  constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
1205 
1206  const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
1207  const auto nb_row_base_functions = row_data.getN().size2();
1208 
1209  auto t_c_dstrain =
1210  getFTensor2SymmetricFromMat<SPACE_DIM>(commonDataPtr->resCdStrain);
1211  auto t_c_dplastic_strain =
1212  getFTensor2SymmetricFromMat<SPACE_DIM>(commonDataPtr->resCdPlasticStrain);
1213 
1214  auto next = [&]() {
1215  ++t_c_dstrain;
1216  ++t_c_dplastic_strain;
1217  };
1218 
1220 
1221  auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
1222  auto t_row_base = row_data.getFTensor0N();
1223  for (auto gg = 0; gg != nb_integration_pts; ++gg) {
1224  const double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
1225  ++t_w;
1226 
1228  t_res_vec(L) =
1229  t_L(i, j, L) * (t_c_dplastic_strain(i, j) - t_c_dstrain(i, j));
1230  next();
1231 
1234  size_t rr = 0;
1235  for (; rr != AssemblyDomainEleOp::nbRows; ++rr) {
1236  const auto row_base = alpha * t_row_base;
1237  auto t_col_base = col_data.getFTensor0N(gg, 0);
1238  for (size_t cc = 0; cc != AssemblyDomainEleOp::nbCols / size_symm; cc++) {
1239  t_mat(L) += (row_base * t_col_base) * t_res_vec(L);
1240  ++t_mat;
1241  ++t_col_base;
1242  }
1243  ++t_row_base;
1244  }
1245  for (; rr != nb_row_base_functions; ++rr)
1246  ++t_row_base;
1247  }
1248 
1250 }
1251 
1252 template <typename AssemblyDomainEleOp>
1254  : public AssemblyDomainEleOp {
1256  const std::string row_field_name, const std::string col_field_name,
1257  boost::shared_ptr<CommonData> common_data_ptr);
1258  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
1259  EntitiesFieldData::EntData &col_data);
1260 
1261 private:
1262  boost::shared_ptr<CommonData> commonDataPtr;
1263 };
1264 
1265 template <typename AssemblyDomainEleOp>
1268  const std::string row_field_name, const std::string col_field_name,
1269  boost::shared_ptr<CommonData> common_data_ptr)
1270  : AssemblyDomainEleOp(row_field_name, col_field_name,
1271  DomainEleOp::OPROWCOL),
1272  commonDataPtr(common_data_ptr) {
1273  AssemblyDomainEleOp::sYmm = false;
1274 }
1275 
1276 template <typename AssemblyDomainEleOp>
1279  EntitiesFieldData::EntData &row_data,
1280  EntitiesFieldData::EntData &col_data) {
1282 
1283  const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
1284  const auto nb_row_base_functions = row_data.getN().size2();
1285 
1286  auto t_res_c_dtau = getFTensor0FromVec(commonDataPtr->resCdTau);
1287  auto next = [&]() { ++t_res_c_dtau; };
1288 
1289  auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
1290  auto t_row_base = row_data.getFTensor0N();
1291  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1292  const double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
1293  ++t_w;
1294 
1295  const auto res = alpha * (t_res_c_dtau);
1296  next();
1297 
1298  auto mat_ptr = AssemblyDomainEleOp::locMat.data().begin();
1299  size_t rr = 0;
1300  for (; rr != AssemblyDomainEleOp::nbRows; ++rr) {
1301  auto t_col_base = col_data.getFTensor0N(gg, 0);
1302  for (size_t cc = 0; cc != AssemblyDomainEleOp::nbCols; ++cc) {
1303  *mat_ptr += t_row_base * t_col_base * res;
1304  ++t_col_base;
1305  ++mat_ptr;
1306  }
1307  ++t_row_base;
1308  }
1309  for (; rr < nb_row_base_functions; ++rr)
1310  ++t_row_base;
1311  }
1312 
1314 }
1315 }; // namespace PlasticOps
PlasticOps::diff_plastic_flow_dstress
auto diff_plastic_flow_dstress(long double f, FTensor::Tensor2_symmetric< T, DIM > &t_flow, FTensor::Ddg< double, 3, DIM > &&t_diff_deviator)
Definition: PlasticOpsGeneric.hpp:208
PlasticOps::OpCalculateConstraintsRhsImpl
Definition: PlasticOps.hpp:133
NOSPACE
@ NOSPACE
Definition: definitions.h:83
PlasticOps::CommonData::BISO
@ BISO
Definition: PlasticOps.hpp:52
PlasticOps::OpCalculateConstraintsLhs_dTAUImpl
Definition: PlasticOps.hpp:163
PlasticOps::OpCalculatePlasticityImpl< DIM, GAUSS, DomainEleOp >::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOpsGeneric.hpp:464
ContactOps::sign
double sign(double x)
Definition: ContactOps.hpp:551
PlasticOps::get_mat_tensor_sym_dtensor_sym
static auto get_mat_tensor_sym_dtensor_sym(size_t rr, MatrixDouble &mat, FTensor::Number< 2 >)
Definition: PlasticOpsGeneric.hpp:958
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
PlasticOps::OpCalculateConstraintsRhsImpl< GAUSS, AssemblyDomainEleOp >::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOpsGeneric.hpp:885
PlasticOps::OpCalculateConstraintsLhs_dUImpl< DIM, GAUSS, AssemblyDomainEleOp >::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOpsGeneric.hpp:1059
EshelbianPlasticity::size_symm
constexpr static auto size_symm
Definition: EshelbianOperators.cpp:47
DEPRECATED
#define DEPRECATED
Definition: definitions.h:17
PlasticOps::w
double w(double eqiv, double dot_tau, double f, double sigma_y, double sigma_Y)
Definition: PlasticOpsGeneric.hpp:276
PlasticOps::J
FTensor::Index< 'J', 3 > J
Definition: PlasticOps.hpp:116
PlasticOps::OpCalculatePlasticFlowRhsImpl< DIM, GAUSS, AssemblyDomainEleOp >::mDPtr
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOpsGeneric.hpp:818
PlasticOps::OpPlasticStressImpl< DIM, GAUSS, DomainEleOp >::mDPtr
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOpsGeneric.hpp:752
PlasticOps::OpCalculateConstraintsLhs_dEPImpl< DIM, GAUSS, AssemblyDomainEleOp >::mDPtr
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOpsGeneric.hpp:1169
MoFEM::OpBaseImpl::locMat
MatrixDouble locMat
local entity block matrix
Definition: FormsIntegrators.hpp:247
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
PlasticOps::diff_constrain_ddot_tau
double diff_constrain_ddot_tau(double sign, double eqiv, double dot_tau, double vis_H, double sigma_Y)
Definition: PlasticOpsGeneric.hpp:300
PlasticOps::OpCalculateConstraintsLhs_dEPImpl< DIM, GAUSS, AssemblyDomainEleOp >::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOpsGeneric.hpp:1168
PlasticOps::diff_deviator
auto diff_deviator(FTensor::Ddg< double, DIM, DIM > &&t_diff_stress, FTensor::Number< DIM >)
Definition: PlasticOpsGeneric.hpp:125
PlasticOps::diff_constrain_df
auto diff_constrain_df(double sign)
Definition: PlasticOpsGeneric.hpp:310
PlasticOps::diff_tensor
auto diff_tensor(FTensor::Number< DIM >)
[Lambda functions]
Definition: PlasticOpsGeneric.hpp:10
PlasticOps::CommonData::C1_k
@ C1_k
Definition: PlasticOps.hpp:53
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
PlasticOps::equivalent_strain_dot
auto equivalent_strain_dot(FTensor::Tensor2_symmetric< T, DIM > &t_plastic_strain_dot)
Definition: PlasticOpsGeneric.hpp:338
PlasticOps::diff_equivalent_strain_dot
auto diff_equivalent_strain_dot(const T1 eqiv, T2 &t_plastic_strain_dot, T3 &t_diff_plastic_strain, FTensor::Number< DIM >)
Definition: PlasticOpsGeneric.hpp:349
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
iso_hardening_dtau
double iso_hardening_dtau(double tau, double H, double Qinf, double b_iso)
Definition: plastic.cpp:77
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:13
PlasticOps::diff_constrain_dstress
auto diff_constrain_dstress(double diff_constrain_df, FTensor::Tensor2_symmetric< T, DIM > &t_plastic_flow)
Definition: PlasticOpsGeneric.hpp:316
zeta
double zeta
Viscous hardening.
Definition: plastic.cpp:126
FTensor::Tensor2
Definition: Tensor2_value.hpp:16
PlasticOps::N
FTensor::Index< 'N', 3 > N
Definition: PlasticOps.hpp:118
MoFEM::OpBaseImpl::locF
VectorDouble locF
local entity vector
Definition: FormsIntegrators.hpp:249
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
PlasticOps::OpPlasticStressImpl< DIM, GAUSS, DomainEleOp >::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOpsGeneric.hpp:753
PlasticOps::M
FTensor::Index< 'M', 3 > M
Definition: PlasticOps.hpp:117
PlasticOps::CommonData::H
@ H
Definition: PlasticOps.hpp:49
PlasticOps::diff_constrain_deqiv
double diff_constrain_deqiv(double sign, double eqiv, double dot_tau, double sigma_Y)
Definition: PlasticOpsGeneric.hpp:305
PlasticOps::platsic_surface
double platsic_surface(FTensor::Tensor2_symmetric< double, 3 > &&t_stress_deviator)
Definition: PlasticOpsGeneric.hpp:189
MoFEM::OpBaseImpl
Definition: FormsIntegrators.hpp:178
FTensor::Number
Definition: Number.hpp:11
PlasticOps::get_mat_scalar_dtensor_sym
auto get_mat_scalar_dtensor_sym(MatrixDouble &mat, FTensor::Number< 2 >)
Definition: PlasticOpsGeneric.hpp:1184
PlasticOps::OpCalculatePlasticityImpl
Definition: PlasticOps.hpp:124
a
constexpr double a
Definition: approx_sphere.cpp:30
PlasticOps::CommonData::activityData
static std::array< int, 5 > activityData
Definition: PlasticOps.hpp:107
kinematic_hardening
auto kinematic_hardening(FTensor::Tensor2_symmetric< T, DIM > &t_plastic_strain, double C1_k)
Definition: plastic.cpp:91
convert.type
type
Definition: convert.py:64
PlasticOps::diff_constrain_dstrain
auto diff_constrain_dstrain(T1 &t_D, T2 &&t_diff_constrain_dstress)
Definition: PlasticOpsGeneric.hpp:326
PlasticOps::constrain_abs
double constrain_abs(double x)
Definition: PlasticOpsGeneric.hpp:265
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
cn1
double cn1
Definition: plastic.cpp:132
PlasticOps::deviator
auto deviator(FTensor::Tensor2_symmetric< T, DIM > &t_stress, double trace, FTensor::Tensor2_symmetric< double, DIM > &t_alpha, FTensor::Number< DIM >)
Definition: PlasticOpsGeneric.hpp:92
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
PlasticOps::OpCalculatePlasticFlowLhs_dTAUImpl< DIM, GAUSS, AssemblyDomainEleOp >::mDPtr
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOpsGeneric.hpp:1075
PlasticOps::OpCalculatePlasticFlowLhs_dTAUImpl< DIM, GAUSS, AssemblyDomainEleOp >::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOpsGeneric.hpp:1074
cn0
double cn0
Definition: plastic.cpp:131
MoFEM::L
VectorDouble L
Definition: Projection10NodeCoordsOnField.cpp:124
PlasticOps::OpCalculateConstraintsLhs_dTAUImpl< GAUSS, AssemblyDomainEleOp >::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOpsGeneric.hpp:1262
PlasticOps::get_mat_tensor_sym_dvector
static auto get_mat_tensor_sym_dvector(size_t rr, MatrixDouble &mat, FTensor::Number< 2 >)
[Lambda functions]
Definition: PlasticOpsGeneric.hpp:366
FTensor::Tensor4
Definition: Tensor4_value.hpp:18
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FTensor::Index< 'i', DIM >
PlasticOps::CommonData::QINF
@ QINF
Definition: PlasticOps.hpp:51
PlasticOps::OpCalculatePlasticFlowRhsImpl< DIM, GAUSS, AssemblyDomainEleOp >::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOpsGeneric.hpp:817
PlasticOps::plastic_flow
auto plastic_flow(long double f, FTensor::Tensor2_symmetric< double, 3 > &&t_dev_stress, FTensor::Ddg< double, 3, DIM > &&t_diff_deviator)
Definition: PlasticOpsGeneric.hpp:195
convert.n
n
Definition: convert.py:82
PlasticOps::OpCalculatePlasticFlowLhs_dTAUImpl
Definition: PlasticOps.hpp:151
MoFEM::OpBaseImpl::nbCols
int nbCols
number if dof on column
Definition: FormsIntegrators.hpp:235
DomainEleOp
PlasticOps
Definition: PlasticNaturalBCs.hpp:13
MoFEM::OpBaseImpl::nbRows
int nbRows
number of dofs on rows
Definition: FormsIntegrators.hpp:234
PlasticOps::OpCalculateConstraintsRhsImpl< GAUSS, AssemblyDomainEleOp >::mDPtr
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOpsGeneric.hpp:886
FTensor::Dg
Definition: Dg_value.hpp:9
PlasticOps::OpCalculatePlasticityImpl< DIM, GAUSS, DomainEleOp >::mDPtr
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOpsGeneric.hpp:465
PlasticOps::symm_L_tensor
auto symm_L_tensor(FTensor::Number< DIM >)
Definition: PlasticOpsGeneric.hpp:21
iso_hardening
double iso_hardening(double tau, double H, double Qinf, double b_iso, double sigmaY)
Definition: plastic.cpp:72
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
PlasticOps::OpCalculatePlasticFlowLhs_dEPImpl
Definition: PlasticOps.hpp:148
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
FTensor::Ddg
Definition: Ddg_value.hpp:7
PlasticOps::get_mat_tensor_sym_dscalar
static auto get_mat_tensor_sym_dscalar(size_t rr, MatrixDouble &mat, FTensor::Number< 2 >)
Definition: PlasticOpsGeneric.hpp:1090
PlasticOps::OpCalculateConstraintsLhs_dEPImpl
Definition: PlasticOps.hpp:160
PlasticOps::OpCalculatePlasticSurfaceImpl
Definition: PlasticOps.hpp:121
PlasticOps::OpCalculateConstraintsLhs_dUImpl< DIM, GAUSS, AssemblyDomainEleOp >::mDPtr
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOpsGeneric.hpp:1060
PlasticOps::constraint
double constraint(double eqiv, double dot_tau, double f, double sigma_y, double abs_w, double vis_H, double sigma_Y)
Definition: PlasticOpsGeneric.hpp:292
PlasticOps::diff_symmetrize
auto diff_symmetrize(FTensor::Number< DIM >)
Definition: PlasticOpsGeneric.hpp:43
PlasticOps::CommonData::VIS_H
@ VIS_H
Definition: PlasticOps.hpp:50
PlasticOps::OpCalculatePlasticFlowLhs_dEPImpl< DIM, GAUSS, AssemblyDomainEleOp >::mDPtr
boost::shared_ptr< MatrixDouble > mDPtr
Definition: PlasticOpsGeneric.hpp:943
PlasticOps::constrian_sign
double constrian_sign(double x)
Definition: PlasticOpsGeneric.hpp:251
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
PlasticOps::trace
double trace(FTensor::Tensor2_symmetric< T, 2 > &t_stress)
Definition: PlasticOpsGeneric.hpp:80
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
third
constexpr double third
Definition: EshelbianADOL-C.cpp:14
PlasticOps::diff_plastic_flow_dstrain
auto diff_plastic_flow_dstrain(FTensor::Ddg< T, DIM, DIM > &t_D, FTensor::Ddg< double, DIM, DIM > &&t_diff_plastic_flow_dstress)
Definition: PlasticOpsGeneric.hpp:224
PlasticOps::I
FTensor::Index< 'I', 3 > I
[Common data]
Definition: PlasticOps.hpp:115
PlasticOps::OpPlasticStressImpl
Definition: PlasticOps.hpp:127
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
PlasticOps::OpCalculatePlasticFlowLhs_dEPImpl< DIM, GAUSS, AssemblyDomainEleOp >::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOpsGeneric.hpp:942
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
PlasticOps::diff_constrain_dsigma_y
auto diff_constrain_dsigma_y(double sign)
Definition: PlasticOpsGeneric.hpp:312
PlasticOps::OpCalculatePlasticFlowRhsImpl
Definition: PlasticOps.hpp:130
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
PlasticOps::OpCalculatePlasticSurfaceImpl< DIM, GAUSS, DomainEleOp >::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: PlasticOpsGeneric.hpp:394
PlasticOps::CommonData::SIGMA_Y
@ SIGMA_Y
Definition: PlasticOps.hpp:48
PlasticOps::OpCalculateConstraintsLhs_dUImpl
Definition: PlasticOps.hpp:154