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