GetFEM  5.4.3
getfem_plasticity.cc
1 /*===========================================================================
2 
3  Copyright (C) 2002-2020 Amandine Cottaz, Yves Renard
4  Copyright (C) 2014-2020 Konstantinos Poulios
5 
6  This file is a part of GetFEM
7 
8  GetFEM is free software; you can redistribute it and/or modify it
9  under the terms of the GNU Lesser General Public License as published
10  by the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version along with the GCC Runtime Library
12  Exception either version 3.1 or (at your option) any later version.
13  This program is distributed in the hope that it will be useful, but
14  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  License and GCC Runtime Library Exception for more details.
17  You should have received a copy of the GNU Lesser General Public License
18  along with this program; if not, write to the Free Software Foundation,
19  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
20 
21 ===========================================================================*/
22 
23 
24 #include "getfem/getfem_models.h"
25 #include "getfem/getfem_plasticity.h"
26 #include "getfem/getfem_interpolation.h"
27 #include "getfem/getfem_generic_assembly.h"
28 #include "gmm/gmm_dense_matrix_functions.h"
29 #include <iomanip>
30 
31 namespace getfem {
32 
33  //=========================================================================
34  //
35  // Specific nonlinear operators of the high-level generic assembly
36  // language, useful for plasticity modeling
37  //
38  //=========================================================================
39 
40  static void ga_init_scalar(bgeot::multi_index &mi) { mi.resize(0); }
41  static void ga_init_vector(bgeot::multi_index &mi, size_type N)
42  { mi.resize(1); mi[0] = N; }
43  static void ga_init_matrix(bgeot::multi_index &mi, size_type M, size_type N)
44  { mi.resize(2); mi[0] = M; mi[1] = N; }
45  static void ga_init_square_matrix(bgeot::multi_index &mi, size_type N)
46  { mi.resize(2); mi[0] = mi[1] = N; }
47 
48 
49  bool expm(const base_matrix &a_, base_matrix &aexp, scalar_type tol=1e-15) {
50 
51  const size_type itmax = 40;
52  base_matrix a(a_);
53  // scale input matrix a
54  int e;
55  frexp(gmm::mat_norminf(a), &e);
56  e = std::max(0, std::min(1023, e));
57  gmm::scale(a, pow(scalar_type(2),-scalar_type(e)));
58 
59  base_matrix atmp(a), an(a);
60  gmm::copy(a, aexp);
61  gmm::add(gmm::identity_matrix(), aexp);
62  scalar_type factn(1);
63  bool success(false);
64  for (size_type n=2; n < itmax; ++n) {
65  factn /= scalar_type(n);
66  gmm::mult(an, a, atmp);
67  gmm::copy(atmp, an);
68  gmm::scale(atmp, factn);
69  gmm::add(atmp, aexp);
70  if (gmm::mat_euclidean_norm(atmp) < tol) {
71  success = true;
72  break;
73  }
74  }
75  // unscale result
76  for (int i=0; i < e; ++i) {
77  gmm::mult(aexp, aexp, atmp);
78  gmm::copy(atmp, aexp);
79  }
80  return success;
81  }
82 
83  bool expm_deriv(const base_matrix &a_, base_tensor &daexp,
84  base_matrix *paexp=NULL, scalar_type tol=1e-15) {
85 
86  const size_type itmax = 40;
87  size_type N = gmm::mat_nrows(a_);
88  size_type N2 = N*N;
89 
90  base_matrix a(a_);
91  // scale input matrix a
92  int e;
93  frexp(gmm::mat_norminf(a), &e);
94  e = std::max(0, std::min(1023, e));
95  scalar_type scale = pow(scalar_type(2),-scalar_type(e));
96  gmm::scale(a, scale);
97 
98  base_vector factnn(itmax);
99  base_matrix atmp(a), an(a), aexp(a);
100  base_tensor ann(bgeot::multi_index(N,N,itmax));
101  gmm::add(gmm::identity_matrix(), aexp);
102  gmm::copy(gmm::identity_matrix(), atmp);
103  std::copy(atmp.begin(), atmp.end(), ann.begin());
104  factnn[1] = 1;
105  std::copy(a.begin(), a.end(), ann.begin()+N2);
106  size_type n;
107  bool success(false);
108  for (n=2; n < itmax; ++n) {
109  factnn[n] = factnn[n-1]/scalar_type(n);
110  gmm::mult(an, a, atmp);
111  gmm::copy(atmp, an);
112  std::copy(an.begin(), an.end(), ann.begin()+n*N2);
113  gmm::scale(atmp, factnn[n]);
114  gmm::add(atmp, aexp);
115  if (gmm::mat_euclidean_norm(atmp) < tol) {
116  success = true;
117  break;
118  }
119  }
120 
121  if (!success)
122  return false;
123 
124  gmm::clear(daexp.as_vector());
125  gmm::scale(factnn, scale);
126  for (--n; n >= 1; --n) {
127  scalar_type factn = factnn[n];
128  for (size_type m=1; m <= n; ++m)
129  for (size_type l=0; l < N; ++l)
130  for (size_type k=0; k < N; ++k)
131  for (size_type j=0; j < N; ++j)
132  for (size_type i=0; i < N; ++i)
133  daexp(i,j,k,l) += factn*ann(i,k,m-1)*ann(l,j,n-m);
134  }
135 
136  // unscale result
137  base_matrix atmp1(a), atmp2(a);
138  for (int i=0; i < e; ++i) {
139  for (size_type l=0; l < N; ++l)
140  for (size_type k=0; k < N; ++k) {
141  std::copy(&daexp(0,0,k,l), &daexp(0,0,k,l)+N*N, atmp.begin());
142  gmm::mult(atmp, aexp, atmp1);
143  gmm::mult(aexp, atmp, atmp2);
144  gmm::add(atmp1, atmp2, atmp);
145  std::copy(atmp.begin(), atmp.end(), &daexp(0,0,k,l));
146  }
147  gmm::mult(aexp, aexp, atmp);
148  gmm::copy(atmp, aexp);
149  }
150 
151  if (paexp) gmm::copy(aexp, *paexp);
152  return true;
153  }
154 
155  // numerical differantiation of logm
156  // not used because it caused some issues and was slower than
157  // simply inverting the derivative of expm
158  bool logm_deriv(const base_matrix &a, base_tensor &dalog,
159  base_matrix *palog=NULL) {
160 
161  size_type N = gmm::mat_nrows(a);
162  base_matrix a1(a), alog1(a), alog(a);
163  logm(a, alog);
164  scalar_type eps(1e-8);
165  for (size_type k=0; k < N; ++k)
166  for (size_type l=0; l < N; ++l) {
167  gmm::copy(a, a1);
168  a1(k,l) += eps;
169  gmm::logm(a1, alog1);
170  for (size_type i=0; i < N; ++i)
171  for (size_type j=0; j < N; ++j)
172  dalog(i,j,k,l) = (alog1(i,j) - alog(i,j))/eps;
173  }
174  if (palog) gmm::copy(alog, *palog);
175  return true;
176  }
177 
178 
179  // Matrix exponential
180  struct matrix_exponential_operator : public ga_nonlinear_operator {
181  bool result_size(const arg_list &args, bgeot::multi_index &sizes) const {
182  if (args.size() != 1 || args[0]->sizes().size() != 2
183  || args[0]->sizes()[0] != args[0]->sizes()[1]) return false;
184  ga_init_matrix(sizes, args[0]->sizes()[0], args[0]->sizes()[1]);
185  return true;
186  }
187 
188  // Value:
189  void value(const arg_list &args, base_tensor &result) const {
190  size_type N = args[0]->sizes()[0];
191  base_matrix inpmat(N,N), outmat(N,N);
192  gmm::copy(args[0]->as_vector(), inpmat.as_vector());
193  bool info = expm(inpmat, outmat);
194  GMM_ASSERT1(info, "Matrix exponential calculation "
195  "failed to converge");
196  gmm::copy(outmat.as_vector(), result.as_vector());
197  }
198 
199  // Derivative:
200  void derivative(const arg_list &args, size_type /*nder*/,
201  base_tensor &result) const {
202  size_type N = args[0]->sizes()[0];
203  base_matrix inpmat(N,N);
204  gmm::copy(args[0]->as_vector(), inpmat.as_vector());
205  bool info = expm_deriv(inpmat, result);
206  GMM_ASSERT1(info, "Matrix exponential derivative calculation "
207  "failed to converge");
208  }
209 
210  // Second derivative : not implemented
211  void second_derivative(const arg_list &, size_type, size_type,
212  base_tensor &) const {
213  GMM_ASSERT1(false, "Sorry, second derivative not implemented");
214  }
215  };
216 
217 
218  // Matrix logarithm
219  struct matrix_logarithm_operator : public ga_nonlinear_operator {
220  bool result_size(const arg_list &args, bgeot::multi_index &sizes) const {
221  if (args.size() != 1 || args[0]->sizes().size() != 2
222  || args[0]->sizes()[0] != args[0]->sizes()[1]) return false;
223  ga_init_matrix(sizes, args[0]->sizes()[0], args[0]->sizes()[1]);
224  return true;
225  }
226 
227  // Value:
228  void value(const arg_list &args, base_tensor &result) const {
229  size_type N = args[0]->sizes()[0];
230  base_matrix inpmat(N,N), outmat(N,N);
231  gmm::copy(args[0]->as_vector(), inpmat.as_vector());
232  gmm::logm(inpmat, outmat);
233  gmm::copy(outmat.as_vector(), result.as_vector());
234  }
235 
236  // Derivative:
237  void derivative(const arg_list &args, size_type /*nder*/,
238  base_tensor &result) const {
239  size_type N = args[0]->sizes()[0];
240  base_matrix inpmat(N,N), outmat(N,N), tmpmat(N*N,N*N);
241  gmm::copy(args[0]->as_vector(), inpmat.as_vector());
242  gmm::logm(inpmat, outmat);
243  bool info = expm_deriv(outmat, result);
244  if (info) {
245  gmm::copy(result.as_vector(), tmpmat.as_vector());
246  scalar_type det = gmm::lu_inverse(tmpmat);
247  if (det <= 0) gmm::copy(gmm::identity_matrix(), tmpmat);
248  gmm::copy(tmpmat.as_vector(), result.as_vector());
249  }
250  GMM_ASSERT1(info, "Matrix logarithm derivative calculation "
251  "failed to converge");
252  }
253 
254  // Second derivative : not implemented
255  void second_derivative(const arg_list &, size_type, size_type,
256  base_tensor &) const {
257  GMM_ASSERT1(false, "Sorry, second derivative not implemented");
258  }
259  };
260 
261 
262  // Normalized vector/matrix operator : Vector/matrix divided by its Frobenius norm
263  struct normalized_operator : public ga_nonlinear_operator {
264  bool result_size(const arg_list &args, bgeot::multi_index &sizes) const {
265  if (args.size() != 1 || args[0]->sizes().size() > 2
266  || args[0]->sizes().size() < 1) return false;
267  if (args[0]->sizes().size() == 1)
268  ga_init_vector(sizes, args[0]->sizes()[0]);
269  else
270  ga_init_matrix(sizes, args[0]->sizes()[0], args[0]->sizes()[1]);
271  return true;
272  }
273 
274  // Value : u/|u|
275  void value(const arg_list &args, base_tensor &result) const {
276 # if 1
277  const base_tensor &t = *args[0];
278  scalar_type eps = 1E-25;
279  scalar_type no = ::sqrt(gmm::vect_norm2_sqr(t.as_vector())+gmm::sqr(eps));
280  gmm::copy(gmm::scaled(t.as_vector(), scalar_type(1)/no),
281  result.as_vector());
282 # else
283  scalar_type no = gmm::vect_norm2(args[0]->as_vector());
284  if (no < 1E-15)
285  gmm::clear(result.as_vector());
286  else
287  gmm::copy(gmm::scaled(args[0]->as_vector(), scalar_type(1)/no),
288  result.as_vector());
289 # endif
290  }
291 
292  // Derivative : (|u|^2 Id - u x u)/|u|^3
293  void derivative(const arg_list &args, size_type,
294  base_tensor &result) const {
295  const base_tensor &t = *args[0];
296  size_type N = t.size();
297 # if 1
298  scalar_type eps = 1E-25;
299  scalar_type no = ::sqrt(gmm::vect_norm2_sqr(t.as_vector())+gmm::sqr(eps));
300  scalar_type no3 = no*no*no;
301 
302  gmm::clear(result.as_vector());
303  for (size_type i = 0; i < N; ++i) {
304  result[i*N+i] += scalar_type(1)/no;
305  for (size_type j = 0; j < N; ++j)
306  result[j*N+i] -= t[i]*t[j] / no3;
307  }
308 # else
309  scalar_type no = gmm::vect_norm2(t.as_vector());
310 
311  gmm::clear(result.as_vector());
312  if (no >= 1E-15) {
313  scalar_type no3 = no*no*no;
314  for (size_type i = 0; i < N; ++i) {
315  result[i*N+i] += scalar_type(1)/no;
316  for (size_type j = 0; j < N; ++j)
317  result[j*N+i] -= t[i]*t[j] / no3;
318  }
319  }
320 # endif
321  }
322 
323  // Second derivative : not implemented
324  void second_derivative(const arg_list &/*args*/, size_type, size_type,
325  base_tensor &/*result*/) const {
326  GMM_ASSERT1(false, "Sorry, second derivative not implemented");
327  }
328  };
329 
330 
331  // Ball Projection operator.
332  struct Ball_projection_operator : public ga_nonlinear_operator {
333  bool result_size(const arg_list &args, bgeot::multi_index &sizes) const {
334  if (args.size() != 2 || args[0]->sizes().size() > 2
335  || (args[0]->sizes().size() < 1 && args[0]->size() != 1)
336  || args[1]->size() != 1) return false;
337  if (args[0]->sizes().size() < 1)
338  ga_init_scalar(sizes);
339  else if (args[0]->sizes().size() == 1)
340  ga_init_vector(sizes, args[0]->sizes()[0]);
341  else
342  ga_init_matrix(sizes, args[0]->sizes()[0], args[0]->sizes()[1]);
343  return true;
344  }
345 
346  // Value : ru/|u| if |u| > r, else u
347  void value(const arg_list &args, base_tensor &result) const {
348  const base_tensor &t = *args[0];
349  scalar_type r = (*args[1])[0];
350  scalar_type no = gmm::vect_norm2(t.as_vector());
351  if (no > r)
352  gmm::copy(gmm::scaled(t.as_vector(), r/no), result.as_vector());
353  else
354  gmm::copy(t.as_vector(), result.as_vector());
355  }
356 
357  // Derivative
358  void derivative(const arg_list &args, size_type n,
359  base_tensor &result) const {
360  const base_tensor &t = *args[0];
361  size_type N = t.size();
362  scalar_type r = (*args[1])[0];
363  scalar_type no = gmm::vect_norm2(t.as_vector()), rno3 = r/(no*no*no);
364 
365  gmm::clear(result.as_vector());
366 
367  switch(n) {
368 
369  case 1 : // derivative with respect to u
370  if (r > 0) {
371  if (no <= r) {
372  for (size_type i = 0; i < N; ++i)
373  result[i*N+i] += scalar_type(1);
374  } else {
375  for (size_type i = 0; i < N; ++i) {
376  result[i*N+i] += r/no;
377  for (size_type j = 0; j < N; ++j)
378  result[j*N+i] -= t[i]*t[j]*rno3;
379  }
380  }
381  }
382  break;
383  case 2 : // derivative with respect to r
384  if (r > 0 && no > r) {
385  for (size_type i = 0; i < N; ++i)
386  result[i] = t[i]/no;
387  }
388  break;
389  default : GMM_ASSERT1(false, "Wrong derivative number");
390  }
391  }
392 
393  // Second derivative : not implemented
394  void second_derivative(const arg_list &/*args*/, size_type, size_type,
395  base_tensor &/*result*/) const {
396  GMM_ASSERT1(false, "Sorry, second derivative not implemented");
397  }
398  };
399 
400 
401  // Normalized_reg vector/matrix operator : Regularized Vector/matrix divided by its Frobenius norm
402  struct normalized_reg_operator : public ga_nonlinear_operator {
403  bool result_size(const arg_list &args, bgeot::multi_index &sizes) const {
404  if (args.size() != 2 || args[0]->sizes().size() > 2
405  || args[0]->sizes().size() < 1 || args[1]->size() != 1) return false;
406  if (args[0]->sizes().size() == 1)
407  ga_init_vector(sizes, args[0]->sizes()[0]);
408  else
409  ga_init_matrix(sizes, args[0]->sizes()[0], args[0]->sizes()[1]);
410  return true;
411  }
412 
413  // Value : u/(sqrt([u|^2+\eps^2))
414  void value(const arg_list &args, base_tensor &result) const {
415  const base_tensor &t = *args[0];
416  scalar_type eps = (*args[1])[0];
417  scalar_type no = ::sqrt(gmm::vect_norm2_sqr(t.as_vector())+gmm::sqr(eps));
418  gmm::copy(gmm::scaled(t.as_vector(), scalar_type(1)/no),
419  result.as_vector());
420  }
421 
422  // Derivative / u : ((|u|^2+eps^2) Id - u x u)/(|u|^2+eps^2)^(3/2)
423  // Derivative / eps : -eps*u/(|u|^2+eps^2)^(3/2)
424  void derivative(const arg_list &args, size_type nder,
425  base_tensor &result) const {
426  const base_tensor &t = *args[0];
427  scalar_type eps = (*args[1])[0];
428 
429  size_type N = t.size();
430  scalar_type no = ::sqrt(gmm::vect_norm2_sqr(t.as_vector())+gmm::sqr(eps));
431  scalar_type no3 = no*no*no;
432 
433  switch (nder) {
434  case 1:
435  gmm::clear(result.as_vector());
436  for (size_type i = 0; i < N; ++i) {
437  result[i*N+i] += scalar_type(1)/no;
438  for (size_type j = 0; j < N; ++j)
439  result[j*N+i] -= t[i]*t[j] / no3;
440  }
441  break;
442 
443  case 2:
444  gmm::copy(gmm::scaled(t.as_vector(), -scalar_type(eps)/no3),
445  result.as_vector());
446  break;
447  }
448  }
449 
450  // Second derivative : not implemented
451  void second_derivative(const arg_list &/*args*/, size_type, size_type,
452  base_tensor &/*result*/) const {
453  GMM_ASSERT1(false, "Sorry, second derivative not implemented");
454  }
455  };
456 
457 
458  //=================================================================
459  // Von Mises projection
460  //=================================================================
461 
462 
463  struct Von_Mises_projection_operator : public ga_nonlinear_operator {
464  bool result_size(const arg_list &args, bgeot::multi_index &sizes) const {
465  if (args.size() != 2 || args[0]->sizes().size() > 2
466  || args[1]->size() != 1) return false;
467  size_type N = (args[0]->sizes().size() == 2) ? args[0]->sizes()[0] : 1;
468  if (args[0]->sizes().size() == 2 && args[0]->sizes()[1] != N) return false;
469  if (args[0]->sizes().size() != 2 && args[0]->size() != 1) return false;
470  if (N > 1) ga_init_square_matrix(sizes, N); else ga_init_scalar(sizes);
471  return true;
472  }
473 
474  // Value:
475  void value(const arg_list &args, base_tensor &result) const {
476  size_type N = (args[0]->sizes().size() == 2) ? args[0]->sizes()[0] : 1;
477  base_matrix tau(N, N), tau_D(N, N);
478  gmm::copy(args[0]->as_vector(), tau.as_vector());
479  scalar_type s = (*(args[1]))[0];
480 
481  scalar_type tau_m = gmm::mat_trace(tau) / scalar_type(N);
482  gmm::copy(tau, tau_D);
483  for (size_type i = 0; i < N; ++i) tau_D(i,i) -= tau_m;
484 
485  scalar_type norm_tau_D = gmm::mat_euclidean_norm(tau_D);
486 
487  if (norm_tau_D > s) gmm::scale(tau_D, s / norm_tau_D);
488 
489  for (size_type i = 0; i < N; ++i) tau_D(i,i) += tau_m;
490 
491  gmm::copy(tau_D.as_vector(), result.as_vector());
492  }
493 
494  // Derivative:
495  void derivative(const arg_list &args, size_type nder,
496  base_tensor &result) const {
497  size_type N = (args[0]->sizes().size() == 2) ? args[0]->sizes()[0] : 1;
498  base_matrix tau(N, N), tau_D(N, N);
499  gmm::copy(args[0]->as_vector(), tau.as_vector());
500  scalar_type s = (*(args[1]))[0];
501  scalar_type tau_m = gmm::mat_trace(tau) / scalar_type(N);
502  gmm::copy(tau, tau_D);
503  for (size_type i = 0; i < N; ++i) tau_D(i,i) -= tau_m;
504  scalar_type norm_tau_D = gmm::mat_euclidean_norm(tau_D);
505 
506  if (norm_tau_D != scalar_type(0))
507  gmm::scale(tau_D, scalar_type(1)/norm_tau_D);
508 
509  switch(nder) {
510  case 1:
511  if (norm_tau_D <= s) {
512  gmm::clear(result.as_vector());
513  for (size_type i = 0; i < N; ++i)
514  for (size_type j = 0; j < N; ++j)
515  result(i,j,i,j) = scalar_type(1);
516  } else {
517  for (size_type i = 0; i < N; ++i)
518  for (size_type j = 0; j < N; ++j)
519  for (size_type m = 0; m < N; ++m)
520  for (size_type n = 0; n < N; ++n)
521  result(i,j,m,n)
522  = s * (-tau_D(i,j) * tau_D(m,n)
523  + ((i == m && j == n) ? scalar_type(1) : scalar_type(0))
524  - ((i == j && m == n) ? scalar_type(1)/scalar_type(N)
525  : scalar_type(0))) / norm_tau_D;
526  for (size_type i = 0; i < N; ++i)
527  for (size_type j = 0; j < N; ++j)
528  result(i,i,j,j) += scalar_type(1)/scalar_type(N);
529  }
530  break;
531  case 2:
532  if (norm_tau_D < s)
533  gmm::clear(result.as_vector());
534  else
535  gmm::copy(tau_D.as_vector(), result.as_vector());
536  break;
537  }
538  }
539 
540  // Second derivative : not implemented
541  void second_derivative(const arg_list &, size_type, size_type,
542  base_tensor &) const {
543  GMM_ASSERT1(false, "Sorry, second derivative not implemented");
544  }
545  };
546 
547  static bool init_predef_operators(void) {
548 
549  ga_predef_operator_tab &PREDEF_OPERATORS
550  = dal::singleton<ga_predef_operator_tab>::instance();
551 
552  PREDEF_OPERATORS.add_method("Expm",
553  std::make_shared<matrix_exponential_operator>());
554  PREDEF_OPERATORS.add_method("Logm",
555  std::make_shared<matrix_logarithm_operator>());
556  PREDEF_OPERATORS.add_method("Normalized",
557  std::make_shared<normalized_operator>());
558  PREDEF_OPERATORS.add_method("Normalized_reg",
559  std::make_shared<normalized_reg_operator>());
560  PREDEF_OPERATORS.add_method("Ball_projection",
561  std::make_shared<Ball_projection_operator>());
562  PREDEF_OPERATORS.add_method("Von_Mises_projection",
563  std::make_shared<Von_Mises_projection_operator>());
564  return true;
565  }
566 
567  // declared in getfem_generic_assembly.cc
568  bool predef_operators_plasticity_initialized
569  = init_predef_operators();
570 
571 
572  // ======================================
573  //
574  // Small strain plasticity brick
575  //
576  // ======================================
577 
578 
579  // Assembly strings for isotropic perfect elastoplasticity with Von-Mises
580  // criterion (Prandtl-Reuss model). With the use of a plastic multiplier
581  void build_isotropic_perfect_elastoplasticity_expressions_mult
582  (model &md, const std::string &dispname, const std::string &xi,
583  const std::string &Previous_Ep, const std::string &lambda,
584  const std::string &mu, const std::string &sigma_y,
585  const std::string &theta, const std::string &dt,
586  std::string &sigma_np1, std::string &Epnp1, std::string &compcond,
587  std::string &sigma_after, std::string &von_mises) {
588 
589  const mesh_fem *mfu = md.pmesh_fem_of_variable(dispname);
590  size_type N = mfu->linked_mesh().dim();
591  GMM_ASSERT1(mfu && mfu->get_qdim() == N, "The small strain "
592  "elastoplasticity brick can only be applied on a fem "
593  "variable of the same dimension as the mesh");
594 
595  GMM_ASSERT1(!(md.is_data(xi)) && md.pmesh_fem_of_variable(xi),
596  "The provided name '" << xi << "' for the plastic multiplier, "
597  "should be defined as a fem variable");
598 
599  GMM_ASSERT1(md.is_data(Previous_Ep) &&
600  (md.pim_data_of_variable(Previous_Ep) ||
601  md.pmesh_fem_of_variable(Previous_Ep)),
602  "The provided name '" << Previous_Ep << "' for the plastic "
603  "strain tensor at the previous timestep, should be defined "
604  "either as fem or as im data");
605 
606  bgeot::multi_index Ep_size(N, N);
607  GMM_ASSERT1((md.pim_data_of_variable(Previous_Ep) &&
608  md.pim_data_of_variable(Previous_Ep)->tensor_size() == Ep_size)
609  ||
610  (md.pmesh_fem_of_variable(Previous_Ep) &&
611  md.pmesh_fem_of_variable(Previous_Ep)->get_qdims() == Ep_size),
612  "Wrong size of " << Previous_Ep);
613 
614  std::map<std::string, std::string> dict;
615  dict["Grad_u"] = "Grad_"+dispname; dict["xi"] = xi;
616  dict["Previous_xi"] = "Previous_"+xi;
617  dict["Grad_Previous_u"] = "Grad_Previous_"+dispname;
618  dict["theta"] = theta; dict["dt"] = dt; dict["Epn"] = Previous_Ep;
619  dict["lambda"] = lambda; dict["mu"] = mu; dict["sigma_y"] = sigma_y;
620 
621  dict["Enp1"] = ga_substitute("Sym(Grad_u)", dict);
622  dict["En"] = ga_substitute("Sym(Grad_Previous_u)", dict);
623  dict["zetan"] = ga_substitute
624  ("((Epn)+(1-(theta))*(2*(mu)*(dt)*(Previous_xi))*(Deviator(En)-(Epn)))",
625  dict);
626  Epnp1 = ga_substitute
627  ("((zetan)+(1-1/(1+(theta)*2*(mu)*(dt)*(xi)))*(Deviator(Enp1)-(zetan)))",
628  dict);
629  dict["Epnp1"] = Epnp1;
630  sigma_np1 = ga_substitute
631  ("((lambda)*Trace(Enp1)*Id(meshdim)+2*(mu)*((Enp1)-(Epnp1)))", dict);
632  dict["fbound"] = ga_substitute
633  ("(2*(mu)*Norm(Deviator(Enp1)-(Epnp1))-sqrt(2/3)*(sigma_y))", dict);
634  dict["sigma_after"] = sigma_after = ga_substitute
635  ("((lambda)*Trace(Enp1)*Id(meshdim)+2*(mu)*((Enp1)-(Epn)))", dict);
636  compcond = ga_substitute
637  ("((mu)*xi-pos_part((mu)*xi+100*(fbound)/(mu)))", dict);
638  von_mises = ga_substitute
639  ("sqrt(3/2)*Norm(Deviator(sigma_after))", dict);
640  }
641 
642 
643  // Assembly strings for isotropic perfect elastoplasticity with Von-Mises
644  // criterion (Prandtl-Reuss model). Without the use of a plastic multiplier
645  void build_isotropic_perfect_elastoplasticity_expressions_no_mult
646  (model &md, const std::string &dispname, const std::string &xi,
647  const std::string &Previous_Ep, const std::string &lambda,
648  const std::string &mu, const std::string &sigma_y,
649  const std::string &theta, const std::string &dt,
650  std::string &sigma_np1, std::string &Epnp1, std::string &xi_np1,
651  std::string &sigma_after, std::string &von_mises) {
652 
653  const mesh_fem *mfu = md.pmesh_fem_of_variable(dispname);
654  size_type N = mfu->linked_mesh().dim();
655  GMM_ASSERT1(mfu && mfu->get_qdim() == N, "The small strain "
656  "elastoplasticity brick can only be applied on a fem "
657  "variable of the same dimension as the mesh");
658 
659  GMM_ASSERT1(md.is_data(xi) &&
660  (md.pim_data_of_variable(xi) ||
661  md.pmesh_fem_of_variable(xi)),
662  "The provided name '" << xi << "' for the plastic multiplier, "
663  "should be defined either as fem data or as im data");
664 
665  GMM_ASSERT1(md.is_data(Previous_Ep) &&
666  (md.pim_data_of_variable(Previous_Ep) ||
667  md.pmesh_fem_of_variable(Previous_Ep)),
668  "The provided name '" << Previous_Ep << "' for the plastic "
669  "strain tensor at the previous timestep, should be defined "
670  "either as fem or as im data");
671 
672  bgeot::multi_index Ep_size(N, N);
673  GMM_ASSERT1((md.pim_data_of_variable(Previous_Ep) &&
674  md.pim_data_of_variable(Previous_Ep)->tensor_size() == Ep_size)
675  ||
676  (md.pmesh_fem_of_variable(Previous_Ep) &&
677  md.pmesh_fem_of_variable(Previous_Ep)->get_qdims() == Ep_size),
678  "Wrong size of " << Previous_Ep);
679 
680  std::map<std::string, std::string> dict;
681  dict["Grad_u"] = "Grad_"+dispname; dict["xi"] = xi;
682  dict["Previous_xi"] = "Previous_"+xi;
683  dict["Grad_Previous_u"] = "Grad_Previous_"+dispname;
684  dict["theta"] = theta; dict["dt"] = dt; dict["Epn"] = Previous_Ep;
685  dict["lambda"] = lambda; dict["mu"] = mu; dict["sigma_y"] = sigma_y;
686 
687  dict["Enp1"] = ga_substitute("Sym(Grad_u)", dict);
688  dict["En"] = ga_substitute("Sym(Grad_Previous_u)", dict) ;
689 
690 
691  dict["zetan"] = ga_substitute
692  ("(Epn)+(1-(theta))*(2*(mu)*(dt)*(Previous_xi))*(Deviator(En)-(Epn))",
693  dict);
694  dict["B"] = ga_substitute("Deviator(Enp1)-(zetan)", dict);
695  Epnp1 = ga_substitute
696  ("(zetan)+pos_part(1-sqrt(2/3)*(sigma_y)/(2*(mu)*Norm(B)+1e-40))*(B)",
697  dict);
698  dict["Epnp1"] = Epnp1;
699 
700  sigma_np1 = ga_substitute
701  ("(lambda)*Trace(Enp1)*Id(meshdim)+2*(mu)*((Enp1)-(Epnp1))", dict);
702  dict["sigma_after"] = sigma_after = ga_substitute
703  ("(lambda)*Trace(Enp1)*Id(meshdim)+2*(mu)*((Enp1)-(Epn))", dict);
704  xi_np1 = ga_substitute
705  ("pos_part(sqrt(3/2)*Norm(B)/(sigma_y)-1/(2*(mu)))/((theta)*(dt))", dict);
706  von_mises = ga_substitute
707  ("sqrt(3/2)*Norm(Deviator(sigma_after))", dict);
708  }
709 
710  // Assembly strings for isotropic perfect elastoplasticity with Von-Mises
711  // criterion (Prandtl-Reuss model). With the use of a plastic multiplier
712  // and plane strain version.
713  void build_isotropic_perfect_elastoplasticity_expressions_mult_ps
714  (model &md, const std::string &dispname, const std::string &xi,
715  const std::string &Previous_Ep, const std::string &lambda,
716  const std::string &mu, const std::string &sigma_y,
717  const std::string &theta, const std::string &dt,
718  std::string &sigma_np1, std::string &Epnp1, std::string &compcond,
719  std::string &sigma_after, std::string &von_mises) {
720 
721  const mesh_fem *mfu = md.pmesh_fem_of_variable(dispname);
722  size_type N = mfu->linked_mesh().dim();
723  GMM_ASSERT1(N == 2, "This plastic law is restricted to 2D");
724 
725  GMM_ASSERT1(mfu && mfu->get_qdim() == N, "The small strain "
726  "elastoplasticity brick can only be applied on a fem "
727  "variable of the same dimension as the mesh");
728 
729  GMM_ASSERT1(!(md.is_data(xi)) && md.pmesh_fem_of_variable(xi),
730  "The provided name '" << xi << "' for the plastic multiplier, "
731  "should be defined as a fem variable");
732 
733  GMM_ASSERT1(md.is_data(Previous_Ep) &&
734  (md.pim_data_of_variable(Previous_Ep) ||
735  md.pmesh_fem_of_variable(Previous_Ep)),
736  "The provided name '" << Previous_Ep << "' for the plastic "
737  "strain tensor at the previous timestep, should be defined "
738  "either as fem or as im data");
739 
740  bgeot::multi_index Ep_size(N, N);
741  GMM_ASSERT1((md.pim_data_of_variable(Previous_Ep) &&
742  md.pim_data_of_variable(Previous_Ep)->tensor_size() == Ep_size)
743  ||
744  (md.pmesh_fem_of_variable(Previous_Ep) &&
745  md.pmesh_fem_of_variable(Previous_Ep)->get_qdims() == Ep_size),
746  "Wrong size of " << Previous_Ep);
747 
748  std::map<std::string, std::string> dict;
749  dict["Grad_u"] = "Grad_"+dispname; dict["xi"] = xi;
750  dict["Previous_xi"] = "Previous_"+xi;
751  dict["Grad_Previous_u"] = "Grad_Previous_"+dispname;
752  dict["theta"] = theta; dict["dt"] = dt; dict["Epn"] = Previous_Ep;
753  dict["lambda"] = lambda; dict["mu"] = mu; dict["sigma_y"] = sigma_y;
754 
755  dict["Enp1"] = ga_substitute("Sym(Grad_u)", dict);
756  dict["En"] = ga_substitute("Sym(Grad_Previous_u)", dict);
757  dict["Dev_En"]= ga_substitute("(En-(Trace(En)/3)*Id(meshdim))", dict);
758  dict["Dev_Enp1"]= ga_substitute("(Enp1-(Trace(Enp1)/3)*Id(meshdim))", dict);
759  dict["zetan"] = ga_substitute
760  ("((Epn)+(1-(theta))*(2*(mu)*(dt)*(Previous_xi))*((Dev_En)-(Epn)))",
761  dict);
762  Epnp1 = ga_substitute
763  ("((zetan)+(1-1/(1+(theta)*2*(mu)*(dt)*(xi)))*((Dev_Enp1)-(zetan)))",
764  dict);
765  dict["Epnp1"] = Epnp1;
766  sigma_np1 = ga_substitute
767  ("((lambda)*Trace(Enp1)*Id(meshdim)+2*(mu)*((Enp1)-(Epnp1)))", dict);
768  dict["fbound"] = ga_substitute
769  ("(2*(mu)*sqrt(Norm_sqr(Dev_Enp1-(Epnp1))"
770  "+sqr(Trace(Enp1)/3-Trace(Epnp1)))-sqrt(2/3)*(sigma_y))", dict);
771 
772  sigma_after = ga_substitute
773  ("((lambda)*Trace(Enp1)*Id(meshdim)+2*(mu)*((Enp1)-(Epn)))", dict);
774  compcond = ga_substitute
775  ("((mu)*xi-pos_part((mu)*xi+100*(fbound)/(mu)))", dict);
776  von_mises = ga_substitute
777  ("sqrt(3/2)*sqrt(Norm_sqr((2*(mu))*(Dev_En)-(2*(mu))*(Epn))"
778  "+sqr(2*(mu)*Trace(En)/3-(2*(mu))*Trace(Epn)))", dict);
779  }
780 
781 
782  // Assembly strings for isotropic perfect elastoplasticity with Von-Mises
783  // criterion (Prandtl-Reuss model). Without the use of a plastic multiplier
784  // and plane strain version.
785  void build_isotropic_perfect_elastoplasticity_expressions_no_mult_ps
786  (model &md, const std::string &dispname, const std::string &xi,
787  const std::string &Previous_Ep, const std::string &lambda,
788  const std::string &mu, const std::string &sigma_y,
789  const std::string &theta, const std::string &dt,
790  std::string &sigma_np1, std::string &Epnp1,
791  std::string &xi_np1, std::string &sigma_after, std::string &von_mises) {
792 
793  const mesh_fem *mfu = md.pmesh_fem_of_variable(dispname);
794  size_type N = mfu->linked_mesh().dim();
795  GMM_ASSERT1(N == 2, "This plastic law is restricted to 2D");
796 
797  GMM_ASSERT1(mfu && mfu->get_qdim() == N, "The small strain "
798  "elastoplasticity brick can only be applied on a fem "
799  "variable of the same dimension as the mesh");
800 
801  GMM_ASSERT1(md.is_data(xi) &&
802  (md.pim_data_of_variable(xi) ||
803  md.pmesh_fem_of_variable(xi)),
804  "The provided name '" << xi << "' for the plastic multiplier, "
805  "should be defined either as fem data or as im data");
806 
807  GMM_ASSERT1(md.is_data(Previous_Ep) &&
808  (md.pim_data_of_variable(Previous_Ep) ||
809  md.pmesh_fem_of_variable(Previous_Ep)),
810  "The provided name '" << Previous_Ep << "' for the plastic "
811  "strain tensor at the previous timestep, should be defined "
812  "either as fem or as im data");
813 
814  bgeot::multi_index Ep_size(N, N);
815  GMM_ASSERT1((md.pim_data_of_variable(Previous_Ep) &&
816  md.pim_data_of_variable(Previous_Ep)->tensor_size() == Ep_size)
817  ||
818  (md.pmesh_fem_of_variable(Previous_Ep) &&
819  md.pmesh_fem_of_variable(Previous_Ep)->get_qdims() == Ep_size),
820  "Wrong size of " << Previous_Ep);
821 
822  std::map<std::string, std::string> dict;
823  dict["Grad_u"] = "Grad_"+dispname; dict["xi"] = xi;
824  dict["Previous_xi"] = "Previous_"+xi;
825  dict["Grad_Previous_u"] = "Grad_Previous_"+dispname;
826  dict["theta"] = theta; dict["dt"] = dt; dict["Epn"] = Previous_Ep;
827  dict["lambda"] = lambda; dict["mu"] = mu; dict["sigma_y"] = sigma_y;
828 
829  dict["Enp1"] = ga_substitute("Sym(Grad_u)", dict);
830  dict["En"] = ga_substitute("Sym(Grad_Previous_u)", dict) ;
831  dict["Dev_En"]= ga_substitute("(En-(Trace(En)/3)*Id(meshdim))", dict);
832  dict["Dev_Enp1"]= ga_substitute("(Enp1-(Trace(Enp1)/3)*Id(meshdim))", dict);
833  dict["zetan"] = ga_substitute
834  ("((Epn)+(1-(theta))*(2*(mu)*(dt)*(Previous_xi))*(Dev_En-(Epn)))",
835  dict);
836  dict["B"] = ga_substitute("(Dev_Enp1)-(zetan)", dict);
837  Epnp1 = ga_substitute
838  ("(zetan)+pos_part(1-sqrt(2/3)*(sigma_y)/(2*(mu)*(sqrt(Norm_sqr(B)+"
839  "sqr(Trace(Enp1)/3-Trace(zetan))))+1e-25))*(B)", dict);
840  dict["Epnp1"] = Epnp1;
841 
842  sigma_np1 = ga_substitute
843  ("(lambda)*Trace(Enp1)*Id(meshdim)+2*(mu)*((Enp1)-(Epnp1))", dict);
844  sigma_after = ga_substitute
845  ("(lambda)*Trace(Enp1)*Id(meshdim)+2*(mu)*((Enp1)-(Epn))", dict);
846  xi_np1 = ga_substitute
847  ("pos_part(sqrt(3/2)*Norm(B)/(sigma_y)-1/(2*(mu)))/((theta)*(dt))", dict);
848  von_mises = ga_substitute
849  ("sqrt(3/2)*sqrt(Norm_sqr((2*(mu))*(Dev_En)-(2*(mu))*(Epn))"
850  "+sqr(2*(mu)*Trace(En)/3-(2*(mu))*Trace(Epn)))", dict);
851  }
852 
853  // Assembly strings for isotropic elastoplasticity with Von-Mises
854  // criterion (Prandtl-Reuss model) and linear isotropic and kinematic
855  // hardening. With the use of a plastic multiplier
856  void build_isotropic_perfect_elastoplasticity_expressions_hard_mult
857  (model &md, const std::string &dispname, const std::string &xi,
858  const std::string &Previous_Ep, const std::string &alpha,
859  const std::string &lambda, const std::string &mu,
860  const std::string &sigma_y, const std::string &Hk, const std::string &Hi,
861  const std::string &theta, const std::string &dt,
862  std::string &sigma_np1, std::string &Epnp1, std::string &compcond,
863  std::string &sigma_after, std::string &von_mises, std::string &alphanp1) {
864 
865  const mesh_fem *mfu = md.pmesh_fem_of_variable(dispname);
866  size_type N = mfu->linked_mesh().dim();
867  GMM_ASSERT1(mfu && mfu->get_qdim() == N, "The small strain "
868  "elastoplasticity brick can only be applied on a fem "
869  "variable of the same dimension as the mesh");
870 
871  GMM_ASSERT1(!(md.is_data(xi)) && md.pmesh_fem_of_variable(xi),
872  "The provided name '" << xi << "' for the plastic multiplier, "
873  "should be defined as a fem variable");
874 
875  GMM_ASSERT1(md.is_data(Previous_Ep) &&
876  (md.pim_data_of_variable(Previous_Ep) ||
877  md.pmesh_fem_of_variable(Previous_Ep)),
878  "The provided name '" << Previous_Ep << "' for the plastic "
879  "strain tensor at the previous timestep, should be defined "
880  "either as fem or as im data");
881 
882  bgeot::multi_index Ep_size(N, N);
883  GMM_ASSERT1((md.pim_data_of_variable(Previous_Ep) &&
884  md.pim_data_of_variable(Previous_Ep)->tensor_size() == Ep_size)
885  ||
886  (md.pmesh_fem_of_variable(Previous_Ep) &&
887  md.pmesh_fem_of_variable(Previous_Ep)->get_qdims() == Ep_size),
888  "Wrong size of " << Previous_Ep);
889 
890  std::map<std::string, std::string> dict;
891  dict["Hk"] = Hk; dict["Hi"] = Hi; dict["alphan"] = alpha;
892  dict["Grad_u"] = "Grad_"+dispname; dict["xi"] = xi;
893  dict["Previous_xi"] = "Previous_"+xi;
894  dict["Grad_Previous_u"] = "Grad_Previous_"+dispname;
895  dict["theta"] = theta; dict["dt"] = dt; dict["Epn"] = Previous_Ep;
896  dict["lambda"] = lambda; dict["mu"] = mu; dict["sigma_y"] = sigma_y;
897 
898  dict["Enp1"] = ga_substitute("Sym(Grad_u)", dict);
899  dict["En"] = ga_substitute("Sym(Grad_Previous_u)", dict);
900  dict["zetan"] = ga_substitute
901  ("((Epn)+(1-(theta))*((dt)*(Previous_xi))*((2*(mu))*Deviator(En)"
902  "-(2*(mu)+2/3*(Hk))*(Epn)))", dict);
903  dict["etan"] = ga_substitute
904  ("((alphan)+sqrt(2/3)*(1-(theta))*((dt)*(Previous_xi))*"
905  "Norm((2*(mu))*Deviator(En)-(2*(mu)+2/3*(Hk))*(Epn)))", dict);
906  dict["B"] = ga_substitute
907  ("((2*(mu))*Deviator(Enp1)-(2*(mu)+2/3*(Hk))*(zetan))", dict);
908  dict["beta"] = ga_substitute
909  ("((theta)*(dt)*(xi)/(1+(2*(mu)+2/3*(Hk))*(theta)*(dt)*(xi)))", dict);
910  dict["Epnp1"] = Epnp1 = ga_substitute("((zetan)+(beta)*(B))", dict);
911  alphanp1 = ga_substitute("((etan)+sqrt(2/3)*(beta)*Norm(B))", dict);
912  dict["alphanp1"] = alphanp1;
913  sigma_np1 = ga_substitute
914  ("((lambda)*Trace(Enp1)*Id(meshdim)+2*(mu)*((Enp1)-(Epnp1)))", dict);
915  dict["fbound"] = ga_substitute
916  ("(Norm((2*(mu))*Deviator(Enp1)-(2*(mu)+2/3*(Hk))*(Epnp1))"
917  "-sqrt(2/3)*(sigma_y+(Hi)*(alphanp1)))", dict);
918  dict["sigma_after"] = sigma_after = ga_substitute
919  ("((lambda)*Trace(Enp1)*Id(meshdim)+2*(mu)*((Enp1)-(Epn)))", dict);
920  compcond = ga_substitute
921  ("((mu)*xi-pos_part((mu)*xi+100*(fbound)/(mu)))", dict);
922  von_mises = ga_substitute
923  ("sqrt(3/2)*Norm(Deviator(sigma_after))", dict);
924  }
925 
926  // Assembly strings for isotropic elastoplasticity with Von-Mises
927  // criterion (Prandtl-Reuss model) and linear isotropic and kinematic
928  // hardening. Without the use of a plastic multiplier
929  void build_isotropic_perfect_elastoplasticity_expressions_hard_no_mult
930  (model &md, const std::string &dispname, const std::string &xi,
931  const std::string &Previous_Ep, const std::string &alpha,
932  const std::string &lambda, const std::string &mu,
933  const std::string &sigma_y, const std::string &Hk, const std::string &Hi,
934  const std::string &theta, const std::string &dt,
935  std::string &sigma_np1, std::string &Epnp1, std::string &xi_np1,
936  std::string &sigma_after, std::string &von_mises, std::string &alphanp1) {
937 
938  const mesh_fem *mfu = md.pmesh_fem_of_variable(dispname);
939  size_type N = mfu->linked_mesh().dim();
940  GMM_ASSERT1(mfu && mfu->get_qdim() == N, "The small strain "
941  "elastoplasticity brick can only be applied on a fem "
942  "variable of the same dimension as the mesh");
943 
944  GMM_ASSERT1(md.is_data(xi) &&
945  (md.pim_data_of_variable(xi) ||
946  md.pmesh_fem_of_variable(xi)),
947  "The provided name '" << xi << "' for the plastic multiplier, "
948  "should be defined either as fem data or as im data");
949 
950  GMM_ASSERT1(md.is_data(Previous_Ep) &&
951  (md.pim_data_of_variable(Previous_Ep) ||
952  md.pmesh_fem_of_variable(Previous_Ep)),
953  "The provided name '" << Previous_Ep << "' for the plastic "
954  "strain tensor at the previous timestep, should be defined "
955  "either as fem or as im data");
956 
957  bgeot::multi_index Ep_size(N, N);
958  GMM_ASSERT1((md.pim_data_of_variable(Previous_Ep) &&
959  md.pim_data_of_variable(Previous_Ep)->tensor_size() == Ep_size)
960  ||
961  (md.pmesh_fem_of_variable(Previous_Ep) &&
962  md.pmesh_fem_of_variable(Previous_Ep)->get_qdims() == Ep_size),
963  "Wrong size of " << Previous_Ep);
964 
965  std::map<std::string, std::string> dict;
966  dict["Hk"] = Hk; dict["Hi"] = Hi; dict["alphan"] = alpha;
967  dict["Grad_u"] = "Grad_"+dispname; dict["xi"] = xi;
968  dict["Previous_xi"] = "Previous_"+xi;
969  dict["Grad_Previous_u"] = "Grad_Previous_"+dispname;
970  dict["theta"] = theta; dict["dt"] = dt; dict["Epn"] = Previous_Ep;
971  dict["lambda"] = lambda; dict["mu"] = mu; dict["sigma_y"] = sigma_y;
972 
973  dict["Enp1"] = ga_substitute("Sym(Grad_u)", dict);
974  dict["En"] = ga_substitute("Sym(Grad_Previous_u)", dict);
975  dict["zetan"] = ga_substitute
976  ("((Epn)+(1-(theta))*((dt)*(Previous_xi))*((2*(mu))*Deviator(En)"
977  "-(2*(mu)+2/3*(Hk))*(Epn)))", dict);
978  dict["etan"] = ga_substitute
979  ("((alphan)+sqrt(2/3)*(1-(theta))*((dt)*(Previous_xi))*"
980  "Norm((2*(mu))*Deviator(En)-(2*(mu)+2/3*(Hk))*(Epn)))", dict);
981 
982  dict["B"] = ga_substitute
983  ("((2*(mu))*Deviator(Enp1)-(2*(mu)+2/3*(Hk))*(zetan))", dict);
984  dict["beta"] =
985  ga_substitute("(1/((Norm(B)+1e-40)*(2*(mu)+2/3*(Hk)+(2/3)*(Hi))))"
986  "*pos_part(Norm(B)-sqrt(2/3)*((sigma_y)+(Hi)*(etan)))", dict);
987  dict["Epnp1"] = Epnp1 = ga_substitute("((zetan)+(beta)*(B))", dict);
988  alphanp1 = ga_substitute("((etan)+sqrt(2/3)*(beta)*Norm(B))", dict);
989  dict["alphanp1"] = alphanp1;
990 
991  sigma_np1 = ga_substitute
992  ("(lambda)*Trace(Enp1)*Id(meshdim)+2*(mu)*((Enp1)-(Epnp1))", dict);
993  dict["sigma_after"] = sigma_after = ga_substitute
994  ("(lambda)*Trace(Enp1)*Id(meshdim)+2*(mu)*((Enp1)-(Epn))", dict);
995  xi_np1 = ga_substitute
996  ("(((beta)/(1-(2*(mu)+2/3*(Hk))*(beta)))/((theta)*(dt)))", dict);
997  von_mises = ga_substitute
998  ("sqrt(3/2)*Norm(Deviator(sigma_after))", dict);
999  }
1000 
1001  // Assembly strings for isotropic elastoplasticity with Von-Mises
1002  // criterion (Prandtl-Reuss model) and linear isotropic and kinematic
1003  // hardening. With the use of a plastic multiplier and plane strain version.
1004  void build_isotropic_perfect_elastoplasticity_expressions_hard_mult_ps
1005  (model &md, const std::string &dispname, const std::string &xi,
1006  const std::string &Previous_Ep, const std::string &alpha,
1007  const std::string &lambda, const std::string &mu,
1008  const std::string &sigma_y, const std::string &Hk, const std::string &Hi,
1009  const std::string &theta, const std::string &dt,
1010  std::string &sigma_np1, std::string &Epnp1, std::string &compcond,
1011  std::string &sigma_after, std::string &von_mises, std::string &alphanp1) {
1012 
1013  const mesh_fem *mfu = md.pmesh_fem_of_variable(dispname);
1014  size_type N = mfu->linked_mesh().dim();
1015  GMM_ASSERT1(N == 2, "This plastic law is restricted to 2D");
1016 
1017  GMM_ASSERT1(mfu && mfu->get_qdim() == N, "The small strain "
1018  "elastoplasticity brick can only be applied on a fem "
1019  "variable of the same dimension as the mesh");
1020 
1021  GMM_ASSERT1(!(md.is_data(xi)) && md.pmesh_fem_of_variable(xi),
1022  "The provided name '" << xi << "' for the plastic multiplier, "
1023  "should be defined as a fem variable");
1024 
1025  GMM_ASSERT1(md.is_data(Previous_Ep) &&
1026  (md.pim_data_of_variable(Previous_Ep) ||
1027  md.pmesh_fem_of_variable(Previous_Ep)),
1028  "The provided name '" << Previous_Ep << "' for the plastic "
1029  "strain tensor at the previous timestep, should be defined "
1030  "either as fem or as im data");
1031 
1032  bgeot::multi_index Ep_size(N, N);
1033  GMM_ASSERT1((md.pim_data_of_variable(Previous_Ep) &&
1034  md.pim_data_of_variable(Previous_Ep)->tensor_size() == Ep_size)
1035  ||
1036  (md.pmesh_fem_of_variable(Previous_Ep) &&
1037  md.pmesh_fem_of_variable(Previous_Ep)->get_qdims() == Ep_size),
1038  "Wrong size of " << Previous_Ep);
1039 
1040  std::map<std::string, std::string> dict;
1041  dict["Hk"] = Hk; dict["Hi"] = Hi; dict["alphan"] = alpha;
1042  dict["Grad_u"] = "Grad_"+dispname; dict["xi"] = xi;
1043  dict["Previous_xi"] = "Previous_"+xi;
1044  dict["Grad_Previous_u"] = "Grad_Previous_"+dispname;
1045  dict["theta"] = theta; dict["dt"] = dt; dict["Epn"] = Previous_Ep;
1046  dict["lambda"] = lambda; dict["mu"] = mu; dict["sigma_y"] = sigma_y;
1047 
1048  dict["Enp1"] = ga_substitute("Sym(Grad_u)", dict);
1049  dict["En"] = ga_substitute("Sym(Grad_Previous_u)", dict);
1050  dict["Dev_En"]= ga_substitute("(En-(Trace(En)/3)*Id(meshdim))", dict);
1051  dict["Dev_Enp1"]= ga_substitute("(Enp1-(Trace(Enp1)/3)*Id(meshdim))", dict);
1052  dict["zetan"] = ga_substitute
1053  ("((Epn)+(1-(theta))*((dt)*(Previous_xi))*((2*(mu))*(Dev_En)"
1054  "-(2*(mu)+2/3*(Hk))*(Epn)))", dict);
1055  dict["etan"] = ga_substitute
1056  ("((alphan)+sqrt(2/3)*(1-(theta))*((dt)*(Previous_xi))*"
1057  "sqrt(Norm_sqr((2*(mu))*(Dev_En)-(2*(mu)+2/3*(Hk))*(Epn))"
1058  "+sqr(2*(mu)*Trace(En)/3-(2*(mu)+2/3*(Hk))*Trace(Epn))))", dict);
1059  dict["B"] = ga_substitute("((2*(mu))*(Dev_Enp1)-(2*(mu)+2/3*(Hk))*(zetan))",
1060  dict);
1061  dict["Norm_B"] = ga_substitute("sqrt(Norm_sqr(B)+sqr(2*(mu)*Trace(Enp1)/3"
1062  "-(2*(mu)+2/3*(Hk))*Trace(zetan)))", dict);
1063 
1064  dict["beta"] = ga_substitute
1065  ("((theta)*(dt)*(xi)/(1+(2*(mu)+2/3*(Hk))*(theta)*(dt)*(xi)))", dict);
1066  dict["Epnp1"] = Epnp1 = ga_substitute("((zetan)+(beta)*(B))", dict);
1067  alphanp1 = ga_substitute("((etan)+sqrt(2/3)*(beta)*(Norm_B))", dict);
1068  dict["alphanp1"] = alphanp1;
1069  sigma_np1 = ga_substitute
1070  ("((lambda)*Trace(Enp1)*Id(meshdim)+2*(mu)*((Enp1)-(Epnp1)))", dict);
1071  dict["fbound"] = ga_substitute
1072  ("(sqrt(Norm_sqr((2*(mu))*(Dev_Enp1)-(2*(mu)+2/3*(Hk))*(Epnp1))"
1073  "+sqr(2*(mu)*Trace(Enp1)/3-(2*(mu)+2/3*(Hk))*Trace(Epnp1)))"
1074  "-sqrt(2/3)*(sigma_y+(Hi)*(alphanp1)))", dict);
1075  sigma_after = ga_substitute
1076  ("((lambda)*Trace(Enp1)*Id(meshdim)+2*(mu)*((Enp1)-(Epn)))", dict);
1077  compcond = ga_substitute
1078  ("((mu)*xi-pos_part((mu)*xi+100*(fbound)/(mu)))", dict);
1079  von_mises = ga_substitute
1080  ("sqrt(3/2)*sqrt(Norm_sqr((2*(mu))*(Dev_En)-(2*(mu)+2/3*(Hk))*(Epn))"
1081  "+sqr(2*(mu)*Trace(En)/3-(2*(mu)+2/3*(Hk))*Trace(Epn)))", dict);
1082  }
1083 
1084  // Assembly strings for isotropic elastoplasticity with Von-Mises
1085  // criterion (Prandtl-Reuss model) and linear isotropic and kinematic
1086  // hardening.
1087  // Without the use of a plastic multiplier and plane strain version.
1088  void build_isotropic_perfect_elastoplasticity_expressions_hard_no_mult_ps
1089  (model &md, const std::string &dispname, const std::string &xi,
1090  const std::string &Previous_Ep, const std::string &alpha,
1091  const std::string &lambda, const std::string &mu,
1092  const std::string &sigma_y, const std::string &Hk, const std::string &Hi,
1093  const std::string &theta, const std::string &dt,
1094  std::string &sigma_np1, std::string &Epnp1, std::string &xi_np1,
1095  std::string &sigma_after, std::string &von_mises, std::string &alphanp1) {
1096 
1097  const mesh_fem *mfu = md.pmesh_fem_of_variable(dispname);
1098  size_type N = mfu->linked_mesh().dim();
1099  GMM_ASSERT1(N == 2, "This plastic law is restricted to 2D");
1100 
1101  GMM_ASSERT1(mfu && mfu->get_qdim() == N, "The small strain "
1102  "elastoplasticity brick can only be applied on a fem "
1103  "variable of the same dimension as the mesh");
1104 
1105  GMM_ASSERT1(md.is_data(xi) &&
1106  (md.pim_data_of_variable(xi) ||
1107  md.pmesh_fem_of_variable(xi)),
1108  "The provided name '" << xi << "' for the plastic multiplier, "
1109  "should be defined either as fem data or as im data");
1110 
1111  GMM_ASSERT1(md.is_data(Previous_Ep) &&
1112  (md.pim_data_of_variable(Previous_Ep) ||
1113  md.pmesh_fem_of_variable(Previous_Ep)),
1114  "The provided name '" << Previous_Ep << "' for the plastic "
1115  "strain tensor at the previous timestep, should be defined "
1116  "either as fem or as im data");
1117 
1118  bgeot::multi_index Ep_size(N, N);
1119  GMM_ASSERT1((md.pim_data_of_variable(Previous_Ep) &&
1120  md.pim_data_of_variable(Previous_Ep)->tensor_size() == Ep_size)
1121  ||
1122  (md.pmesh_fem_of_variable(Previous_Ep) &&
1123  md.pmesh_fem_of_variable(Previous_Ep)->get_qdims() == Ep_size),
1124  "Wrong size of " << Previous_Ep);
1125 
1126  std::map<std::string, std::string> dict;
1127  dict["Hk"] = Hk; dict["Hi"] = Hi; dict["alphan"] = alpha;
1128  dict["Grad_u"] = "Grad_"+dispname; dict["xi"] = xi;
1129  dict["Previous_xi"] = "Previous_"+xi;
1130  dict["Grad_Previous_u"] = "Grad_Previous_"+dispname;
1131  dict["theta"] = theta; dict["dt"] = dt; dict["Epn"] = Previous_Ep;
1132  dict["lambda"] = lambda; dict["mu"] = mu; dict["sigma_y"] = sigma_y;
1133 
1134  dict["Enp1"] = ga_substitute("Sym(Grad_u)", dict);
1135  dict["En"] = ga_substitute("Sym(Grad_Previous_u)", dict);
1136  dict["Dev_En"]= ga_substitute("(En-(Trace(En)/3)*Id(meshdim))", dict);
1137  dict["Dev_Enp1"]= ga_substitute("(Enp1-(Trace(Enp1)/3)*Id(meshdim))", dict);
1138  dict["zetan"] = ga_substitute
1139  ("((Epn)+(1-(theta))*((dt)*(Previous_xi))*((2*(mu))*(Dev_En)"
1140  "-(2*(mu)+2/3*(Hk))*(Epn)))", dict);
1141  dict["etan"] = ga_substitute
1142  ("((alphan)+sqrt(2/3)*(1-(theta))*((dt)*(Previous_xi))*"
1143  "sqrt(Norm_sqr((2*(mu))*(Dev_En)-(2*(mu)+2/3*(Hk))*(Epn))"
1144  "+sqr(2*(mu)*Trace(En)/3-(2*(mu)+2/3*(Hk))*Trace(Epn))))", dict);
1145 
1146  dict["B"] = ga_substitute
1147  ("((2*(mu))*(Dev_Enp1)-(2*(mu)+2/3*(Hk))*(zetan))", dict);
1148  dict["Norm_B"] = ga_substitute("sqrt(Norm_sqr(B)+sqr(2*(mu)*Trace(Enp1)/3"
1149  "-(2*(mu)+2/3*(Hk))*Trace(zetan)))", dict);
1150  dict["beta"] =
1151  ga_substitute("(1/(((Norm_B)+1e-40)*(2*(mu)+2/3*(Hk)+(2/3)*(Hi))))"
1152  "*pos_part((Norm_B)-sqrt(2/3)*((sigma_y)+(Hi)*(etan)))", dict);
1153  dict["Epnp1"] = Epnp1 = ga_substitute("((zetan)+(beta)*(B))", dict);
1154  alphanp1 = ga_substitute("((etan)+sqrt(2/3)*(beta)*(Norm_B))", dict);
1155  dict["alphanp1"] = alphanp1;
1156 
1157  sigma_np1 = ga_substitute
1158  ("(lambda)*Trace(Enp1)*Id(meshdim)+2*(mu)*((Enp1)-(Epnp1))", dict);
1159  sigma_after = ga_substitute
1160  ("(lambda)*Trace(Enp1)*Id(meshdim)+2*(mu)*((Enp1)-(Epn))", dict);
1161  xi_np1 = ga_substitute
1162  ("(((beta)/(1-(2*(mu)+2/3*(Hk))*(beta)))/((theta)*(dt)))", dict);
1163  von_mises = ga_substitute
1164  ("sqrt(3/2)*sqrt(Norm_sqr((2*(mu))*(Dev_En)-(2*(mu)+2/3*(Hk))*(Epn))"
1165  "+sqr(2*(mu)*Trace(En)/3-(2*(mu)+2/3*(Hk))*Trace(Epn)))", dict);
1166  }
1167 
1168  void build_isotropic_perfect_elastoplasticity_expressions_generic
1169  (model &md, const std::string &lawname,
1170  plasticity_unknowns_type unknowns_type,
1171  const std::vector<std::string> &varnames,
1172  const std::vector<std::string> &params,
1173  std::string &sigma_np1, std::string &Epnp1,
1174  std::string &compcond, std::string &xi_np1,
1175  std::string &sigma_after, std::string &von_mises,
1176  std::string &alphanp1) {
1177 
1178  GMM_ASSERT1(unknowns_type == DISPLACEMENT_ONLY ||
1179  unknowns_type == DISPLACEMENT_AND_PLASTIC_MULTIPLIER,
1180  "Not supported type of unknowns");
1181  bool plastic_multiplier_is_var
1182  = (unknowns_type == DISPLACEMENT_AND_PLASTIC_MULTIPLIER);
1183 
1184  bool hardening = (lawname.find("_hardening") != std::string::npos);
1185  size_type nhard = hardening ? 2 : 0;
1186 
1187  GMM_ASSERT1(varnames.size() == (hardening ? 4 : 3),
1188  "Incorrect number of variables: " << varnames.size());
1189  GMM_ASSERT1(params.size() >= 3+nhard &&
1190  params.size() <= 5+nhard,
1191  "Incorrect number of parameters: " << params.size());
1192  const std::string &dispname = sup_previous_and_dot_to_varname(varnames[0]);
1193  const std::string &xi = sup_previous_and_dot_to_varname(varnames[1]);
1194  const std::string &Previous_Ep = varnames[2];
1195  const std::string &alpha = hardening ? varnames[3] : "";
1196  const std::string &lambda = params[0];
1197  const std::string &mu = params[1];
1198  const std::string &sigma_y = params[2];
1199  const std::string &Hk = hardening ? params[3] : "";
1200  const std::string &Hi = hardening ? params[4] : "";
1201  const std::string &theta = (params.size() >= 4+nhard)
1202  ? params[3+nhard] : "1";
1203  const std::string &dt = (params.size() >= 5+nhard)
1204  ? params[4+nhard] : "timestep";
1205 
1206  sigma_np1 = Epnp1 = compcond = xi_np1 = "";;
1207  sigma_after = von_mises = alphanp1 = "";
1208 
1209  if (lawname.compare("isotropic_perfect_plasticity") == 0 ||
1210  lawname.compare("prandtl_reuss") == 0) {
1211  if (plastic_multiplier_is_var) {
1212  build_isotropic_perfect_elastoplasticity_expressions_mult
1213  (md, dispname, xi, Previous_Ep, lambda, mu, sigma_y, theta, dt,
1214  sigma_np1, Epnp1, compcond, sigma_after, von_mises);
1215  } else {
1216  build_isotropic_perfect_elastoplasticity_expressions_no_mult
1217  (md, dispname, xi, Previous_Ep, lambda, mu, sigma_y, theta, dt,
1218  sigma_np1, Epnp1, xi_np1, sigma_after, von_mises);
1219  }
1220  } else if (lawname.compare("plane_strain_isotropic_perfect_plasticity")
1221  == 0 ||
1222  lawname.compare("plane_strain_prandtl_reuss") == 0) {
1223  if (plastic_multiplier_is_var) {
1224  build_isotropic_perfect_elastoplasticity_expressions_mult_ps
1225  (md, dispname, xi, Previous_Ep, lambda, mu, sigma_y, theta, dt,
1226  sigma_np1, Epnp1, compcond, sigma_after, von_mises);
1227  } else {
1228  build_isotropic_perfect_elastoplasticity_expressions_no_mult_ps
1229  (md, dispname, xi, Previous_Ep, lambda, mu, sigma_y, theta, dt,
1230  sigma_np1, Epnp1, xi_np1, sigma_after, von_mises);
1231  }
1232  } else if (lawname.compare("isotropic_plasticity_linear_hardening") == 0 ||
1233  lawname.compare("prandtl_reuss_linear_hardening") == 0) {
1234  if (plastic_multiplier_is_var) {
1235  build_isotropic_perfect_elastoplasticity_expressions_hard_mult
1236  (md, dispname, xi, Previous_Ep, alpha,
1237  lambda, mu, sigma_y, Hk, Hi, theta, dt,
1238  sigma_np1, Epnp1, compcond, sigma_after, von_mises, alphanp1);
1239  } else {
1240  build_isotropic_perfect_elastoplasticity_expressions_hard_no_mult
1241  (md, dispname, xi, Previous_Ep, alpha,
1242  lambda, mu, sigma_y, Hk, Hi, theta, dt,
1243  sigma_np1, Epnp1, xi_np1, sigma_after, von_mises, alphanp1);
1244  }
1245  } else if
1246  (lawname.compare("plane_strain_isotropic_plasticity_linear_hardening")
1247  == 0 ||
1248  lawname.compare("plane_strain_prandtl_reuss_linear_hardening") == 0) {
1249  if (plastic_multiplier_is_var) {
1250  build_isotropic_perfect_elastoplasticity_expressions_hard_mult_ps
1251  (md, dispname, xi, Previous_Ep, alpha,
1252  lambda, mu, sigma_y, Hk, Hi, theta, dt,
1253  sigma_np1, Epnp1, compcond, sigma_after, von_mises, alphanp1);
1254  } else {
1255  build_isotropic_perfect_elastoplasticity_expressions_hard_no_mult_ps
1256  (md, dispname, xi, Previous_Ep, alpha,
1257  lambda, mu, sigma_y, Hk, Hi, theta, dt,
1258  sigma_np1, Epnp1, xi_np1, sigma_after, von_mises, alphanp1);
1259  }
1260  } else
1261  GMM_ASSERT1(false, lawname << " is not an implemented elastoplastic law");
1262  }
1263 
1264  static void filter_lawname(std::string &lawname) {
1265  for (auto &c : lawname)
1266  { if (c == ' ') c = '_'; if (c >= 'A' && c <= 'Z') c = char(c+'a'-'A'); }
1267  }
1268 
1269  size_type add_small_strain_elastoplasticity_brick
1270  (model &md, const mesh_im &mim,
1271  std::string lawname, plasticity_unknowns_type unknowns_type,
1272  const std::vector<std::string> &varnames,
1273  const std::vector<std::string> &params, size_type region) {
1274 
1275  filter_lawname(lawname);
1276  std::string sigma_np1, compcond;
1277  {
1278  std::string dum2, dum4, dum5, dum6, dum7;
1279  build_isotropic_perfect_elastoplasticity_expressions_generic
1280  (md, lawname, unknowns_type, varnames, params,
1281  sigma_np1, dum2, compcond, dum4, dum5, dum6, dum7);
1282  }
1283 
1284  const std::string dispname=sup_previous_and_dot_to_varname(varnames[0]);
1285  const std::string xi =sup_previous_and_dot_to_varname(varnames[1]);
1286 
1287  if (unknowns_type == DISPLACEMENT_AND_PLASTIC_MULTIPLIER) {
1288  std::string expr = ("("+sigma_np1+"):Grad_Test_"+dispname
1289  + "+("+compcond+")*Test_"+xi);
1290  return add_nonlinear_term
1291  (md, mim, expr, region, false, false,
1292  "Small strain isotropic perfect elastoplasticity brick");
1293  } else {
1294  return add_nonlinear_term
1295  (md, mim, "("+sigma_np1+"):Grad_Test_"+dispname, region, true, false,
1296  "Small strain isotropic perfect elastoplasticity brick");
1297  }
1298  }
1299 
1300  void small_strain_elastoplasticity_next_iter
1301  (model &md, const mesh_im &mim,
1302  std::string lawname, plasticity_unknowns_type unknowns_type,
1303  const std::vector<std::string> &varnames,
1304  const std::vector<std::string> &params, size_type region) {
1305 
1306  filter_lawname(lawname);
1307  std::string Epnp1, xi_np1, alphanp1;
1308  {
1309  std::string dum1, dum3, dum5, dum6;
1310  build_isotropic_perfect_elastoplasticity_expressions_generic
1311  (md, lawname, unknowns_type, varnames, params,
1312  dum1, Epnp1, dum3, xi_np1, dum5, dum6, alphanp1);
1313  }
1314 
1315  std::string disp = sup_previous_and_dot_to_varname(varnames[0]);
1316  std::string xi = sup_previous_and_dot_to_varname(varnames[1]);
1317  std::string Previous_Ep = varnames[2];
1318 
1319  std::string Previous_alpha;
1320  base_vector tmpv_alpha;
1321  if (alphanp1.size()) { // Interpolation of the accumulated plastic strain
1322  Previous_alpha = varnames[3];
1323  tmpv_alpha.resize(gmm::vect_size(md.real_variable(Previous_alpha)));
1324  const im_data *pimd = md.pim_data_of_variable(Previous_alpha);
1325  if (pimd)
1326  ga_interpolation_im_data(md, alphanp1, *pimd, tmpv_alpha, region);
1327  else {
1328  const mesh_fem *pmf = md.pmesh_fem_of_variable(Previous_alpha);
1329  GMM_ASSERT1(pmf, "Provided data " << Previous_alpha
1330  << " should be defined on a im_data or a mesh_fem object");
1331  ga_local_projection(md, mim, alphanp1, *pmf, tmpv_alpha, region);
1332  }
1333  }
1334 
1335  base_vector tmpv_xi;
1336  if (xi_np1.size()) { // Interpolation of the plastic multiplier for the
1337  // theta-scheme and the case without multiplier (return mapping)
1338  // Not really necessary for the Backward Euler scheme.
1339  tmpv_xi.resize(gmm::vect_size(md.real_variable(xi)));
1340  const im_data *pimd = md.pim_data_of_variable(xi);
1341  if (pimd)
1342  ga_interpolation_im_data(md, xi_np1, *pimd, tmpv_xi, region);
1343  else {
1344  const mesh_fem *pmf = md.pmesh_fem_of_variable(xi);
1345  GMM_ASSERT1(pmf, "Provided data " << xi
1346  << " should be defined on a im_data or a mesh_fem object");
1347  ga_local_projection(md, mim, xi_np1, *pmf, tmpv_xi, region);
1348  }
1349  }
1350 
1351  base_vector tmpv_ep(gmm::vect_size(md.real_variable(Previous_Ep)));
1352  const im_data *pimd = md.pim_data_of_variable(Previous_Ep);
1353  if (pimd)
1354  ga_interpolation_im_data(md, Epnp1, *pimd, tmpv_ep, region);
1355  else {
1356  const mesh_fem *pmf = md.pmesh_fem_of_variable(Previous_Ep);
1357  GMM_ASSERT1(pmf, "Provided data " << Previous_Ep
1358  << " should be defined on a im_data or a mesh_fem object");
1359  ga_local_projection(md, mim, Epnp1, *pmf, tmpv_ep, region);
1360  }
1361 
1362  if (xi_np1.size())
1363  gmm::copy(tmpv_xi, md.set_real_variable(xi));
1364  if (alphanp1.size())
1365  gmm::copy(tmpv_alpha, md.set_real_variable(Previous_alpha));
1366  gmm::copy(tmpv_ep, md.set_real_variable(Previous_Ep));
1367  gmm::copy(md.real_variable(disp), md.set_real_variable("Previous_"+disp));
1368  gmm::copy(md.real_variable(xi), md.set_real_variable("Previous_"+xi));
1369  }
1370 
1371  // To be called after next_iter, not before
1372  void compute_small_strain_elastoplasticity_Von_Mises
1373  (model &md, const mesh_im &mim,
1374  std::string lawname, plasticity_unknowns_type unknowns_type,
1375  const std::vector<std::string> &varnames,
1376  const std::vector<std::string> &params,
1377  const mesh_fem &mf_vm, model_real_plain_vector &VM, size_type region) {
1378 
1379  GMM_ASSERT1(mf_vm.get_qdim() == 1,
1380  "Von mises stress can only be approximated on a scalar fem");
1381  VM.resize(mf_vm.nb_dof());
1382 
1383  filter_lawname(lawname);
1384 
1385  std::string sigma_after, von_mises;
1386  {
1387  std::string dum1, dum2, dum3, dum4, dum7;
1388  build_isotropic_perfect_elastoplasticity_expressions_generic
1389  (md, lawname, unknowns_type, varnames, params,
1390  dum1, dum2, dum3, dum4, sigma_after, von_mises, dum7);
1391  }
1392 
1393  size_type n_ep = 2; // Index of the plastic strain variable
1394 
1395  const im_data *pimd = md.pim_data_of_variable(varnames[n_ep]);
1396  if (pimd) {
1397  ga_local_projection(md, mim, von_mises, mf_vm, VM, region);
1398  }
1399  else {
1400  const mesh_fem *pmf = md.pmesh_fem_of_variable(varnames[n_ep]);
1401  GMM_ASSERT1(pmf, "Provided data " << varnames[n_ep]
1402  << " should be defined on a im_data or a mesh_fem object");
1403  ga_interpolation_Lagrange_fem(md, von_mises, mf_vm, VM, region);
1404  }
1405  }
1406 
1407 
1408 
1409  // ==============================
1410  //
1411  // Finite strain elastoplasticity
1412  //
1413  // ==============================
1414 
1415  const std::string _TWOTHIRD_("0.6666666666666666667");
1416  const std::string _FIVETHIRD_("1.6666666666666666667");
1417  const std::string _SQRTTHREEHALF_("1.2247448713915889407");
1418 
1419  void ga_define_linear_hardening_function
1420  (const std::string &name, scalar_type sigma_y0, scalar_type H, bool frobenius)
1421  {
1422  if (frobenius) {
1423  sigma_y0 *= sqrt(2./3.);
1424  H *= 2./3.;
1425  }
1426  std::stringstream expr, der;
1427  expr << std::setprecision(17) << sigma_y0 << "+" << H << "*t";
1428  der << std::setprecision(17) << H;
1429  ga_define_function(name, 1, expr.str(), der.str());
1430  }
1431 
1432  void ga_define_Ramberg_Osgood_hardening_function
1433  (const std::string &name,
1434  scalar_type sigma_ref, scalar_type eps_ref, scalar_type n, bool frobenius)
1435  {
1436  scalar_type coef = sigma_ref / pow(eps_ref, 1./n);
1437  if (frobenius)
1438  coef *= pow(2./3., 0.5 + 0.5/n); // = sqrt(2/3) * sqrt(2/3)^(1/n)
1439 
1440  std::stringstream expr, der;
1441  expr << std::setprecision(17) << coef << "*pow(t+1e-12," << 1./n << ")";
1442  der << std::setprecision(17) << coef/n << "*pow(t+1e-12," << 1./n-1 << ")";
1443  ga_define_function(name, 1, expr.str(), der.str());
1444  }
1445 
1446  // Simo-Miehe
1447  void build_Simo_Miehe_elastoplasticity_expressions
1448  (model &md, plasticity_unknowns_type unknowns_type,
1449  const std::vector<std::string> &varnames,
1450  const std::vector<std::string> &params,
1451  std::string &expr, std::string &plaststrain, std::string &invCp, std::string &vm)
1452  {
1453  GMM_ASSERT1(unknowns_type == DISPLACEMENT_AND_PLASTIC_MULTIPLIER ||
1454  unknowns_type == DISPLACEMENT_AND_PLASTIC_MULTIPLIER_AND_PRESSURE,
1455  "Not supported type of unknowns for this type of plasticity law");
1456  bool has_pressure_var(unknowns_type ==
1457  DISPLACEMENT_AND_PLASTIC_MULTIPLIER_AND_PRESSURE);
1458  GMM_ASSERT1(varnames.size() == (has_pressure_var ? 5 : 4),
1459  "Wrong number of variables.");
1460  GMM_ASSERT1(params.size() == 3, "Wrong number of parameters, "
1461  << params.size() << " != 3.");
1462  const std::string &dispname = sup_previous_and_dot_to_varname(varnames[0]);
1463  const std::string &multname = sup_previous_and_dot_to_varname(varnames[1]);
1464  const std::string &pressname = has_pressure_var ? varnames[2] : "";
1465  const std::string &plaststrain0 = varnames[has_pressure_var ? 3 : 2];
1466  const std::string &invCp0 = varnames[has_pressure_var ? 4 : 3];
1467  const std::string &K = params[0];
1468  const std::string &G = params[1];
1469  const std::string &sigma_y = params[2];
1470 
1471  const mesh_fem *mfu = md.pmesh_fem_of_variable(dispname);
1472  GMM_ASSERT1(mfu, "The provided displacement variable " << dispname <<
1473  " has to be defined on a mesh_fem");
1474  size_type N = mfu->linked_mesh().dim();
1475  GMM_ASSERT1(N >= 2 && N <= 3,
1476  "Finite strain elastoplasticity brick works only in 2D or 3D");
1477  GMM_ASSERT1(mfu && mfu->get_qdim() == N, "The finite strain "
1478  "elastoplasticity brick can only be applied on a fem "
1479  "variable of the same dimension as the mesh");
1480  const mesh_fem *mfmult = md.pmesh_fem_of_variable(multname);
1481  GMM_ASSERT1(mfmult && mfmult->get_qdim() == 1, "The plastic multiplier "
1482  "for the finite strain elastoplasticity brick has to be a "
1483  "scalar fem variable");
1484  bool mixed(!pressname.empty());
1485  const mesh_fem *mfpress = mixed ? md.pmesh_fem_of_variable(pressname) : 0;
1486  GMM_ASSERT1(!mixed || (mfpress && mfpress->get_qdim() == 1),
1487  "The hydrostatic pressure multiplier for the finite strain "
1488  "elastoplasticity brick has to be a scalar fem variable");
1489 
1490  GMM_ASSERT1(ga_function_exists(sigma_y), "The provided isotropic "
1491  "hardening function name '" << sigma_y << "' is not defined");
1492 
1493  GMM_ASSERT1(md.is_data(plaststrain0) &&
1494  (md.pim_data_of_variable(plaststrain0) ||
1495  md.pmesh_fem_of_variable(plaststrain0)),
1496  "The provided name '" << plaststrain0 << "' for the plastic "
1497  "strain field at the previous timestep, should be defined "
1498  "either as fem or as im data");
1499  GMM_ASSERT1((md.pim_data_of_variable(plaststrain0) &&
1500  md.pim_data_of_variable(plaststrain0)->nb_tensor_elem() == 1) ||
1501  (md.pmesh_fem_of_variable(plaststrain0) &&
1502  md.pmesh_fem_of_variable(plaststrain0)->get_qdim() == 1),
1503  "Wrong size of " << plaststrain0);
1504  GMM_ASSERT1(md.is_data(invCp0) &&
1505  (md.pim_data_of_variable(invCp0) ||
1506  md.pmesh_fem_of_variable(invCp0)),
1507  "The provided name '" << invCp0 << "' for the inverse of the "
1508  "plastic right Cauchy-Green tensor field at the previous "
1509  "timestep, should be defined either as fem or as im data");
1510  bgeot::multi_index Cp_size(1);
1511  Cp_size[0] = 4 + (N==3)*2;
1512  GMM_ASSERT1((md.pim_data_of_variable(invCp0) &&
1513  md.pim_data_of_variable(invCp0)->tensor_size() == Cp_size) ||
1514  (md.pmesh_fem_of_variable(invCp0) &&
1515  md.pmesh_fem_of_variable(invCp0)->get_qdims() == Cp_size),
1516  "Wrong size of " << invCp0);
1517 
1518  const std::string _U_ = sup_previous_and_dot_to_varname(dispname);
1519  const std::string _KSI_ = sup_previous_and_dot_to_varname(multname);
1520  const std::string _I_(N == 2 ? "Id(2)" : "Id(3)");
1521  const std::string _F_("("+_I_+"+Grad_"+_U_+")");
1522  const std::string _J_("Det"+_F_); // in 2D assumes plane strain
1523 
1524  std::string _P_;
1525  if (mixed)
1526  _P_ = "-"+pressname+"*"+_J_;
1527  else
1528  _P_ = "("+K+")*log("+_J_+")";
1529 
1530  std::string _INVCP0_, _F3d_, _DEVLOGBETR_, _DEVLOGBETR_3D_;
1531  if (N == 2) { // plane strain
1532  _INVCP0_ = "([[[1,0,0],[0,0,0],[0,0,0]],"
1533  "[[0,0,0],[0,1,0],[0,0,0]],"
1534  "[[0,0,0],[0,0,0],[0,0,1]],"
1535  "[[0,1,0],[1,0,0],[0,0,0]]]."+invCp0+")";
1536  _F3d_ = "(Id(3)+[[1,0,0],[0,1,0]]*Grad_"+_U_+"*[[1,0,0],[0,1,0]]')";
1537  _DEVLOGBETR_3D_ = "(Deviator(Logm("+_F3d_+"*"+_INVCP0_+"*"+_F3d_+"')))";
1538  _DEVLOGBETR_ = "([[[[1,0],[0,0]],[[0,1],[0,0]],[[0,0],[0,0]]],"
1539  "[[[0,0],[1,0]],[[0,0],[0,1]],[[0,0],[0,0]]],"
1540  "[[[0,0],[0,0]],[[0,0],[0,0]],[[0,0],[0,0]]]]"
1541  ":"+_DEVLOGBETR_3D_+")";
1542  } else { // 3D
1543  _INVCP0_ = "([[[1,0,0],[0,0,0],[0,0,0]],"
1544  "[[0,0,0],[0,1,0],[0,0,0]],"
1545  "[[0,0,0],[0,0,0],[0,0,1]],"
1546  "[[0,1,0],[1,0,0],[0,0,0]],"
1547  "[[0,0,1],[0,0,0],[1,0,0]],"
1548  "[[0,0,0],[0,0,1],[0,1,0]]]."+invCp0+")";
1549  _F3d_ = _F_;
1550  _DEVLOGBETR_3D_ =
1551  _DEVLOGBETR_ = "(Deviator(Logm("+_F_+"*"+_INVCP0_+"*"+_F_+"')))";
1552  }
1553  const std::string _DEVTAUTR_("("+G+"*"+_DEVLOGBETR_+")");
1554  const std::string _DEVTAUTR_3D_("("+G+"*"+_DEVLOGBETR_3D_+")");
1555  const std::string _DEVTAU_("((1-2*"+_KSI_+")*"+_DEVTAUTR_+")");
1556  const std::string _DEVTAU_3D_("((1-2*"+_KSI_+")*"+_DEVTAUTR_3D_+")");
1557  const std::string _TAU_("("+_P_+"*"+_I_+"+"+_DEVTAU_+")");
1558 
1559  const std::string _PLASTSTRAIN_("("+plaststrain0+"+"+_KSI_+"*Norm"+_DEVLOGBETR_3D_+")");
1560  const std::string _SIGMA_Y_(sigma_y+"("+_PLASTSTRAIN_+")");
1561 
1562  // results
1563  expr = _TAU_+":(Grad_Test_"+_U_+"*Inv"+_F_+")";
1564  if (mixed)
1565  expr += "+("+pressname+"/("+K+")+log("+_J_+")/"+_J_+")*Test_"+pressname;
1566  expr += "+(Norm"+_DEVTAU_+
1567  "-min("+_SIGMA_Y_+",Norm"+_DEVTAUTR_+"+1e-12*"+_KSI_+"))*Test_"+_KSI_;
1568 
1569  plaststrain = _PLASTSTRAIN_;
1570 
1571  if (N==2) invCp = "[[[1,0,0,0.0],[0,0,0,0.5],[0,0,0,0]],"
1572  "[[0,0,0,0.5],[0,1,0,0.0],[0,0,0,0]],"
1573  "[[0,0,0,0.0],[0,0,0,0.0],[0,0,1,0]]]";
1574  else invCp = "[[[1.0,0,0,0,0,0],[0,0,0,0.5,0,0],[0,0,0,0,0.5,0]],"
1575  "[[0,0,0,0.5,0,0],[0,1.0,0,0,0,0],[0,0,0,0,0,0.5]],"
1576  "[[0,0,0,0,0.5,0],[0,0,0,0,0,0.5],[0,0,1.0,0,0,0]]]";
1577  invCp += ":((Inv"+_F3d_+"*Expm(-"+_KSI_+"*"+_DEVLOGBETR_3D_+")*"+_F3d_+")*"+_INVCP0_+
1578  "*(Inv"+_F3d_+"*Expm(-"+_KSI_+"*"+_DEVLOGBETR_3D_+")*"+_F3d_+")')";
1579 
1580  vm = _SQRTTHREEHALF_+"*Norm("+_DEVTAU_+")/"+_J_;
1581  }
1582 
1583  size_type add_finite_strain_elastoplasticity_brick
1584  (model &md, const mesh_im &mim,
1585  std::string lawname, plasticity_unknowns_type unknowns_type,
1586  const std::vector<std::string> &varnames,
1587  const std::vector<std::string> &params, size_type region)
1588  {
1589  filter_lawname(lawname);
1590  if (lawname.compare("simo_miehe") == 0 ||
1591  lawname.compare("eterovic_bathe") == 0) {
1592  std::string expr, dummy1, dummy2, dummy3;
1593  build_Simo_Miehe_elastoplasticity_expressions
1594  (md, unknowns_type, varnames, params, expr, dummy1, dummy2, dummy3);
1595  return add_nonlinear_term
1596  (md, mim, expr, region, true, false, "Simo Miehe elastoplasticity brick");
1597  } else
1598  GMM_ASSERT1(false, lawname << " is not a known elastoplastic law");
1599  }
1600 
1601  // Updates any state variables included in params for the given lawname
1602  void finite_strain_elastoplasticity_next_iter
1603  (model &md, const mesh_im &mim,
1604  std::string lawname, plasticity_unknowns_type unknowns_type,
1605  const std::vector<std::string> &varnames,
1606  const std::vector<std::string> &params, size_type region) {
1607 
1608  filter_lawname(lawname);
1609  if (lawname.compare("simo_miehe") == 0 ||
1610  lawname.compare("eterovic_bathe") == 0) {
1611  std::string dummy1, dummy2, plaststrain, invCp;
1612  build_Simo_Miehe_elastoplasticity_expressions
1613  (md, unknowns_type, varnames, params, dummy1, plaststrain, invCp, dummy2);
1614 
1615  bool has_pressure_var(unknowns_type ==
1616  DISPLACEMENT_AND_PLASTIC_MULTIPLIER_AND_PRESSURE);
1617  const std::string &multname = sup_previous_and_dot_to_varname(varnames[1]);
1618  const std::string &plaststrain0 = varnames[has_pressure_var ? 3 : 2];
1619  const std::string &invCp0 = varnames[has_pressure_var ? 4 : 3];
1620  { // update plaststrain0
1621  model_real_plain_vector tmpvec(gmm::vect_size(md.real_variable(plaststrain0)));
1622  const im_data *pimd = md.pim_data_of_variable(plaststrain0);
1623  if (pimd)
1624  ga_interpolation_im_data(md, plaststrain, *pimd, tmpvec, region);
1625  else {
1626  const mesh_fem *pmf = md.pmesh_fem_of_variable(plaststrain0);
1627  GMM_ASSERT1(pmf, "Provided data " << plaststrain0 << " should be defined "
1628  "either on a im_data or a mesh_fem object");
1629  //ga_interpolation_Lagrange_fem(md, plaststrain, *pmf, tmpvec, region);
1630  ga_local_projection(md, mim, plaststrain, *pmf, tmpvec, region);
1631  }
1632  gmm::copy(tmpvec, md.set_real_variable(plaststrain0));
1633  }
1634 
1635  { // update invCp0
1636  model_real_plain_vector tmpvec(gmm::vect_size(md.real_variable(invCp0)));
1637  const im_data *pimd = md.pim_data_of_variable(invCp0);
1638  if (pimd)
1639  ga_interpolation_im_data(md, invCp, *pimd, tmpvec, region);
1640  else {
1641  const mesh_fem *pmf = md.pmesh_fem_of_variable(invCp0);
1642  GMM_ASSERT1(pmf, "Provided data " << invCp0 << " should be defined "
1643  "either on a im_data or a mesh_fem object");
1644  //ga_interpolation_Lagrange_fem(md, invCp, *pmf, tmpvec, region);
1645  ga_local_projection(md, mim, invCp, *pmf, tmpvec, region);
1646  }
1647  gmm::copy(tmpvec, md.set_real_variable(invCp0));
1648  }
1649 
1650  gmm::clear(md.set_real_variable(multname));
1651 
1652  } else
1653  GMM_ASSERT1(false, lawname << " is not a known elastoplastic law");
1654  }
1655 
1656  void compute_finite_strain_elastoplasticity_Von_Mises
1657  (model &md, const mesh_im &mim,
1658  std::string lawname, plasticity_unknowns_type unknowns_type,
1659  const std::vector<std::string> &varnames,
1660  const std::vector<std::string> &params,
1661  const mesh_fem &mf_vm, model_real_plain_vector &VM,
1662  size_type region) {
1663 
1664  filter_lawname(lawname);
1665  if (lawname.compare("simo_miehe") == 0 ||
1666  lawname.compare("eterovic_bathe") == 0) {
1667  std::string dummy1, dummy2, dummy3, von_mises;
1668  build_Simo_Miehe_elastoplasticity_expressions
1669  (md, unknowns_type, varnames, params, dummy1, dummy2, dummy3, von_mises);
1670  VM.resize(mf_vm.nb_dof());
1671 
1672  bool has_pressure_var(unknowns_type ==
1673  DISPLACEMENT_AND_PLASTIC_MULTIPLIER_AND_PRESSURE);
1674  const std::string &plaststrain0 = varnames[has_pressure_var ? 3 : 2];
1675  const std::string &invCp0 = varnames[has_pressure_var ? 4 : 3];
1676  bool im_data1 = md.pim_data_of_variable(plaststrain0) != 0;
1677  bool im_data2 = md.pim_data_of_variable(invCp0) != 0;
1678  bool fem_data1 = md.pmesh_fem_of_variable(plaststrain0) != 0;
1679  bool fem_data2 = md.pmesh_fem_of_variable(invCp0) != 0;
1680  GMM_ASSERT1(im_data1 || fem_data1,
1681  "Provided data " << plaststrain0 <<
1682  " should be defined on a im_data or a mesh_fem object");
1683  GMM_ASSERT1(im_data2 || fem_data2,
1684  "Provided data " << invCp0 <<
1685  " should be defined on a im_data or a mesh_fem object");
1686  if (im_data1 || im_data2) {
1687  ga_local_projection(md, mim, von_mises, mf_vm, VM, region);
1688  } else {
1689  ga_interpolation_Lagrange_fem(md, von_mises, mf_vm, VM, region);
1690  }
1691 
1692  } else
1693  GMM_ASSERT1(false, lawname << " is not a known elastoplastic law");
1694 
1695  }
1696 
1697 
1698 
1699 
1700 
1701 
1702 
1703 
1704 
1705 
1706 
1707 
1708 
1709 
1710 
1711 
1712 
1713  //=================================================================
1714  //
1715  // Old version of an elastoplasticity Brick for isotropic perfect
1716  // plasticity with the low level generic assembly.
1717  // Particularity of this brick: the flow rule is integrated on
1718  // finite element nodes (not on Gauss points).
1719  //
1720  //=================================================================
1721 
1722  enum elastoplasticity_nonlinear_term_version { PROJ,
1723  GRADPROJ,
1724  PLAST
1725  };
1726 
1727  /** Compute the projection of D*e + sigma_bar_
1728  on the dof of sigma. */
1729  class elastoplasticity_nonlinear_term : public nonlinear_elem_term {
1730 
1731  protected:
1732  const mesh_im &mim;
1733  const mesh_fem &mf_u;
1734  const mesh_fem &mf_sigma;
1735  const mesh_fem *pmf_data;
1736  model_real_plain_vector U_n, U_np1;
1737  model_real_plain_vector Sigma_n;
1738  model_real_plain_vector threshold, lambda, mu;
1739  const abstract_constraints_projection &t_proj;
1740  const size_type option;
1741  const size_type flag_proj;
1742  const bool store_sigma;
1743 
1744  bgeot::multi_index sizes_;
1745 
1746  size_type N, size_proj;
1747 
1748  // temporary variables
1749  base_vector params;
1750  size_type current_cv;
1751  model_real_plain_vector convex_coeffs, interpolated_val;
1752 
1753  // storage variables
1754  model_real_plain_vector cumulated_sigma; // either the projected stress (option==PROJ)
1755  // or the plastic stress (option==PLAST)
1756  model_real_plain_vector cumulated_count;
1757 
1758  fem_precomp_pool fppool;
1759 
1760 
1761  // computes stresses or stress projections on all sigma dofs of a convex
1762  void compute_convex_coeffs(size_type cv) {
1763 
1764  current_cv = cv;
1765 
1766  pfem pf_sigma = mf_sigma.fem_of_element(cv);
1767  size_type nbd_sigma = pf_sigma->nb_dof(cv);
1768  size_type qdim_sigma = mf_sigma.get_qdim();
1769 
1770  gmm::resize(convex_coeffs, size_proj*nbd_sigma);
1771 
1772  base_matrix G;
1773  bgeot::vectors_to_base_matrix
1774  (G, mf_u.linked_mesh().points_of_convex(cv));
1775  bgeot::pgeometric_trans pgt =
1776  mf_u.linked_mesh().trans_of_convex(cv);
1777 
1778  // if the Lame coefficient are vector fields
1779  base_vector coeff_data;
1780  pfem pf_data;
1781  fem_interpolation_context ctx_data;
1782  if (pmf_data) {
1783  pf_data = pmf_data->fem_of_element(cv);
1784  size_type nbd_data = pf_data->nb_dof(cv);
1785  coeff_data.resize(nbd_data*3);
1786 
1787  // Definition of the Lame coeff
1788  mesh_fem::ind_dof_ct::const_iterator itdof
1789  = pmf_data->ind_basic_dof_of_element(cv).begin();
1790  for (size_type k = 0; k < nbd_data; ++k, ++itdof) {
1791  coeff_data[k*3] = lambda[*itdof];
1792  coeff_data[k*3+1] = mu[*itdof];
1793  coeff_data[k*3+2] = threshold[*itdof];
1794  }
1795  GMM_ASSERT1(pf_data->target_dim() == 1,
1796  "won't interpolate on a vector FEM... ");
1797 
1798  pfem_precomp pfp_data = fppool(pf_data, pf_sigma->node_tab(cv));
1799  ctx_data = fem_interpolation_context
1800  (pgt, pfp_data, size_type(-1), G, cv, short_type(-1));
1801  }
1802 
1803  // Definition of the coeff for du = u_n-u_np1 and optionally for u_np1
1804  size_type cvnbdof_u = mf_u.nb_basic_dof_of_element(cv);
1805  model_real_plain_vector coeff_du(cvnbdof_u);
1806  model_real_plain_vector coeff_u_np1(cvnbdof_u);
1807  mesh_fem::ind_dof_ct::const_iterator itdof
1808  = mf_u.ind_basic_dof_of_element(cv).begin();
1809  for (size_type k = 0; k < cvnbdof_u; ++k, ++itdof) {
1810  coeff_du[k] = U_np1[*itdof] - U_n[*itdof];
1811  coeff_u_np1[k] = U_np1[*itdof];
1812  }
1813 
1814  pfem pf_u = mf_u.fem_of_element(cv);
1815  pfem_precomp pfp_u = fppool(pf_u, pf_sigma->node_tab(cv));
1816  fem_interpolation_context
1817  ctx_u(pgt, pfp_u, size_type(-1), G, cv, short_type(-1));
1818 
1819  size_type qdim = mf_u.get_qdim();
1820  base_matrix G_du(qdim, qdim), G_u_np1(qdim, qdim); // G_du = G_u_np1 - G_u_n
1821 
1822  for (size_type ii = 0; ii < nbd_sigma; ++ii) {
1823 
1824  if (pmf_data) {
1825  // interpolation of the data on sigma dof
1826  ctx_data.set_ii(ii);
1827  pf_data->interpolation(ctx_data, coeff_data, params, 3);
1828  }
1829 
1830  // interpolation of the gradient of du and u_np1 on sigma dof
1831  ctx_u.set_ii(ii);
1832  pf_u->interpolation_grad(ctx_u, coeff_du, G_du, dim_type(qdim));
1833  if (option == PLAST)
1834  pf_u->interpolation_grad(ctx_u, coeff_u_np1, G_u_np1, dim_type(qdim));
1835 
1836  // Compute lambda*(tr(eps_np1)-tr(eps_n)) and lambda*tr(eps_np1)
1837  scalar_type ltrace_deps = params[0]*gmm::mat_trace(G_du);
1838  scalar_type ltrace_eps_np1 = (option == PLAST) ?
1839  params[0]*gmm::mat_trace(G_u_np1) : 0.;
1840 
1841  // Compute sigma_hat = D*(eps_np1 - eps_n) + sigma_n
1842  // where D represents the elastic stiffness tensor
1843  base_matrix sigma_hat(qdim, qdim);
1844  size_type sigma_dof = mf_sigma.ind_basic_dof_of_element(cv)[ii*qdim_sigma];
1845  for (dim_type j = 0; j < qdim; ++j) {
1846  for (dim_type i = 0; i < qdim; ++i)
1847  sigma_hat(i,j) = Sigma_n[sigma_dof++]
1848  + params[1]*(G_du(i,j) + G_du(j,i));
1849  sigma_hat(j,j) += ltrace_deps;
1850  }
1851 
1852  // Compute the projection or its grad
1853  base_matrix proj;
1854  t_proj.do_projection(sigma_hat, params[2], proj, flag_proj);
1855 
1856  // Compute the plastic part if required
1857  if (option == PLAST)
1858  for (dim_type i = 0; i < qdim; ++i) {
1859  for (dim_type j = 0; j < qdim; ++j)
1860  proj(i,j) -= params[1]*(G_u_np1(i,j) + G_u_np1(j,i));
1861  proj(i,i) -= ltrace_eps_np1;
1862  }
1863 
1864  // Fill in convex_coeffs with sigma or its grad
1865  std::copy(proj.begin(), proj.end(),
1866  convex_coeffs.begin() + proj.size() * ii);
1867 
1868  // Store the projected or plastic sigma
1869  if (store_sigma) {
1870  sigma_dof = mf_sigma.ind_basic_dof_of_element(cv)[ii*qdim_sigma];
1871  for (dim_type j = 0; j < qdim; ++j) {
1872  for (dim_type i = 0; i < qdim; ++i) {
1873  cumulated_count[sigma_dof] += 1;
1874  cumulated_sigma[sigma_dof++] += proj(i,j);
1875  }
1876  }
1877  }
1878 
1879  } // ii = 0:nbd_sigma-1
1880 
1881  }
1882 
1883  public:
1884 
1885  // constructor
1886  elastoplasticity_nonlinear_term
1887  (const mesh_im &mim_,
1888  const mesh_fem &mf_u_,
1889  const mesh_fem &mf_sigma_,
1890  const mesh_fem *pmf_data_,
1891  const model_real_plain_vector &U_n_,
1892  const model_real_plain_vector &U_np1_,
1893  const model_real_plain_vector &Sigma_n_,
1894  const model_real_plain_vector &threshold_,
1895  const model_real_plain_vector &lambda_,
1896  const model_real_plain_vector &mu_,
1897  const abstract_constraints_projection &t_proj_,
1898  size_type option_,
1899  bool store_sigma_) :
1900  mim(mim_), mf_u(mf_u_), mf_sigma(mf_sigma_), pmf_data(pmf_data_),
1901  Sigma_n(Sigma_n_), t_proj(t_proj_), option(option_),
1902  flag_proj(option == GRADPROJ ? 1 : 0),
1903  store_sigma(option == GRADPROJ ? false : store_sigma_) {
1904 
1905  params.resize(3);
1906  N = mf_u.linked_mesh().dim();
1907 
1908  sizes_ = (flag_proj == 0 ? bgeot::multi_index(N,N)
1909  : bgeot::multi_index(N,N,N,N));
1910 
1911  // size_proj is different if we compute the projection
1912  // or the gradient of the projection
1913  size_proj = (flag_proj == 0 ? N*N : N*N*N*N);
1914 
1915  gmm::resize(U_n, mf_u.nb_basic_dof());
1916  gmm::resize(U_np1, mf_u.nb_basic_dof());
1917  gmm::resize(Sigma_n, mf_sigma.nb_basic_dof());
1918  mf_u.extend_vector(gmm::sub_vector(U_n_,
1919  gmm::sub_interval(0,mf_u.nb_dof())),
1920  U_n);
1921  mf_u.extend_vector(gmm::sub_vector(U_np1_,
1922  gmm::sub_interval(0,mf_u.nb_dof())),
1923  U_np1);
1924  mf_sigma.extend_vector(gmm::sub_vector(Sigma_n_,
1925  gmm::sub_interval(0,mf_sigma.nb_dof())),
1926  Sigma_n);
1927 
1928  if (pmf_data) {
1929  gmm::resize(mu, pmf_data->nb_basic_dof());
1930  gmm::resize(lambda, pmf_data->nb_basic_dof());
1931  gmm::resize(threshold, pmf_data->nb_basic_dof());
1932  pmf_data->extend_vector(threshold_, threshold);
1933  pmf_data->extend_vector(lambda_, lambda);
1934  pmf_data->extend_vector(mu_, mu);
1935  } else {
1936  gmm::resize(mu, 1); mu[0] = mu_[0];
1937  gmm::resize(lambda, 1); lambda[0] = lambda_[0];
1938  gmm::resize(threshold, 1); threshold[0] = threshold_[0];
1939  params[0] = lambda[0];
1940  params[1] = mu[0];
1941  params[2] = threshold[0];
1942  }
1943  GMM_ASSERT1(mf_u.get_qdim() == N,
1944  "wrong qdim for the mesh_fem");
1945 
1946  gmm::resize(interpolated_val, size_proj);
1947 
1948  if (store_sigma) {
1949  cumulated_sigma.resize(mf_sigma.nb_dof());
1950  cumulated_count.resize(mf_sigma.nb_dof());
1951  }
1952 
1953  // used to know if the current element is different
1954  // than the previous one and so if a new computation
1955  // is necessary or not.
1956  current_cv = size_type(-1);
1957 
1958  }
1959 
1960 
1961  const bgeot::multi_index &sizes(size_type) const { return sizes_; }
1962 
1963 
1964  // method from nonlinear_elem_term, gives on output the tensor
1965  virtual void compute(fem_interpolation_context& ctx,
1966  bgeot::base_tensor &t) {
1967  size_type cv = ctx.convex_num(); //index of current element
1968  pfem pf_sigma = ctx.pf();
1969  GMM_ASSERT1(pf_sigma->is_lagrange(),
1970  "Sorry, works only for Lagrange fems");
1971 
1972  // if the current element is different than the previous one
1973  if (cv != current_cv)
1974  compute_convex_coeffs(cv);
1975 
1976  // interpolation of the sigma or its grad on sigma dof
1977  pf_sigma->interpolation(ctx, convex_coeffs, interpolated_val, dim_type(size_proj));
1978 
1979  // copy the result into the returned tensor t
1980  t.adjust_sizes(sizes_);
1981  std::copy(interpolated_val.begin(), interpolated_val.end(), t.begin());
1982  }
1983 
1984 
1985  // method to get the averaged sigma stored during the assembly
1986  void get_averaged_sigmas(model_real_plain_vector &sigma) {
1987  model_real_plain_vector glob_cumulated_count(mf_sigma.nb_dof());
1988  MPI_SUM_VECTOR(cumulated_sigma, sigma);
1989  MPI_SUM_VECTOR(cumulated_count, glob_cumulated_count);
1990  size_type imax = mf_sigma.nb_dof();
1991  for (size_type i = 0; i < imax; ++i)
1992  sigma[i] /= glob_cumulated_count[i];
1993  }
1994 
1995 };
1996 
1997 
1998 
1999  /**
2000  Right hand side vector for elastoplasticity
2001  @ingroup asm
2002  */
2003  void asm_elastoplasticity_rhs
2004  (model_real_plain_vector &V,
2005  model_real_plain_vector *saved_sigma,
2006  const mesh_im &mim,
2007  const mesh_fem &mf_u,
2008  const mesh_fem &mf_sigma,
2009  const mesh_fem &mf_data,
2010  const model_real_plain_vector &u_n,
2011  const model_real_plain_vector &u_np1,
2012  const model_real_plain_vector &sigma_n,
2013  const model_real_plain_vector &lambda,
2014  const model_real_plain_vector &mu,
2015  const model_real_plain_vector &threshold,
2016  const abstract_constraints_projection &t_proj,
2017  size_type option_sigma,
2018  const mesh_region &rg = mesh_region::all_convexes()) {
2019 
2020  GMM_ASSERT1(mf_u.get_qdim() == mf_u.linked_mesh().dim(),
2021  "wrong qdim for the mesh_fem");
2022  GMM_ASSERT1(option_sigma == PROJ || option_sigma == PLAST,
2023  "wrong option parameter");
2024 
2025  elastoplasticity_nonlinear_term plast(mim, mf_u, mf_sigma, &mf_data,
2026  u_n, u_np1, sigma_n,
2027  threshold, lambda, mu,
2028  t_proj, option_sigma, (saved_sigma != NULL));
2029 
2030  generic_assembly assem("V(#1) + =comp(NonLin(#2).vGrad(#1))(i,j,:,i,j);");
2031 
2032  assem.push_mi(mim);
2033  assem.push_mf(mf_u);
2034  assem.push_mf(mf_sigma);
2035  assem.push_nonlinear_term(&plast);
2036  assem.push_vec(V);
2037  assem.assembly(rg);
2038 
2039  if (saved_sigma)
2040  plast.get_averaged_sigmas(*saved_sigma);
2041  }
2042 
2043 
2044  /**
2045  Tangent matrix for elastoplasticity
2046  @ingroup asm
2047  */
2048  void asm_elastoplasticity_tangent_matrix
2049  (model_real_sparse_matrix &H,
2050  const mesh_im &mim,
2051  const mesh_fem &mf_u,
2052  const mesh_fem &mf_sigma,
2053  const mesh_fem *pmf_data,
2054  const model_real_plain_vector &u_n,
2055  const model_real_plain_vector &u_np1,
2056  const model_real_plain_vector &sigma_n,
2057  const model_real_plain_vector &lambda,
2058  const model_real_plain_vector &mu,
2059  const model_real_plain_vector &threshold,
2060  const abstract_constraints_projection &t_proj,
2061  const mesh_region &rg = mesh_region::all_convexes()) {
2062 
2063  GMM_ASSERT1(mf_u.get_qdim() == mf_u.linked_mesh().dim(),
2064  "wrong qdim for the mesh_fem");
2065 
2066  elastoplasticity_nonlinear_term gradplast(mim, mf_u, mf_sigma, pmf_data,
2067  u_n, u_np1, sigma_n,
2068  threshold, lambda, mu,
2069  t_proj, GRADPROJ, false);
2070 
2071  generic_assembly assem;
2072 
2073  if (pmf_data)
2074  assem.set("lambda=data$1(#3); mu=data$2(#3);"
2075  "t=comp(NonLin(#2).vGrad(#1).vGrad(#1).Base(#3))(i,j,:,:,:,:,:,:,i,j,:);"
2076  "M(#1,#1)+= sym(t(k,l,:,l,k,:,m).mu(m)+t(k,l,:,k,l,:,m).mu(m)+t(k,k,:,l,l,:,m).lambda(m))");
2077  else
2078  assem.set("lambda=data$1(1); mu=data$2(1);"
2079  "t=comp(NonLin(#2).vGrad(#1).vGrad(#1))(i,j,:,:,:,:,:,:,i,j);"
2080  "M(#1,#1)+= sym(t(k,l,:,l,k,:).mu(1)+t(k,l,:,k,l,:).mu(1)+t(k,k,:,l,l,:).lambda(1))");
2081 
2082  assem.push_mi(mim);
2083  assem.push_mf(mf_u);
2084  assem.push_mf(mf_sigma);
2085  if (pmf_data)
2086  assem.push_mf(*pmf_data);
2087  assem.push_data(lambda);
2088  assem.push_data(mu);
2089  assem.push_nonlinear_term(&gradplast);
2090  assem.push_mat(H);
2091  assem.assembly(rg);
2092 
2093  }
2094 
2095 
2096 
2097  //=================================================================
2098  // Elastoplasticity Brick
2099  //=================================================================
2100 
2101  struct elastoplasticity_brick : public virtual_brick {
2102 
2103  pconstraints_projection t_proj;
2104 
2105  virtual void asm_real_tangent_terms(const model &md,
2106  size_type /* ib */,
2107  const model::varnamelist &vl,
2108  const model::varnamelist &dl,
2109  const model::mimlist &mims,
2110  model::real_matlist &matl,
2111  model::real_veclist &vecl,
2112  model::real_veclist &,
2113  size_type region,
2114  build_version version)const {
2115 
2116  GMM_ASSERT1(mims.size() == 1,
2117  "Elastoplasticity brick need a single mesh_im");
2118  GMM_ASSERT1(vl.size() == 1,
2119  "Elastoplasticity brick need one variable");
2120  /** vl[0] = u */
2121 
2122  GMM_ASSERT1(dl.size() == 5,
2123  "Wrong number of data for elastoplasticity brick, "
2124  << dl.size() << " should be 4.");
2125  GMM_ASSERT1(matl.size() == 1, "Wrong number of terms for "
2126  "elastoplasticity brick");
2127 
2128  const model_real_plain_vector &u_np1 = md.real_variable(vl[0]);
2129  const model_real_plain_vector &u_n = md.real_variable(dl[4]);
2130  const mesh_fem &mf_u = *(md.pmesh_fem_of_variable(vl[0]));
2131  GMM_ASSERT1(&mf_u == md.pmesh_fem_of_variable(dl[4]),
2132  "The previous displacement data have to be defined on the "
2133  "same mesh_fem as the displacement variable");
2134 
2135  const model_real_plain_vector &lambda = md.real_variable(dl[0]);
2136  const model_real_plain_vector &mu = md.real_variable(dl[1]);
2137  const model_real_plain_vector &threshold = md.real_variable(dl[2]);
2138  const mesh_fem *mf_data = md.pmesh_fem_of_variable(dl[0]);
2139 
2140  const model_real_plain_vector &sigma_n = md.real_variable(dl[3]);
2141  const mesh_fem &mf_sigma = *(md.pmesh_fem_of_variable(dl[3]));
2142  GMM_ASSERT1(!(mf_sigma.is_reduced()),
2143  "Works only for pure Lagrange fems");
2144 
2145  const mesh_im &mim = *mims[0];
2146  mesh_region rg(region);
2147  mim.linked_mesh().intersect_with_mpi_region(rg);
2148 
2149  if (version & model::BUILD_MATRIX) {
2150  gmm::clear(matl[0]);
2151  asm_elastoplasticity_tangent_matrix
2152  (matl[0], mim, mf_u, mf_sigma, mf_data, u_n,
2153  u_np1, sigma_n, lambda, mu, threshold, *t_proj, rg);
2154  }
2155 
2156  if (version & model::BUILD_RHS) {
2157  model_real_plain_vector *dummy = 0;
2158  asm_elastoplasticity_rhs(vecl[0], dummy,
2159  mim, mf_u, mf_sigma, *mf_data,
2160  u_n, u_np1, sigma_n,
2161  lambda, mu, threshold, *t_proj, PROJ, rg);
2162  gmm::scale(vecl[0], scalar_type(-1));
2163  }
2164 
2165  }
2166 
2167  // constructor
2168  elastoplasticity_brick(const pconstraints_projection &t_proj_)
2169  : t_proj(t_proj_) {
2170  set_flags("Elastoplasticity brick", false /* is linear*/,
2171  true /* is symmetric */, false /* is coercive */,
2172  true /* is real */, false /* is complex */);
2173  }
2174 
2175  };
2176 
2177 
2178  //=================================================================
2179  // Add a elastoplasticity brick
2180  //=================================================================
2181 
2182  size_type add_elastoplasticity_brick
2183  (model &md,
2184  const mesh_im &mim,
2185  const pconstraints_projection &ACP,
2186  const std::string &varname,
2187  const std::string &data_previous_disp,
2188  const std::string &datalambda,
2189  const std::string &datamu,
2190  const std::string &datathreshold,
2191  const std::string &datasigma,
2192  size_type region) {
2193 
2194  pbrick pbr = std::make_shared<elastoplasticity_brick>(ACP);
2195 
2196  model::termlist tl{model::term_description(varname, varname, true)};
2197  model::varnamelist
2198  dl{datalambda, datamu, datathreshold, datasigma, data_previous_disp},
2199  vl{varname};
2200 
2201  return md.add_brick(pbr, vl, dl, tl,
2202  model::mimlist(1,&mim), region);
2203  }
2204 
2205 
2206  //=================================================================
2207  // New stress constraints values computation and saved
2208  //=================================================================
2209  // Update of u and sigma on time iterates :
2210  // u_np1 -> u_n sigma_np1 -> sigma_n
2211 
2212  void elastoplasticity_next_iter(model &md,
2213  const mesh_im &mim,
2214  const std::string &varname,
2215  const std::string &data_previous_disp,
2216  const pconstraints_projection &ACP,
2217  const std::string &datalambda,
2218  const std::string &datamu,
2219  const std::string &datathreshold,
2220  const std::string &datasigma) {
2221 
2222  const model_real_plain_vector &u_np1 = md.real_variable(varname);
2223  model_real_plain_vector &u_n = md.set_real_variable(data_previous_disp);
2224  const mesh_fem &mf_u = *(md.pmesh_fem_of_variable(varname));
2225 
2226  const model_real_plain_vector &lambda = md.real_variable(datalambda);
2227  const model_real_plain_vector &mu = md.real_variable(datamu);
2228  const model_real_plain_vector &threshold = md.real_variable(datathreshold);
2229  const mesh_fem *mf_data = md.pmesh_fem_of_variable(datalambda);
2230 
2231  const model_real_plain_vector &sigma_n = md.real_variable(datasigma);
2232  const mesh_fem &mf_sigma = *(md.pmesh_fem_of_variable(datasigma));
2233 
2234  // dim_type N = mf_sigma.linked_mesh().dim();
2235 
2236  mesh_region rg = mim.linked_mesh().get_mpi_region();
2237 
2238  model_real_plain_vector sigma_np1(mf_sigma.nb_dof());
2239  model_real_plain_vector dummyV(mf_u.nb_dof());
2240  asm_elastoplasticity_rhs(dummyV, &sigma_np1,
2241  mim, mf_u, mf_sigma, *mf_data,
2242  u_n, u_np1, sigma_n,
2243  lambda, mu, threshold, *ACP, PROJ, rg);
2244 
2245  // upload sigma and u : u_np1 -> u_n, sigma_np1 -> sigma_n
2246  // be careful to use this function
2247  // only if the computation is over
2248  gmm::copy(sigma_np1, md.set_real_variable(datasigma));
2249  gmm::copy(u_np1, u_n);
2250 
2251  }
2252 
2253 
2254 
2255  //=================================================================
2256  // Von Mises or Tresca stress computation for elastoplasticity
2257  //=================================================================
2258 
2259  void compute_elastoplasticity_Von_Mises_or_Tresca
2260  (model &md,
2261  const std::string &datasigma,
2262  const mesh_fem &mf_vm,
2263  model_real_plain_vector &VM,
2264  bool tresca) {
2265 
2266  GMM_ASSERT1(gmm::vect_size(VM) == mf_vm.nb_dof(),
2267  "The vector has not the right size");
2268 
2269  const model_real_plain_vector &sigma_np1 = md.real_variable(datasigma);
2270  const mesh_fem &mf_sigma = md.mesh_fem_of_variable(datasigma);
2271 
2272  // dimension of the finite element used
2273  dim_type N = mf_sigma.linked_mesh().dim();
2274 
2275  GMM_ASSERT1(mf_vm.get_qdim() == 1,
2276  "Target dimension of mf_vm should be 1");
2277 
2278  base_matrix sigma(N, N), Id(N, N);
2279  base_vector eig(N);
2280  base_vector sigma_vm(mf_vm.nb_dof()*N*N);
2281 
2282  gmm::copy(gmm::identity_matrix(), Id);
2283 
2284  interpolation(mf_sigma, mf_vm, sigma_np1, sigma_vm);
2285 
2286  // for each dof we compute the Von Mises or Tresca stress
2287  for (size_type ii = 0; ii < mf_vm.nb_dof(); ++ii) {
2288 
2289  /* we retrieve the matrix sigma_vm on this dof */
2290  std::copy(sigma_vm.begin()+ii*N*N, sigma_vm.begin()+(ii+1)*N*N,
2291  sigma.begin());
2292 
2293  if (!tresca) {
2294  /* von mises: norm(deviator(sigma)) */
2295  gmm::add(gmm::scaled(Id, -gmm::mat_trace(sigma) / N), sigma);
2296 
2297  /* von mises stress=sqrt(3/2)* norm(sigma) */
2298  VM[ii] = sqrt(3.0/2.)*gmm::mat_euclidean_norm(sigma);
2299  } else {
2300  /* else compute the tresca criterion */
2301  gmm::symmetric_qr_algorithm(sigma, eig);
2302  std::sort(eig.begin(), eig.end());
2303  VM[ii] = eig.back() - eig.front();
2304  }
2305  }
2306  }
2307 
2308 
2309 
2310  //=================================================================
2311  // Compute the plastic part
2312  //=================================================================
2313 
2314  void compute_plastic_part(model &md,
2315  const mesh_im &mim,
2316  const mesh_fem &mf_pl,
2317  const std::string &varname,
2318  const std::string &data_previous_disp,
2319  const pconstraints_projection &ACP,
2320  const std::string &datalambda,
2321  const std::string &datamu,
2322  const std::string &datathreshold,
2323  const std::string &datasigma,
2324  model_real_plain_vector &plast) {
2325 
2326  const model_real_plain_vector &u_np1 = md.real_variable(varname);
2327  const model_real_plain_vector &u_n = md.real_variable(data_previous_disp);
2328  const mesh_fem &mf_u = *(md.pmesh_fem_of_variable(varname));
2329 
2330  const model_real_plain_vector &lambda = md.real_variable(datalambda);
2331  const model_real_plain_vector &mu = md.real_variable(datamu);
2332  const model_real_plain_vector &threshold = md.real_variable(datathreshold);
2333  const mesh_fem *pmf_data = md.pmesh_fem_of_variable(datalambda);
2334 
2335  const model_real_plain_vector &sigma_n = md.real_variable(datasigma);
2336  const mesh_fem &mf_sigma = *(md.pmesh_fem_of_variable(datasigma));
2337 
2338  dim_type N = mf_sigma.linked_mesh().dim();
2339 
2340  mesh_region rg = mim.linked_mesh().get_mpi_region();
2341 
2342  model_real_plain_vector dummyV(mf_u.nb_dof());
2343  model_real_plain_vector saved_plast(mf_sigma.nb_dof());
2344  asm_elastoplasticity_rhs(dummyV, &saved_plast,
2345  mim, mf_u, mf_sigma, *pmf_data,
2346  u_n, u_np1, sigma_n,
2347  lambda, mu, threshold, *ACP, PLAST, rg);
2348 
2349  /* Retrieve and save the plastic part */
2350  GMM_ASSERT1(gmm::vect_size(plast) == mf_pl.nb_dof(),
2351  "The vector has not the right size");
2352  GMM_ASSERT1(mf_pl.get_qdim() == 1,
2353  "Target dimension of mf_pl should be 1");
2354 
2355  base_vector saved_pl(mf_pl.nb_dof()*N*N);
2356  interpolation(mf_sigma, mf_pl, saved_plast, saved_pl);
2357 
2358  // for each dof we compute the norm of the plastic part
2359  base_matrix plast_tmp(N, N);
2360  for (size_type ii = 0; ii < mf_pl.nb_dof(); ++ii) {
2361 
2362  /* we retrieve the matrix sigma_pl on this dof */
2363  std::copy(saved_pl.begin()+ii*N*N, saved_pl.begin()+(ii+1)*N*N,
2364  plast_tmp.begin());
2365 
2366  plast[ii] = gmm::mat_euclidean_norm(plast_tmp);
2367 
2368  }
2369 
2370  }
2371 
2372 
2373 
2374 
2375 
2376 } /* end of namespace getfem. */
dal::singleton::instance
static T & instance()
Instance from the current thread.
Definition: dal_singleton.h:165
getfem::im_data
im_data provides indexing to the integration points of a mesh im object.
Definition: getfem_im_data.h:69
getfem::mesh_fem
Describe a finite element method linked to a mesh.
Definition: getfem_mesh_fem.h:152
getfem::mesh_fem::get_qdim
virtual dim_type get_qdim() const
Return the Q dimension.
Definition: getfem_mesh_fem.h:316
getfem::mesh_fem::nb_dof
virtual size_type nb_dof() const
Return the total number of degrees of freedom.
Definition: getfem_mesh_fem.h:566
getfem::mesh_im
Describe an integration method linked to a mesh.
Definition: getfem_mesh_im.h:47
getfem::mesh_region::all_convexes
static mesh_region all_convexes()
provide a default value for the mesh_region parameters of assembly procedures etc.
Definition: getfem_mesh_region.h:142
getfem::model
`‘Model’' variables store the variables, the data and the description of a model.
Definition: getfem_models.h:115
getfem::model::pmesh_fem_of_variable
const mesh_fem * pmesh_fem_of_variable(const std::string &name) const
Gives a pointer to the mesh_fem of a variable if any.
Definition: getfem_models.cc:2895
getfem::model::set_real_variable
model_real_plain_vector & set_real_variable(const std::string &name, size_type niter) const
Gives the write access to the vector value of a variable.
Definition: getfem_models.cc:3011
getfem::model::real_variable
const model_real_plain_vector & real_variable(const std::string &name, size_type niter) const
Gives the access to the vector value of a variable.
Definition: getfem_models.cc:2957
getfem_generic_assembly.h
A language for generic assembly of pde boundary value problems.
getfem_interpolation.h
Interpolation of fields from a mesh_fem onto another.
getfem_models.h
Model representation in Getfem.
getfem_plasticity.h
Plasticty bricks.
gmm::copy
void copy(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:978
gmm::vect_norm2
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norm2(const V &v)
Euclidean norm of a vector.
Definition: gmm_blas.h:558
gmm::mat_trace
linalg_traits< M >::value_type mat_trace(const M &m)
Trace of a matrix.
Definition: gmm_blas.h:529
gmm::mat_euclidean_norm
number_traits< typename linalg_traits< M >::value_type >::magnitude_type mat_euclidean_norm(const M &m)
Euclidean norm of a matrix.
Definition: gmm_blas.h:637
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:59
gmm::resize
void resize(V &v, size_type n)
*‍/
Definition: gmm_blas.h:210
gmm::mult
void mult(const L1 &l1, const L2 &l2, L3 &l3)
*‍/
Definition: gmm_blas.h:1664
gmm::add
void add(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:1277
gmm::lu_inverse
void lu_inverse(const DenseMatrixLU &LU, const Pvector &pvector, const DenseMatrix &AInv_)
Given an LU factored matrix, build the inverse of the matrix.
Definition: gmm_dense_lu.h:212
gmm_dense_matrix_functions.h
Common matrix functions for dense matrices.
gmm::logm
void logm(const dense_matrix< T > &A, dense_matrix< T > &LOGMA)
Matrix logarithm (from GNU/Octave)
Definition: gmm_dense_matrix_functions.h:252
getfem::pfem
std::shared_ptr< const getfem::virtual_fem > pfem
type of pointer on a fem description
Definition: getfem_fem.h:244
bgeot
Basic Geometric Tools.
Definition: bgeot_convex_ref.cc:27
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:73
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
bgeot::pgeometric_trans
std::shared_ptr< const bgeot::geometric_trans > pgeometric_trans
pointer type for a geometric transformation
Definition: bgeot_geometric_trans.h:186
bgeot::alpha
size_type alpha(short_type n, short_type d)
Return the value of which is the number of monomials of a polynomial of variables and degree .
Definition: bgeot_poly.cc:47
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:47
getfem::add_nonlinear_term
size_type APIDECL add_nonlinear_term(model &md, const mesh_im &mim, const std::string &expr, size_type region=size_type(-1), bool is_sym=false, bool is_coercive=false, const std::string &brickname="")
Add a nonlinear term given by the weak form language expression expr which will be assembled in regio...
Definition: getfem_models.cc:3676
getfem::ga_define_linear_hardening_function
void ga_define_linear_hardening_function(const std::string &name, scalar_type sigma_y0, scalar_type H, bool frobenius=true)
Add a linear function with the name specified by name to represent linear isotropoc hardening in plas...
Definition: getfem_plasticity.cc:1420
getfem::compute_plastic_part
void compute_plastic_part(model &md, const mesh_im &mim, const mesh_fem &mf_pl, const std::string &varname, const std::string &previous_dep_name, const pconstraints_projection &ACP, const std::string &datalambda, const std::string &datamu, const std::string &datathreshold, const std::string &datasigma, model_real_plain_vector &plast)
This function computes on mf_pl the plastic part, that could appear after a load and an unload,...
getfem::compute_small_strain_elastoplasticity_Von_Mises
void compute_small_strain_elastoplasticity_Von_Mises(model &md, const mesh_im &mim, std::string lawname, plasticity_unknowns_type unknowns_type, const std::vector< std::string > &varnames, const std::vector< std::string > &params, const mesh_fem &mf_vm, model_real_plain_vector &VM, size_type region=size_type(-1))
This function computes the Von Mises stress field with respect to a small strain elastoplasticity ter...
Definition: getfem_plasticity.cc:1373
getfem::finite_strain_elastoplasticity_next_iter
void finite_strain_elastoplasticity_next_iter(model &md, const mesh_im &mim, std::string lawname, plasticity_unknowns_type unknowns_type, const std::vector< std::string > &varnames, const std::vector< std::string > &params, size_type region=size_type(-1))
This function permits to update the state variables for a finite strain elastoplasticity brick,...
getfem::add_finite_strain_elastoplasticity_brick
size_type add_finite_strain_elastoplasticity_brick(model &md, const mesh_im &mim, std::string lawname, plasticity_unknowns_type unknowns_type, const std::vector< std::string > &varnames, const std::vector< std::string > &params, size_type region=size_type(-1))
Add a finite strain elastoplasticity brick to the model.
getfem::elastoplasticity_next_iter
void elastoplasticity_next_iter(model &md, const mesh_im &mim, const std::string &varname, const std::string &previous_dep_name, const pconstraints_projection &ACP, const std::string &datalambda, const std::string &datamu, const std::string &datathreshold, const std::string &datasigma)
This function permits to compute the new stress constraints values supported by the material after a ...
getfem::interpolation
void interpolation(const mesh_fem &mf_source, const mesh_fem &mf_target, const VECTU &U, VECTV &V, int extrapolation=0, double EPS=1E-10, mesh_region rg_source=mesh_region::all_convexes(), mesh_region rg_target=mesh_region::all_convexes())
interpolation/extrapolation of (mf_source, U) on mf_target.
Definition: getfem_interpolation.h:693
getfem::compute_elastoplasticity_Von_Mises_or_Tresca
void compute_elastoplasticity_Von_Mises_or_Tresca(model &md, const std::string &datasigma, const mesh_fem &mf_vm, model_real_plain_vector &VM, bool tresca)
This function computes on mf_vm the Von Mises or Tresca stress of a field for elastoplasticity and re...
getfem::pbrick
std::shared_ptr< const virtual_brick > pbrick
type of pointer on a brick
Definition: getfem_models.h:49
getfem::ga_define_Ramberg_Osgood_hardening_function
void ga_define_Ramberg_Osgood_hardening_function(const std::string &name, scalar_type sigma_ref, scalar_type eps_ref, scalar_type n, bool frobenius=false)
Add a Ramberg-Osgood hardening function with the name specified by name, for reference stress and str...
Definition: getfem_plasticity.cc:1433
getfem::add_small_strain_elastoplasticity_brick
size_type add_small_strain_elastoplasticity_brick(model &md, const mesh_im &mim, std::string lawname, plasticity_unknowns_type unknowns_type, const std::vector< std::string > &varnames, const std::vector< std::string > &params, size_type region=size_type(-1))
Adds a small strain plasticity term to the model md.
Definition: getfem_plasticity.cc:1270
getfem::ga_local_projection
void ga_local_projection(const getfem::model &md, const mesh_im &mim, const std::string &expr, const mesh_fem &mf, base_vector &result, const mesh_region &rg=mesh_region::all_convexes())
Make an elementwise L2 projection of an expression with respect to the mesh_fem mf.
Definition: getfem_generic_assembly_interpolation.cc:437
getfem::compute_finite_strain_elastoplasticity_Von_Mises
void compute_finite_strain_elastoplasticity_Von_Mises(model &md, const mesh_im &mim, std::string lawname, plasticity_unknowns_type unknowns_type, const std::vector< std::string > &varnames, const std::vector< std::string > &params, const mesh_fem &mf_vm, model_real_plain_vector &VM, size_type region=size_type(-1))
This function computes the Von Mises stress field with respect to a finite strain elastoplasticity te...
getfem::small_strain_elastoplasticity_next_iter
void small_strain_elastoplasticity_next_iter(model &md, const mesh_im &mim, std::string lawname, plasticity_unknowns_type unknowns_type, const std::vector< std::string > &varnames, const std::vector< std::string > &params, size_type region=size_type(-1))
Function that allows to pass from a time step to another for the small strain plastic brick.
Definition: getfem_plasticity.cc:1301
getfem::add_elastoplasticity_brick
size_type add_elastoplasticity_brick(model &md, const mesh_im &mim, const pconstraints_projection &ACP, const std::string &varname, const std::string &previous_dep_name, const std::string &datalambda, const std::string &datamu, const std::string &datathreshold, const std::string &datasigma, size_type region=size_type(-1))
Add a nonlinear elastoplasticity term to the model for small deformations and an isotropic material,...