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