GetFEM  5.4.2
getfem_nonlinear_elasticity.cc
1 /*===========================================================================
2 
3  Copyright (C) 2000-2020 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_nonlinear_elasticity.h"
25 #include "getfem/getfem_generic_assembly.h"
26 
27 namespace getfem {
28 
29 
30  /* Useful functions to compute the invariants and their derivatives
31  Note that the second derivative is symmetrized (see the user
32  documentation for more details). The matrix E is assumed to be symmetric.
33  */
34 
35 
36  static scalar_type frobenius_product_trans(const base_matrix &A,
37  const base_matrix &B) {
38  size_type N = gmm::mat_nrows(A);
39  scalar_type res = scalar_type(0);
40  for (size_type i = 0; i < N; ++i)
41  for (size_type j = 0; j < N; ++j)
42  res += A(i, j) * B(j, i);
43  return res;
44  }
45 
46  struct compute_invariants {
47 
48  const base_matrix &E;
49  base_matrix Einv;
50  size_type N;
51  scalar_type i1_, i2_, i3_, j1_, j2_;
52  bool i1_c, i2_c, i3_c, j1_c, j2_c;
53 
54  base_matrix di1, di2, di3, dj1, dj2;
55  bool di1_c, di2_c, di3_c, dj1_c, dj2_c;
56 
57  base_tensor ddi1, ddi2, ddi3, ddj1, ddj2;
58  bool ddi1_c, ddi2_c, ddi3_c, ddj1_c, ddj2_c;
59 
60 
61  /* First invariant tr(E) */
62  void compute_i1() {
63  i1_ = gmm::mat_trace(E);
64  i1_c = true;
65  }
66 
67  void compute_di1() {
68  gmm::resize(di1, N, N);
69  gmm::copy(gmm::identity_matrix(), di1);
70  di1_c = true;
71  }
72 
73  void compute_ddi1() { // not very useful, null tensor
74  ddi1 = base_tensor(N, N, N, N);
75  ddi1_c = true;
76  }
77 
78  inline scalar_type i1()
79  { if (!i1_c) compute_i1(); return i1_; }
80 
81  inline const base_matrix &grad_i1()
82  { if (!di1_c) compute_di1(); return di1; }
83 
84  inline const base_tensor &sym_grad_grad_i1()
85  { if (!ddi1_c) compute_ddi1(); return ddi1; }
86 
87 
88  /* Second invariant (tr(E)^2 - tr(E^2))/2 */
89  void compute_i2() {
90  i2_ = (gmm::sqr(gmm::mat_trace(E))
91  - frobenius_product_trans(E, E)) / scalar_type(2);
92  i2_c = true;
93  }
94 
95  void compute_di2() {
96  gmm::resize(di2, N, N);
97  gmm::copy(gmm::identity_matrix(), di2);
98  gmm::scale(di2, i1());
99  // gmm::add(gmm::scale(gmm::transposed(E), -scalar_type(1)), di2);
100  gmm::add(gmm::scaled(E, -scalar_type(1)), di2);
101  di2_c = true;
102  }
103 
104  void compute_ddi2() {
105  ddi2 = base_tensor(N, N, N, N);
106  for (size_type i = 0; i < N; ++i)
107  for (size_type k = 0; k < N; ++k)
108  ddi2(i,i,k,k) += scalar_type(1);
109  for (size_type i = 0; i < N; ++i)
110  for (size_type j = 0; j < N; ++j) {
111  ddi2(i,j,j,i) -= scalar_type(1)/scalar_type(2);
112  ddi2(j,i,j,i) -= scalar_type(1)/scalar_type(2);
113  }
114  ddi2_c = true;
115  }
116 
117  inline scalar_type i2()
118  { if (!i2_c) compute_i2(); return i2_; }
119 
120  inline const base_matrix &grad_i2()
121  { if (!di2_c) compute_di2(); return di2; }
122 
123  inline const base_tensor &sym_grad_grad_i2()
124  { if (!ddi2_c) compute_ddi2(); return ddi2; }
125 
126  /* Third invariant det(E) */
127  void compute_i3() {
128  Einv = E;
129  i3_ = bgeot::lu_inverse(&(*(Einv.begin())), gmm::mat_nrows(Einv));
130  i3_c = true;
131  }
132 
133  void compute_di3() {
134  scalar_type det = i3();
135  // gmm::resize(di3, N, N);
136  // gmm::copy(gmm::transposed(E), di3);
137  di3 = Einv;
138  // gmm::lu_inverse(di3);
139  gmm::scale(di3, det);
140  di3_c = true;
141  }
142 
143  void compute_ddi3() {
144  ddi3 = base_tensor(N, N, N, N);
145  scalar_type det = i3() / scalar_type(2); // computes also E inverse.
146  for (size_type i = 0; i < N; ++i)
147  for (size_type j = 0; j < N; ++j)
148  for (size_type k = 0; k < N; ++k)
149  for (size_type l = 0; l < N; ++l)
150  ddi3(i,j,k,l) = det*(Einv(j,i)*Einv(l,k) - Einv(j,k)*Einv(l,i)
151  + Einv(i,j)*Einv(l,k) - Einv(i,k)*Einv(l,j));
152  ddi3_c = true;
153  }
154 
155  inline scalar_type i3()
156  { if (!i3_c) compute_i3(); return i3_; }
157 
158  inline const base_matrix &grad_i3()
159  { if (!di3_c) compute_di3(); return di3; }
160 
161  inline const base_tensor &sym_grad_grad_i3()
162  { if (!ddi3_c) compute_ddi3(); return ddi3; }
163 
164  /* Invariant j1(E) = i1(E)*i3(E)^(-1/3) */
165  void compute_j1() {
166  j1_ = i1() * ::pow(gmm::abs(i3()), -scalar_type(1) / scalar_type(3));
167  j1_c = true;
168  }
169 
170  void compute_dj1() {
171  dj1 = grad_i1();
172  gmm::add(gmm::scaled(grad_i3(), -i1() / (scalar_type(3) * i3())), dj1);
173  gmm::scale(dj1, ::pow(gmm::abs(i3()), -scalar_type(1) / scalar_type(3)));
174  dj1_c = true;
175  }
176 
177  void compute_ddj1() {
178  const base_matrix &di1_ = grad_i1();
179  const base_matrix &di3_ = grad_i3();
180  scalar_type coeff1 = scalar_type(1) / (scalar_type(3)*i3());
181  scalar_type coeff2 = scalar_type(4) * coeff1 * coeff1 * i1();
182  ddj1 = sym_grad_grad_i3();
183  gmm::scale(ddj1.as_vector(), -i1() * coeff1);
184 
185  for (size_type i = 0; i < N; ++i)
186  for (size_type j = 0; j < N; ++j)
187  for (size_type k = 0; k < N; ++k)
188  for (size_type l = 0; l < N; ++l)
189  ddj1(i,j,k,l) +=
190  (di3_(i, j) * di3_(k, l)) * coeff2
191  - (di1_(i, j) * di3_(k, l) + di1_(k, l) * di3_(i, j)) * coeff1;
192 
193  gmm::scale(ddj1.as_vector(),
194  ::pow(gmm::abs(i3()), -scalar_type(1)/scalar_type(3)));
195  ddj1_c = true;
196  }
197 
198  inline scalar_type j1()
199  { if (!j1_c) compute_j1(); return j1_; }
200 
201  inline const base_matrix &grad_j1()
202  { if (!dj1_c) compute_dj1(); return dj1; }
203 
204  inline const base_tensor &sym_grad_grad_j1()
205  { if (!ddj1_c) compute_ddj1(); return ddj1; }
206 
207  /* Invariant j2(E) = i2(E)*i3(E)^(-2/3) */
208  void compute_j2() {
209  j2_ = i2() * ::pow(gmm::abs(i3()), -scalar_type(2) / scalar_type(3));
210  j2_c = true;
211  }
212 
213  void compute_dj2() {
214  dj2 = grad_i2();
215  gmm::add(gmm::scaled(grad_i3(), -scalar_type(2) * i2() / (scalar_type(3) * i3())), dj2);
216  gmm::scale(dj2, ::pow(gmm::abs(i3()), -scalar_type(2) / scalar_type(3)));
217  dj2_c = true;
218  }
219 
220  void compute_ddj2() {
221  const base_matrix &di2_ = grad_i2();
222  const base_matrix &di3_ = grad_i3();
223  scalar_type coeff1 = scalar_type(2) / (scalar_type(3)*i3());
224  scalar_type coeff2 = scalar_type(5)*coeff1*coeff1*i2() / scalar_type(2);
225  ddj2 = sym_grad_grad_i2();
226  gmm::add(gmm::scaled(sym_grad_grad_i3().as_vector(), -i2() * coeff1),
227  ddj2.as_vector());
228 
229  for (size_type i = 0; i < N; ++i)
230  for (size_type j = 0; j < N; ++j)
231  for (size_type k = 0; k < N; ++k)
232  for (size_type l = 0; l < N; ++l)
233  ddj2(i,j,k,l) +=
234  (di3_(i, j) * di3_(k, l)) * coeff2
235  - (di2_(i, j) * di3_(k, l) + di2_(k, l) * di3_(i, j)) * coeff1;
236 
237  gmm::scale(ddj2.as_vector(),
238  ::pow(gmm::abs(i3()), -scalar_type(2)/scalar_type(3)));
239  ddj2_c = true;
240  }
241 
242 
243  inline scalar_type j2()
244  { if (!j2_c) compute_j2(); return j2_; }
245 
246  inline const base_matrix &grad_j2()
247  { if (!dj2_c) compute_dj2(); return dj2; }
248 
249  inline const base_tensor &sym_grad_grad_j2()
250  { if (!ddj2_c) compute_ddj2(); return ddj2; }
251 
252 
253  compute_invariants(const base_matrix &EE)
254  : E(EE), i1_c(false), i2_c(false), i3_c(false),
255  j1_c(false), j2_c(false), di1_c(false), di2_c(false), di3_c(false),
256  dj1_c(false), dj2_c(false), ddi1_c(false), ddi2_c(false),
257  ddi3_c(false), ddj1_c(false), ddj2_c(false)
258  { N = gmm::mat_nrows(E); }
259 
260  };
261 
262 
263 
264 
265 
266  /* Symmetry check */
267 
268  int check_symmetry(const base_tensor &t) {
269  int flags = 7; size_type N = 3;
270  for (size_type n = 0; n < N; ++n)
271  for (size_type m = 0; m < N; ++m)
272  for (size_type l = 0; l < N; ++l)
273  for (size_type k = 0; k < N; ++k) {
274  if (gmm::abs(t(n,m,l,k) - t(l,k,n,m))>1e-5) flags &= (~1);
275  if (gmm::abs(t(n,m,l,k) - t(m,n,l,k))>1e-5) flags &= (~2);
276  if (gmm::abs(t(n,m,l,k) - t(n,m,k,l))>1e-5) flags &= (~4);
277  }
278  return flags;
279  }
280 
281  /* Member functions of hyperelastic laws */
282 
283  void abstract_hyperelastic_law::random_E(base_matrix &E) {
284  size_type N = gmm::mat_nrows(E);
285  base_matrix Phi(N,N);
286  scalar_type d;
287  do {
288  gmm::fill_random(Phi);
289  d = bgeot::lu_det(&(*(Phi.begin())), N);
290  } while (d < scalar_type(0.01));
291  gmm::mult(gmm::transposed(Phi),Phi,E);
292  gmm::scale(E,-1.); gmm::add(gmm::identity_matrix(),E);
293  gmm::scale(E,-0.5);
294  }
295 
296  void abstract_hyperelastic_law::test_derivatives
297  (size_type N, scalar_type h, const base_vector& param) const {
298  base_matrix E(N,N), E2(N,N), DE(N,N);
299  bool ok = true;
300 
301  for (size_type count = 0; count < 100; ++count) {
302  random_E(E); random_E(DE);
303  gmm::scale(DE, h);
304  gmm::add(E, DE, E2);
305 
306  base_matrix sigma1(N,N), sigma2(N,N);
307  getfem::base_tensor tdsigma(N,N,N,N);
308  base_matrix dsigma(N,N);
309  gmm::copy(E, E2); gmm::add(DE, E2);
310  sigma(E, sigma1, param, scalar_type(1));
311  sigma(E2, sigma2, param, scalar_type(1));
312 
313  scalar_type d = strain_energy(E2, param, scalar_type(1))
314  - strain_energy(E, param, scalar_type(1));
315  scalar_type d2 = 0;
316  for (size_type i=0; i < N; ++i)
317  for (size_type j=0; j < N; ++j) d2 += sigma1(i,j)*DE(i,j);
318  if (gmm::abs(d-d2)/(gmm::abs(d)+1e-40) > 1e-4) {
319  cout << "Test " << count << " wrong derivative of strain_energy, d="
320  << d/h << ", d2=" << d2/h << endl;
321  ok = false;
322  }
323 
324  grad_sigma(E,tdsigma,param, scalar_type(1));
325  for (size_type i=0; i < N; ++i) {
326  for (size_type j=0; j < N; ++j) {
327  dsigma(i,j) = 0;
328  for (size_type k=0; k < N; ++k) {
329  for (size_type m=0; m < N; ++m) {
330  dsigma(i,j) += tdsigma(i,j,k,m)*DE(k,m);
331  }
332  }
333  sigma2(i,j) -= sigma1(i,j);
334  if (gmm::abs(dsigma(i,j) - sigma2(i,j))
335  /(gmm::abs(dsigma(i,j)) + 1e-40) > 1.5e-4) {
336  cout << "Test " << count << " wrong derivative of sigma, i="
337  << i << ", j=" << j << ", dsigma=" << dsigma(i,j)/h
338  << ", var sigma = " << sigma2(i,j)/h << endl;
339  ok = false;
340  }
341  }
342  }
343  }
344  GMM_ASSERT1(ok, "Derivative test has failed");
345  }
346 
347  void abstract_hyperelastic_law::cauchy_updated_lagrangian
348  (const base_matrix& F, const base_matrix &E,
349  base_matrix &cauchy_stress, const base_vector &params,
350  scalar_type det_trans) const
351  {
352  size_type N = E.ncols();
353  base_matrix PK2(N,N);
354  sigma(E,PK2,params,det_trans);//second Piola-Kirchhoff stress
355  base_matrix aux(N,N);
356  gmm::mult(F,PK2,aux);
357  gmm::mult(aux,gmm::transposed(F),cauchy_stress);
358  gmm::scale(cauchy_stress,scalar_type(1.0/det_trans)); //cauchy = 1/J*F*PK2*F^T
359  }
360 
361 
362  void abstract_hyperelastic_law::grad_sigma_updated_lagrangian
363  (const base_matrix& F, const base_matrix& E,
364  const base_vector &params, scalar_type det_trans,
365  base_tensor &grad_sigma_ul) const
366  {
367  size_type N = E.ncols();
368  base_tensor Cse(N,N,N,N);
369  grad_sigma(E,Cse,params,det_trans);
370  scalar_type mult = 1.0/det_trans;
371  // this is a general transformation for an anisotropic material, very non-efficient;
372  // more effiecient calculations can be overloaded for every specific material
373  for(size_type i = 0; i < N; ++i)
374  for(size_type j = 0; j < N; ++j)
375  for(size_type k = 0; k < N; ++k)
376  for(size_type l = 0; l < N; ++l)
377  {
378  grad_sigma_ul(i,j,k,l) = 0.0;
379  for(size_type m = 0; m < N; ++m)
380  { for(size_type n = 0; n < N; ++n)
381  for(size_type p = 0; p < N; ++p)
382  for(size_type q = 0; q < N; ++q)
383  grad_sigma_ul(i,j,k,l)+=
384  F(i,m)*F(j,n)*F(k,p)*F(l,q)*Cse(m,n,p,q);
385  }
386  grad_sigma_ul(i,j,k,l) *= mult;
387  }
388  }
389 
390  scalar_type SaintVenant_Kirchhoff_hyperelastic_law::strain_energy
391  (const base_matrix &E, const base_vector &params, scalar_type det_trans) const {
392  // should be optimized, maybe deriving sigma from strain energy
393  if (det_trans <= scalar_type(0))
394  { return 1e200; }
395 
396  return gmm::sqr(gmm::mat_trace(E)) * params[0] / scalar_type(2)
397  + gmm::mat_euclidean_norm_sqr(E) * params[1];
398  }
399 
400  void SaintVenant_Kirchhoff_hyperelastic_law::sigma
401  (const base_matrix &E, base_matrix &result,const base_vector &params, scalar_type det_trans) const {
402  gmm::copy(gmm::identity_matrix(), result);
403  gmm::scale(result, params[0] * gmm::mat_trace(E));
404  gmm::add(gmm::scaled(E, 2 * params[1]), result);
405  if (det_trans <= scalar_type(0)) {
406  gmm::add(gmm::scaled(E, 1e200), result);
407  }
408  }
409  void SaintVenant_Kirchhoff_hyperelastic_law::grad_sigma
410  (const base_matrix &E, base_tensor &result,const base_vector &params, scalar_type) const {
411  std::fill(result.begin(), result.end(), scalar_type(0));
412  size_type N = gmm::mat_nrows(E);
413  for (size_type i = 0; i < N; ++i)
414  for (size_type l = 0; l < N; ++l) {
415  result(i, i, l, l) += params[0];
416  result(i, l, i, l) += params[1]/scalar_type(2);
417  result(i, l, l, i) += params[1]/scalar_type(2);
418  result(l, i, i, l) += params[1]/scalar_type(2);
419  result(l, i, l, i) += params[1]/scalar_type(2);
420  }
421  }
422 
423  void SaintVenant_Kirchhoff_hyperelastic_law::grad_sigma_updated_lagrangian(const base_matrix& F,
424  const base_matrix& E,
425  const base_vector &params,
426  scalar_type det_trans,
427  base_tensor &grad_sigma_ul)const
428  {
429  size_type N = E.ncols();
430  base_tensor Cse(N,N,N,N);
431  grad_sigma(E,Cse,params,det_trans);
432  base_matrix Cinv(N,N); // left Cauchy-Green deform. tens.
433  gmm::mult(F,gmm::transposed(F),Cinv);
434  scalar_type mult=1.0/det_trans;
435  for(size_type i = 0; i < N; ++i)
436  for(size_type j = 0; j < N; ++j)
437  for(size_type k = 0; k < N; ++k)
438  for(size_type l = 0; l < N; ++l)
439  grad_sigma_ul(i, j, k, l)= (Cinv(i,j)*Cinv(k,l)*params[0] +
440  params[1]*(Cinv(i,k)*Cinv(j,l) + Cinv(i,l)*Cinv(j,k)))*mult;
441  }
442 
443  SaintVenant_Kirchhoff_hyperelastic_law::SaintVenant_Kirchhoff_hyperelastic_law() {
444  nb_params_ = 2;
445  }
446 
447  scalar_type membrane_elastic_law::strain_energy
448  (const base_matrix & /* E */, const base_vector & /* params */, scalar_type) const {
449  // to be done if needed
450  GMM_ASSERT1(false, "To be done");
451  return 0;
452  }
453 
454  void membrane_elastic_law::sigma
455  (const base_matrix &E, base_matrix &result,const base_vector &params, scalar_type det_trans) const {
456  // should be optimized, maybe deriving sigma from strain energy
457  base_tensor tt(2,2,2,2);
458  size_type N = (gmm::mat_nrows(E) > 2)? 2 : gmm::mat_nrows(E);
459  grad_sigma(E,tt,params, det_trans);
460  for (size_type i = 0; i < N; ++i)
461  for (size_type j = 0; j < N; ++j) {
462  result(i,j)=0.0;
463  for (size_type k = 0; k < N; ++k)
464  for (size_type l = 0; l < N; ++l)
465  result(i,j)+=tt(i,j,k,l)*E(k,l);
466  }
467  // add pretension in X' direction
468  if(params[4]!=0) result(0,0)+=params[4];
469  // add pretension in Y' direction
470  if(params[5]!=0) result(1,1)+=params[5];
471  // cout<<"sigma="<<result<<endl;
472  }
473 
474  void membrane_elastic_law::grad_sigma
475  (const base_matrix & /* E */, base_tensor &result,
476  const base_vector &params, scalar_type) const {
477  // to be optimized!!
478  std::fill(result.begin(), result.end(), scalar_type(0));
479  scalar_type poisonXY=params[0]*params[1]/params[2]; //Ex*vYX=Ey*vXY
480  //if no G entered, compute G=E/(2*(1+v))
481  scalar_type G=( params[3] == 0) ? params[0]/(2*(1+params[1])) : params[3];
482  std::fill(result.begin(), result.end(), scalar_type(0));
483  result(0,0,0,0) = params[0]/(1-params[1]*poisonXY);
484  // result(0,0,0,1) = 0;
485  // result(0,0,1,0) = 0;
486  result(0,0,1,1) = params[1]*params[0]/(1-params[1]*poisonXY);
487  result(1,1,0,0) = params[1]*params[0]/(1-params[1]*poisonXY);
488  // result(1,1,0,1) = 0;
489  // result(1,1,1,0) = 0;
490  result(1,1,1,1) = params[2]/(1-params[1]*poisonXY);
491  // result(0,1,0,0) = 0;
492  result(0,1,0,1) = G;
493  result(0,1,1,0) = G;
494  // result(0,1,1,1) = 0;
495  // result(1,0,0,0) = 0;
496  result(1,0,0,1) = G;
497  result(1,0,1,0) = G;
498  // result(1,0,1,1) = 0;
499  }
500 
501  scalar_type Mooney_Rivlin_hyperelastic_law::strain_energy
502  (const base_matrix &E, const base_vector &params
503  ,scalar_type det_trans) const {
504  //shouldn't negative det_trans be handled here???
505  if (compressible && det_trans <= scalar_type(0)) return 1e200;
506  size_type N = gmm::mat_nrows(E);
507  GMM_ASSERT1(N == 3, "Mooney Rivlin hyperelastic law only defined "
508  "on dimension 3, sorry");
509  base_matrix C = E;
510  gmm::scale(C, scalar_type(2));
511  gmm::add(gmm::identity_matrix(), C);
512  compute_invariants ci(C);
513  size_type i=0;
514  scalar_type C1 = params[i++]; // C10
515  scalar_type W = C1 * (ci.j1() - scalar_type(3));
516  if (!neohookean) {
517  scalar_type C2 = params[i++]; // C01
518  W += C2 * (ci.j2() - scalar_type(3));
519  }
520  if (compressible) {
521  scalar_type D1 = params[i++];
522  W += D1 * gmm::sqr(sqrt(gmm::abs(ci.i3())) - scalar_type(1));
523  }
524  return W;
525  }
526 
527  void Mooney_Rivlin_hyperelastic_law::sigma
528  (const base_matrix &E, base_matrix &result,
529  const base_vector &params, scalar_type det_trans) const {
530  size_type N = gmm::mat_nrows(E);
531  GMM_ASSERT1(N == 3, "Mooney Rivlin hyperelastic law only defined "
532  "on dimension 3, sorry");
533  base_matrix C = E;
534  gmm::scale(C, scalar_type(2));
535  gmm::add(gmm::identity_matrix(), C);
536  compute_invariants ci(C);
537 
538  size_type i=0;
539  scalar_type C1 = params[i++]; // C10
540  gmm::copy(gmm::scaled(ci.grad_j1(), scalar_type(2) * C1), result);
541  if (!neohookean) {
542  scalar_type C2 = params[i++]; // C01
543  gmm::add(gmm::scaled(ci.grad_j2(), scalar_type(2) * C2), result);
544  }
545  if (compressible) {
546  scalar_type D1 = params[i++];
547  scalar_type di3 = D1 - D1 / sqrt(gmm::abs(ci.i3()));
548  gmm::add(gmm::scaled(ci.grad_i3(), scalar_type(2) * di3), result);
549 
550 // shouldn't negative det_trans be handled here???
551  if (det_trans <= scalar_type(0))
552  gmm::add(gmm::scaled(C, 1e200), result);
553  }
554  }
555 
556  void Mooney_Rivlin_hyperelastic_law::grad_sigma
557  (const base_matrix &E, base_tensor &result,
558  const base_vector &params, scalar_type) const {
559  size_type N = gmm::mat_nrows(E);
560  GMM_ASSERT1(N == 3, "Mooney Rivlin hyperelastic law only defined "
561  "on dimension 3, sorry");
562  base_matrix C = E;
563  gmm::scale(C, scalar_type(2));
564  gmm::add(gmm::identity_matrix(), C);
565  compute_invariants ci(C);
566 
567  size_type i=0;
568  scalar_type C1 = params[i++]; // C10
569  gmm::copy(gmm::scaled(ci.sym_grad_grad_j1().as_vector(),
570  scalar_type(4)*C1), result.as_vector());
571  if (!neohookean) {
572  scalar_type C2 = params[i++]; // C01
573  gmm::add(gmm::scaled(ci.sym_grad_grad_j2().as_vector(),
574  scalar_type(4)*C2), result.as_vector());
575  }
576  if (compressible) {
577  scalar_type D1 = params[i++];
578  scalar_type di3 = D1 - D1 / sqrt(gmm::abs(ci.i3()));
579  gmm::add(gmm::scaled(ci.sym_grad_grad_i3().as_vector(),
580  scalar_type(4)*di3), result.as_vector());
581 
582  // second derivatives of W with respect to the third invariant
583  scalar_type A22 = D1 / (scalar_type(2) * pow(gmm::abs(ci.i3()), 1.5));
584  const base_matrix &di = ci.grad_i3();
585  for (size_type l1 = 0; l1 < N; ++l1)
586  for (size_type l2 = 0; l2 < N; ++l2)
587  for (size_type l3 = 0; l3 < N; ++l3)
588  for (size_type l4 = 0; l4 < N; ++l4)
589  result(l1, l2, l3, l4) +=
590  scalar_type(4) * A22 * di(l1, l2) * di(l3, l4);
591  }
592 
593 // GMM_ASSERT1(check_symmetry(result) == 7,
594 // "Fourth order tensor not symmetric : " << result);
595  }
596 
597  Mooney_Rivlin_hyperelastic_law::Mooney_Rivlin_hyperelastic_law
598  (bool compressible_, bool neohookean_)
599  : compressible(compressible_), neohookean(neohookean_)
600  {
601  nb_params_ = 2;
602  if (compressible) ++nb_params_; // D1 != 0
603  if (neohookean) --nb_params_; // C2 == 0
604 
605  }
606 
607 
608 
609 
610  scalar_type Neo_Hookean_hyperelastic_law::strain_energy
611  (const base_matrix &E, const base_vector &params, scalar_type det_trans) const {
612  if (det_trans <= scalar_type(0)) return 1e200;
613  size_type N = gmm::mat_nrows(E);
614  GMM_ASSERT1(N == 3, "Neo Hookean hyperelastic law only defined "
615  "on dimension 3, sorry");
616  base_matrix C = E;
617  gmm::scale(C, scalar_type(2));
618  gmm::add(gmm::identity_matrix(), C);
619  compute_invariants ci(C);
620 
621  scalar_type lambda = params[0];
622  scalar_type mu = params[1];
623  scalar_type logi3 = log(ci.i3());
624  scalar_type W = mu/2 * (ci.i1() - scalar_type(3) - logi3);
625  if (bonet)
626  W += lambda/8 * gmm::sqr(logi3);
627  else // Wriggers
628  W += lambda/4 * (ci.i3() - scalar_type(1) - logi3);
629 
630  return W;
631  }
632 
633  void Neo_Hookean_hyperelastic_law::sigma
634  (const base_matrix &E, base_matrix &result,
635  const base_vector &params , scalar_type det_trans) const {
636  size_type N = gmm::mat_nrows(E);
637  GMM_ASSERT1(N == 3, "Neo Hookean hyperelastic law only defined "
638  "on dimension 3, sorry");
639  base_matrix C = E;
640  gmm::scale(C, scalar_type(2));
641  gmm::add(gmm::identity_matrix(), C);
642  compute_invariants ci(C);
643 
644  scalar_type lambda = params[0];
645  scalar_type mu = params[1];
646  gmm::copy(gmm::scaled(ci.grad_i1(), mu), result);
647  if (bonet)
648  gmm::add(gmm::scaled(ci.grad_i3(),
649  (lambda/2 * log(ci.i3()) - mu) / ci.i3()), result);
650  else
651  gmm::add(gmm::scaled(ci.grad_i3(),
652  lambda/2 - lambda/(2*ci.i3()) - mu / ci.i3()), result);
653  if (det_trans <= scalar_type(0))
654  gmm::add(gmm::scaled(C, 1e200), result);
655  }
656 
657  void Neo_Hookean_hyperelastic_law::grad_sigma
658  (const base_matrix &E, base_tensor &result,
659  const base_vector &params, scalar_type) const {
660  size_type N = gmm::mat_nrows(E);
661  GMM_ASSERT1(N == 3, "Neo Hookean hyperelastic law only defined "
662  "on dimension 3, sorry");
663  base_matrix C = E;
664  gmm::scale(C, scalar_type(2));
665  gmm::add(gmm::identity_matrix(), C);
666  compute_invariants ci(C);
667 
668  scalar_type lambda = params[0];
669  scalar_type mu = params[1];
670 
671  scalar_type coeff;
672  if (bonet) {
673  scalar_type logi3 = log(ci.i3());
674  gmm::copy(gmm::scaled(ci.sym_grad_grad_i3().as_vector(),
675  (lambda * logi3 - 2*mu) / ci.i3()),
676  result.as_vector());
677  coeff = (lambda + 2 * mu - lambda * logi3) / gmm::sqr(ci.i3());
678  } else {
679  gmm::copy(gmm::scaled(ci.sym_grad_grad_i3().as_vector(),
680  lambda - (lambda + 2 * mu) / ci.i3()),
681  result.as_vector());
682  coeff = (lambda + 2 * mu) / gmm::sqr(ci.i3());
683  }
684 
685  const base_matrix &di = ci.grad_i3();
686  for (size_type l1 = 0; l1 < N; ++l1)
687  for (size_type l2 = 0; l2 < N; ++l2)
688  for (size_type l3 = 0; l3 < N; ++l3)
689  for (size_type l4 = 0; l4 < N; ++l4)
690  result(l1, l2, l3, l4) += coeff * di(l1, l2) * di(l3, l4);
691 
692 // GMM_ASSERT1(check_symmetry(result) == 7,
693 // "Fourth order tensor not symmetric : " << result);
694  }
695 
696  Neo_Hookean_hyperelastic_law::Neo_Hookean_hyperelastic_law(bool bonet_)
697  : bonet(bonet_)
698  {
699  nb_params_ = 2;
700  }
701 
702 
703 
704  scalar_type generalized_Blatz_Ko_hyperelastic_law::strain_energy
705  (const base_matrix &E, const base_vector &params, scalar_type det_trans) const {
706  if (det_trans <= scalar_type(0)) return 1e200;
707  scalar_type a = params[0], b = params[1], c = params[2], d = params[3];
708  scalar_type n = params[4];
709  size_type N = gmm::mat_nrows(E);
710  GMM_ASSERT1(N == 3, "Generalized Blatz Ko hyperelastic law only defined "
711  "on dimension 3, sorry");
712  base_matrix C = E;
713  gmm::scale(C, scalar_type(2));
714  gmm::add(gmm::identity_matrix(), C);
715  compute_invariants ci(C);
716 
717  return pow(a*ci.i1() + b*sqrt(gmm::abs(ci.i3()))
718  + c*ci.i2() / ci.i3() + d, n);
719  }
720 
721  void generalized_Blatz_Ko_hyperelastic_law::sigma
722  (const base_matrix &E, base_matrix &result,
723  const base_vector &params, scalar_type det_trans) const {
724  scalar_type a = params[0], b = params[1], c = params[2], d = params[3];
725  scalar_type n = params[4];
726  size_type N = gmm::mat_nrows(E);
727  GMM_ASSERT1(N == 3, "Generalized Blatz Ko hyperelastic law only defined "
728  "on dimension 3, sorry");
729  base_matrix C = E;
730  gmm::scale(C, scalar_type(2));
731  gmm::add(gmm::identity_matrix(), C);
732  compute_invariants ci(C);
733 
734  scalar_type z = a*ci.i1() + b*sqrt(gmm::abs(ci.i3()))
735  + c*ci.i2() / ci.i3() + d;
736  scalar_type nz = n * pow(z, n-1.);
737  scalar_type di1 = nz * a;
738  scalar_type di2 = nz * c / ci.i3();
739  scalar_type di3 = nz *
740  (b / (2. * sqrt(gmm::abs(ci.i3()))) - c * ci.i2() / gmm::sqr(ci.i3()));
741 
742  gmm::copy(gmm::scaled(ci.grad_i1(), di1 * 2.0), result);
743  gmm::add(gmm::scaled(ci.grad_i2(), di2 * 2.0), result);
744  gmm::add(gmm::scaled(ci.grad_i3(), di3 * 2.0), result);
745  if (det_trans <= scalar_type(0))
746  gmm::add(gmm::scaled(C, 1e200), result);
747 
748  }
749 
750  void generalized_Blatz_Ko_hyperelastic_law::grad_sigma
751  (const base_matrix &E, base_tensor &result,
752  const base_vector &params, scalar_type) const {
753  scalar_type a = params[0], b = params[1], c = params[2], d = params[3];
754  scalar_type n = params[4];
755  size_type N = gmm::mat_nrows(E);
756  GMM_ASSERT1(N == 3, "Generalized Blatz Ko hyperelastic law only defined "
757  "on dimension 3, sorry");
758  base_matrix C = E;
759  gmm::scale(C, scalar_type(2));
760  gmm::add(gmm::identity_matrix(), C);
761  compute_invariants ci(C);
762 
763 
764  scalar_type z = a*ci.i1() + b*sqrt(gmm::abs(ci.i3()))
765  + c*ci.i2() / ci.i3() + d;
766  scalar_type nz = n * pow(z, n-1.);
767  scalar_type di1 = nz * a;
768  scalar_type di2 = nz * c / ci.i3();
769  scalar_type y = (b / (2. * sqrt(gmm::abs(ci.i3()))) - c * ci.i2() / gmm::sqr(ci.i3()));
770  scalar_type di3 = nz * y;
771 
772  gmm::copy(gmm::scaled(ci.sym_grad_grad_i1().as_vector(),
773  scalar_type(4)*di1), result.as_vector());
774  gmm::add(gmm::scaled(ci.sym_grad_grad_i2().as_vector(),
775  scalar_type(4)*di2), result.as_vector());
776  gmm::add(gmm::scaled(ci.sym_grad_grad_i3().as_vector(),
777  scalar_type(4)*di3), result.as_vector());
778 
779  scalar_type nnz = n * (n-1.) * pow(z, n-2.);
780  base_matrix A(3, 3); // second derivatives of W with respect to invariants
781  A(0, 0) = nnz * a * a;
782  A(1, 0) = A(0, 1) = nnz * a * c / ci.i3();
783  A(2, 0) = A(0, 2) = nnz * a * y;
784  A(1, 1) = nnz * c * c / gmm::sqr(ci.i3());
785  A(2, 1) = A(1, 2) = nnz * y * c / ci.i3() - nz * c / gmm::sqr(ci.i3());
786  A(2, 2) = nnz * y * y + nz * (2. * c * ci.i2() / pow(ci.i3(), 3.) - b / (4. * pow(ci.i3(), 1.5)));
787 
788  typedef const base_matrix * pointer_base_matrix__;
789  pointer_base_matrix__ di[3];
790  di[0] = &(ci.grad_i1());
791  di[1] = &(ci.grad_i2());
792  di[2] = &(ci.grad_i3());
793 
794  for (size_type j = 0; j < N; ++j)
795  for (size_type k = 0; k < N; ++k) {
796  for (size_type l1 = 0; l1 < N; ++l1)
797  for (size_type l2 = 0; l2 < N; ++l2)
798  for (size_type l3 = 0; l3 < N; ++l3)
799  for (size_type l4 = 0; l4 < N; ++l4)
800  result(l1, l2, l3, l4)
801  += 4. * A(j, k) * (*di[j])(l1, l2) * (*di[k])(l3, l4);
802  }
803 
804 // GMM_ASSERT1(check_symmetry(result) == 7,
805 // "Fourth order tensor not symmetric : " << result);
806  }
807 
808  generalized_Blatz_Ko_hyperelastic_law::generalized_Blatz_Ko_hyperelastic_law() {
809  nb_params_ = 5;
810  base_vector V(5);
811  V[0] = 1.0; V[1] = 1.0, V[2] = 1.5; V[3] = -0.5; V[4] = 1.5;
812  }
813 
814 
815  scalar_type Ciarlet_Geymonat_hyperelastic_law::strain_energy
816  (const base_matrix &E, const base_vector &params, scalar_type det_trans) const {
817  if (det_trans <= scalar_type(0)) return 1e200;
818  size_type N = gmm::mat_nrows(E);
819  scalar_type a = params[2];
820  scalar_type b = params[1]/scalar_type(2) - params[2];
821  scalar_type c = params[0]/scalar_type(4) - params[1]/scalar_type(2)
822  + params[2];
823  scalar_type d = params[0]/scalar_type(2) + params[1];
824  scalar_type e = -(scalar_type(3)*(a+b) + c);
825  base_matrix C(N, N);
826  gmm::copy(gmm::scaled(E, scalar_type(2)), C);
827  gmm::add(gmm::identity_matrix(), C);
828  scalar_type det = bgeot::lu_det(&(*(C.begin())), N);
829  return a * gmm::mat_trace(C)
830  + b * (gmm::sqr(gmm::mat_trace(C)) -
831  gmm::mat_euclidean_norm_sqr(C))/scalar_type(2)
832  + c * det - d * log(det) / scalar_type(2) + e;
833  }
834 
835  void Ciarlet_Geymonat_hyperelastic_law::sigma
836  (const base_matrix &E, base_matrix &result, const base_vector &params, scalar_type det_trans) const {
837  size_type N = gmm::mat_nrows(E);
838  scalar_type a = params[2];
839  scalar_type b = params[1]/scalar_type(2) - params[2];
840  scalar_type c = params[0]/scalar_type(4) - params[1]/scalar_type(2)
841  + params[2];
842  scalar_type d = params[0]/scalar_type(2) + params[1];
843  base_matrix C(N, N);
844  if (a > params[1]/scalar_type(2)
845  || a < params[1]/scalar_type(2) - params[0]/scalar_type(4) || a < 0)
846  GMM_WARNING1("Inconsistent third parameter for Ciarlet-Geymonat "
847  "hyperelastic law");
848  gmm::copy(gmm::scaled(E, scalar_type(2)), C);
849  gmm::add(gmm::identity_matrix(), C);
850  gmm::copy(gmm::identity_matrix(), result);
851  gmm::scale(result, scalar_type(2) * (a + b * gmm::mat_trace(C)));
852  gmm::add(gmm::scaled(C, -scalar_type(2) * b), result);
853  if (det_trans <= scalar_type(0))
854  gmm::add(gmm::scaled(C, 1e200), result);
855  else {
856  scalar_type det = bgeot::lu_inverse(&(*(C.begin())), N);
857  gmm::add(gmm::scaled(C, scalar_type(2) * c * det - d), result);
858  }
859  }
860 
861  void Ciarlet_Geymonat_hyperelastic_law::grad_sigma
862  (const base_matrix &E, base_tensor &result,const base_vector &params, scalar_type) const {
863  size_type N = gmm::mat_nrows(E);
864  // scalar_type a = params[2];
865  scalar_type b2 = params[1] - params[2]*scalar_type(2); // b*2
866  scalar_type c = params[0]/scalar_type(4) - params[1]/scalar_type(2)
867  + params[2];
868  scalar_type d = params[0]/scalar_type(2) + params[1];
869  base_matrix C(N, N);
870  gmm::copy(gmm::scaled(E, scalar_type(2)), C);
871  gmm::add(gmm::identity_matrix(), C);
872  scalar_type det = bgeot::lu_inverse(&(*(C.begin())), N);
873  std::fill(result.begin(), result.end(), scalar_type(0));
874  for (size_type i = 0; i < N; ++i)
875  for (size_type j = 0; j < N; ++j) {
876  result(i, i, j, j) += 2*b2;
877  result(i, j, i, j) -= b2;
878  result(i, j, j, i) -= b2;
879  for (size_type k = 0; k < N; ++k)
880  for (size_type l = 0; l < N; ++l)
881  result(i, j, k, l) +=
882  (C(i, k)*C(l, j) + C(i, l)*C(k, j)) * (d-scalar_type(2)*det*c)
883  + (C(i, j) * C(k, l)) * det*c*scalar_type(4);
884  }
885 
886 // GMM_ASSERT1(check_symmetry(result) == 7,
887 // "Fourth order tensor not symmetric : " << result);
888  }
889 
890 
891  int levi_civita(int i, int j, int k) {
892  int ii=i+1;
893  int jj=j+1;
894  int kk=k+1; //i,j,k from 0 to 2 !
895  return static_cast<int>
896  (int(- 1)*(static_cast<int>(pow(double(ii-jj),2.))%3)
897  * (static_cast<int> (pow(double(ii-kk),double(2)))%3 )
898  * (static_cast<int> (pow(double(jj-kk),double(2)))%3)
899  * (pow(double(jj-(ii%3))-double(0.5),double(2))-double(1.25)));
900  }
901 
902 
903 
904  scalar_type plane_strain_hyperelastic_law::strain_energy
905  (const base_matrix &E, const base_vector &params, scalar_type det_trans) const {
906  GMM_ASSERT1(gmm::mat_nrows(E) == 2, "Plane strain law is for 2D only.");
907  base_matrix E3D(3,3);
908  E3D(0,0)=E(0,0); E3D(1,0)=E(1,0); E3D(0,1)=E(0,1); E3D(1,1)=E(1,1);
909  return pl->strain_energy(E3D, params, det_trans);
910  }
911 
912  void plane_strain_hyperelastic_law::sigma
913  (const base_matrix &E, base_matrix &result,const base_vector &params, scalar_type det_trans) const {
914  GMM_ASSERT1(gmm::mat_nrows(E) == 2, "Plane strain law is for 2D only.");
915  base_matrix E3D(3,3), result3D(3,3);
916  E3D(0,0)=E(0,0); E3D(1,0)=E(1,0); E3D(0,1)=E(0,1); E3D(1,1)=E(1,1);
917  pl->sigma(E3D, result3D, params, det_trans);
918  result(0,0) = result3D(0,0); result(1,0) = result3D(1,0);
919  result(0,1) = result3D(0,1); result(1,1) = result3D(1,1);
920  }
921 
922  void plane_strain_hyperelastic_law::grad_sigma
923  (const base_matrix &E, base_tensor &result,const base_vector &params, scalar_type det_trans) const {
924  GMM_ASSERT1(gmm::mat_nrows(E) == 2, "Plane strain law is for 2D only.");
925  base_matrix E3D(3,3);
926  base_tensor result3D(3,3,3,3);
927  E3D(0,0)=E(0,0); E3D(1,0)=E(1,0); E3D(0,1)=E(0,1); E3D(1,1)=E(1,1);
928  pl->grad_sigma(E3D, result3D, params, det_trans);
929  result(0,0,0,0) = result3D(0,0,0,0); result(1,0,0,0) = result3D(1,0,0,0);
930  result(0,1,0,0) = result3D(0,1,0,0); result(1,1,0,0) = result3D(1,1,0,0);
931  result(0,0,1,0) = result3D(0,0,1,0); result(1,0,1,0) = result3D(1,0,1,0);
932  result(0,1,1,0) = result3D(0,1,1,0); result(1,1,1,0) = result3D(1,1,1,0);
933  result(0,0,0,1) = result3D(0,0,0,1); result(1,0,0,1) = result3D(1,0,0,1);
934  result(0,1,0,1) = result3D(0,1,0,1); result(1,1,0,1) = result3D(1,1,0,1);
935  result(0,0,1,1) = result3D(0,0,1,1); result(1,0,1,1) = result3D(1,0,1,1);
936  result(0,1,1,1) = result3D(0,1,1,1); result(1,1,1,1) = result3D(1,1,1,1);
937  }
938 
939 
940 
941 
942 
943 
944 
945  //=========================================================================
946  //
947  // Nonlinear elasticity (old) Brick
948  //
949  //=========================================================================
950 
951  struct nonlinear_elasticity_brick : public virtual_brick {
952 
953  phyperelastic_law AHL;
954 
955  virtual void asm_real_tangent_terms(const model &md, size_type /* ib */,
956  const model::varnamelist &vl,
957  const model::varnamelist &dl,
958  const model::mimlist &mims,
959  model::real_matlist &matl,
960  model::real_veclist &vecl,
961  model::real_veclist &,
962  size_type region,
963  build_version version) const {
964  GMM_ASSERT1(mims.size() == 1,
965  "Nonlinear elasticity brick need a single mesh_im");
966  GMM_ASSERT1(vl.size() == 1,
967  "Nonlinear elasticity brick need a single variable");
968  GMM_ASSERT1(dl.size() == 1,
969  "Wrong number of data for nonlinear elasticity brick, "
970  << dl.size() << " should be 1 (vector).");
971  GMM_ASSERT1(matl.size() == 1, "Wrong number of terms for nonlinear "
972  "elasticity brick");
973 
974  const model_real_plain_vector &u = md.real_variable(vl[0]);
975  const mesh_fem &mf_u = *(md.pmesh_fem_of_variable(vl[0]));
976 
977  const mesh_fem *mf_params = md.pmesh_fem_of_variable(dl[0]);
978  const model_real_plain_vector &params = md.real_variable(dl[0]);
979  const mesh_im &mim = *mims[0];
980 
981  size_type sl = gmm::vect_size(params);
982  if (mf_params) sl = sl * mf_params->get_qdim() / mf_params->nb_dof();
983  GMM_ASSERT1(sl == AHL->nb_params(), "Wrong number of coefficients for the "
984  "nonlinear constitutive elastic law");
985 
986  mesh_region rg(region);
987  mf_u.linked_mesh().intersect_with_mpi_region(rg);
988 
989  if (version & model::BUILD_MATRIX) {
990  gmm::clear(matl[0]);
991  GMM_TRACE2("Nonlinear elasticity stiffness matrix assembly");
992  asm_nonlinear_elasticity_tangent_matrix
993  (matl[0], mim, mf_u, u, mf_params, params, *AHL, rg);
994  }
995 
996 
997  if (version & model::BUILD_RHS) {
998  asm_nonlinear_elasticity_rhs(vecl[0], mim,
999  mf_u, u, mf_params, params, *AHL, rg);
1000  gmm::scale(vecl[0], scalar_type(-1));
1001  }
1002 
1003  }
1004 
1005  nonlinear_elasticity_brick(const phyperelastic_law &AHL_)
1006  : AHL(AHL_) {
1007  set_flags("Nonlinear elasticity brick", false /* is linear*/,
1008  true /* is symmetric */, true /* is coercive */,
1009  true /* is real */, false /* is complex */);
1010  }
1011 
1012  };
1013 
1014  //=========================================================================
1015  // Add a nonlinear elasticity brick.
1016  //=========================================================================
1017 
1018  // Deprecated brick
1019  size_type add_nonlinear_elasticity_brick
1020  (model &md, const mesh_im &mim, const std::string &varname,
1021  const phyperelastic_law &AHL, const std::string &dataname,
1022  size_type region) {
1023  pbrick pbr = std::make_shared<nonlinear_elasticity_brick>(AHL);
1024 
1025  model::termlist tl;
1026  tl.push_back(model::term_description(varname, varname, true));
1027  model::varnamelist dl(1, dataname);
1028  model::varnamelist vl(1, varname);
1029  return md.add_brick(pbr, vl, dl, tl, model::mimlist(1,&mim), region);
1030  }
1031 
1032  //=========================================================================
1033  // Von Mises or Tresca stress computation.
1034  //=========================================================================
1035 
1036  void compute_Von_Mises_or_Tresca(model &md,
1037  const std::string &varname,
1038  const phyperelastic_law &AHL,
1039  const std::string &dataname,
1040  const mesh_fem &mf_vm,
1041  model_real_plain_vector &VM,
1042  bool tresca) {
1043  GMM_ASSERT1(gmm::vect_size(VM) == mf_vm.nb_dof(),
1044  "The vector has not the good size");
1045  const mesh_fem &mf_u = md.mesh_fem_of_variable(varname);
1046  const model_real_plain_vector &u = md.real_variable(varname);
1047  const mesh_fem *mf_params = md.pmesh_fem_of_variable(dataname);
1048  const model_real_plain_vector &params = md.real_variable(dataname);
1049 
1050  size_type sl = gmm::vect_size(params);
1051  if (mf_params) sl = sl * mf_params->get_qdim() / mf_params->nb_dof();
1052  GMM_ASSERT1(sl == AHL->nb_params(), "Wrong number of coefficients for "
1053  "the nonlinear constitutive elastic law");
1054 
1055  unsigned N = unsigned(mf_u.linked_mesh().dim());
1056  unsigned NP = unsigned(AHL->nb_params()), NFem = mf_u.get_qdim();
1057  model_real_plain_vector GRAD(mf_vm.nb_dof()*NFem*N);
1058  model_real_plain_vector PARAMS(mf_vm.nb_dof()*NP);
1059  if (mf_params) interpolation(*mf_params, mf_vm, params, PARAMS);
1060  compute_gradient(mf_u, mf_vm, u, GRAD);
1061  base_matrix E(N, N), gradphi(NFem,N),gradphit(N,NFem), Id(N, N),
1062  sigmahathat(N,N),aux(NFem,N), sigma(NFem,NFem),
1063  IdNFem(NFem, NFem);
1064  base_vector p(NP);
1065  if (!mf_params) gmm::copy(params, p);
1066  base_vector eig(NFem);
1067  base_vector ez(NFem); // vector normal at deformed surface, (ex X ey)
1068  double normEz(0); //norm of ez
1069  gmm::copy(gmm::identity_matrix(), Id);
1070  gmm::copy(gmm::identity_matrix(), IdNFem);
1071  for (size_type i = 0; i < mf_vm.nb_dof(); ++i) {
1072  gmm::resize(gradphi,NFem,N);
1073  std::copy(GRAD.begin()+i*NFem*N, GRAD.begin()+(i+1)*NFem*N,
1074  gradphit.begin());
1075  gmm::copy(gmm::transposed(gradphit),gradphi);
1076  for (unsigned int alpha = 0; alpha <N; ++alpha)
1077  gradphi(alpha, alpha)+=1;
1078  gmm::mult(gmm::transposed(gradphi), gradphi, E);
1079  gmm::add(gmm::scaled(Id, -scalar_type(1)), E);
1080  gmm::scale(E, scalar_type(1)/scalar_type(2));
1081  if (mf_params)
1082  gmm::copy(gmm::sub_vector(PARAMS, gmm::sub_interval(i*NP,NP)), p);
1083  AHL->sigma(E, sigmahathat, p, scalar_type(1));
1084  if (NFem == 3 && N == 2) {
1085  //jyh : compute ez, normal on deformed surface
1086  for (unsigned int l = 0; l <NFem; ++l) {
1087  ez[l]=0;
1088  for (unsigned int m = 0; m <NFem; ++m)
1089  for (unsigned int n = 0; n <NFem; ++n){
1090  ez[l]+=levi_civita(l,m,n)*gradphi(m,0)*gradphi(n,1);
1091  }
1092  normEz= gmm::vect_norm2(ez);
1093  }
1094  //jyh : end compute ez
1095  }
1096  gmm::mult(gradphi, sigmahathat, aux);
1097  gmm::mult(aux, gmm::transposed(gradphi), sigma);
1098 
1099  /* jyh : complete gradphi for virtual 3rd dim (perpendicular to
1100  deformed surface, same thickness) */
1101  if (NFem == 3 && N == 2) {
1102  gmm::resize(gradphi,NFem,NFem);
1103  for (unsigned int ll = 0; ll <NFem; ++ll)
1104  for (unsigned int ii = 0; ii <NFem; ++ii)
1105  for (unsigned int jj = 0; jj <NFem; ++jj)
1106  gradphi(ll,2)+=(levi_civita(ll,ii,jj)*gradphi(ii,0)
1107  *gradphi(jj,1))/normEz;
1108  //jyh : end complete graphi
1109  }
1110 
1111  gmm::scale(sigma, scalar_type(1) / bgeot::lu_det(&(*(gradphi.begin())), NFem));
1112 
1113  if (!tresca) {
1114  /* von mises: norm(deviator(sigma)) */
1115  gmm::add(gmm::scaled(IdNFem, -gmm::mat_trace(sigma) / NFem), sigma);
1116 
1117  //jyh : von mises stress=sqrt(3/2)* norm(sigma) ?
1118  VM[i] = sqrt(3.0/2)*gmm::mat_euclidean_norm(sigma);
1119  } else {
1120  /* else compute the tresca criterion */
1121  //jyh : to be adapted for membrane if necessary
1122  gmm::symmetric_qr_algorithm(sigma, eig);
1123  std::sort(eig.begin(), eig.end());
1124  VM[i] = eig.back() - eig.front();
1125  }
1126  }
1127  }
1128 
1129 
1130  void compute_sigmahathat(model &md,
1131  const std::string &varname,
1132  const phyperelastic_law &AHL,
1133  const std::string &dataname,
1134  const mesh_fem &mf_sigma,
1135  model_real_plain_vector &SIGMA) {
1136  const mesh_fem &mf_u = md.mesh_fem_of_variable(varname);
1137  const model_real_plain_vector &u = md.real_variable(varname);
1138  const mesh_fem *mf_params = md.pmesh_fem_of_variable(dataname);
1139  const model_real_plain_vector &params = md.real_variable(dataname);
1140 
1141  size_type sl = gmm::vect_size(params);
1142  if (mf_params) sl = sl * mf_params->get_qdim() / mf_params->nb_dof();
1143  GMM_ASSERT1(sl == AHL->nb_params(), "Wrong number of coefficients for "
1144  "the nonlinear constitutive elastic law");
1145 
1146  unsigned N = unsigned(mf_u.linked_mesh().dim());
1147  unsigned NP = unsigned(AHL->nb_params()), NFem = mf_u.get_qdim();
1148  GMM_ASSERT1(mf_sigma.nb_dof() > 0, "Bad mf_sigma");
1149  size_type qqdim = mf_sigma.get_qdim();
1150  size_type ratio = N*N / qqdim;
1151 
1152  GMM_ASSERT1(((ratio == 1) || (ratio == N*N)) &&
1153  (gmm::vect_size(SIGMA) == mf_sigma.nb_dof()*ratio),
1154  "The vector has not the good size");
1155 
1156  model_real_plain_vector GRAD(mf_sigma.nb_dof()*ratio*NFem/N);
1157  model_real_plain_vector PARAMS(mf_sigma.nb_dof()*NP);
1158 
1159 
1160  getfem::mesh_trans_inv mti(mf_sigma.linked_mesh());
1161  if (mf_params) {
1162  for (size_type i = 0; i < mf_sigma.nb_dof(); ++i)
1163  mti.add_point(mf_sigma.point_of_basic_dof(i));
1164  interpolation(*mf_params, mti, params, PARAMS);
1165  }
1166 
1167  compute_gradient(mf_u, mf_sigma, u, GRAD);
1168  base_matrix E(N, N), gradphi(NFem,N),gradphit(N,NFem), Id(N, N),
1169  sigmahathat(N,N),aux(NFem,N), sigma(NFem,NFem),
1170  IdNFem(NFem, NFem);
1171 
1172 
1173  base_vector p(NP);
1174  if (!mf_params) gmm::copy(params, p);
1175  base_vector eig(NFem);
1176  base_vector ez(NFem); // vector normal at deformed surface, (ex X ey)
1177  gmm::copy(gmm::identity_matrix(), Id);
1178  gmm::copy(gmm::identity_matrix(), IdNFem);
1179 
1180  // cout << "GRAD = " << GRAD << endl;
1181  // GMM_ASSERT1(false, "oops");
1182 
1183  for (size_type i = 0; i < mf_sigma.nb_dof()/qqdim; ++i) {
1184 // cout << "GRAD size = " << gmm::vect_size(GRAD) << endl;
1185 // cout << "i = " << i << endl;
1186 // cout << "i*NFem*N = " << i*NFem*N << endl;
1187 
1188  gmm::resize(gradphi,NFem,N);
1189  std::copy(GRAD.begin()+i*NFem*N, GRAD.begin()+(i+1)*NFem*N,
1190  gradphit.begin());
1191  // cout << "GRAD = " << gradphit << endl;
1192  gmm::copy(gmm::transposed(gradphit),gradphi);
1193  for (unsigned int alpha = 0; alpha <N; ++alpha)
1194  gradphi(alpha, alpha) += scalar_type(1);
1195  gmm::mult(gmm::transposed(gradphi), gradphi, E);
1196  gmm::add(gmm::scaled(Id, -scalar_type(1)), E);
1197  gmm::scale(E, scalar_type(1)/scalar_type(2));
1198  if (mf_params)
1199  gmm::copy(gmm::sub_vector(PARAMS, gmm::sub_interval(i*ratio*NP,NP)),p);
1200  // cout << "E = " << E << endl;
1201  AHL->sigma(E, sigmahathat, p, scalar_type(1));
1202  // cout << "ok" << endl;
1203  std::copy(sigmahathat.begin(), sigmahathat.end(), SIGMA.begin()+i*N*N);
1204  }
1205  }
1206 
1207 
1208 
1209  // ----------------------------------------------------------------------
1210  //
1211  // Nonlinear incompressibility brick
1212  //
1213  // ----------------------------------------------------------------------
1214 
1215  struct nonlinear_incompressibility_brick : public virtual_brick {
1216 
1217  virtual void asm_real_tangent_terms(const model &md, size_type,
1218  const model::varnamelist &vl,
1219  const model::varnamelist &dl,
1220  const model::mimlist &mims,
1221  model::real_matlist &matl,
1222  model::real_veclist &vecl,
1223  model::real_veclist &veclsym,
1224  size_type region,
1225  build_version version) const {
1226 
1227  GMM_ASSERT1(matl.size() == 2, "Wrong number of terms for nonlinear "
1228  "incompressibility brick");
1229  GMM_ASSERT1(dl.size() == 0, "Nonlinear incompressibility brick need no "
1230  "data");
1231  GMM_ASSERT1(mims.size() == 1, "Nonlinear incompressibility brick need a "
1232  "single mesh_im");
1233  GMM_ASSERT1(vl.size() == 2, "Wrong number of variables for nonlinear "
1234  "incompressibility brick");
1235 
1236  const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
1237  const mesh_fem &mf_p = md.mesh_fem_of_variable(vl[1]);
1238  const model_real_plain_vector &u = md.real_variable(vl[0]);
1239  const model_real_plain_vector &p = md.real_variable(vl[1]);
1240  const mesh_im &mim = *mims[0];
1241  mesh_region rg(region);
1242  mim.linked_mesh().intersect_with_mpi_region(rg);
1243 
1244  if (version & model::BUILD_MATRIX) {
1245  gmm::clear(matl[0]);
1246  gmm::clear(matl[1]);
1247  asm_nonlinear_incomp_tangent_matrix(matl[0], matl[1],
1248  mim, mf_u, mf_p, u, p, rg);
1249  }
1250 
1251  if (version & model::BUILD_RHS) {
1252  asm_nonlinear_incomp_rhs(vecl[0], veclsym[1], mim, mf_u, mf_p,u,p, rg);
1253  gmm::scale(vecl[0], scalar_type(-1));
1254  gmm::scale(veclsym[1], scalar_type(-1));
1255  }
1256  }
1257 
1258 
1259  nonlinear_incompressibility_brick() {
1260  set_flags("Nonlinear incompressibility brick",
1261  false /* is linear*/,
1262  true /* is symmetric */, false /* is coercive */,
1263  true /* is real */, false /* is complex */);
1264  }
1265 
1266 
1267  };
1268 
1269  size_type add_nonlinear_incompressibility_brick
1270  (model &md, const mesh_im &mim, const std::string &varname,
1271  const std::string &multname, size_type region) {
1272  pbrick pbr = std::make_shared<nonlinear_incompressibility_brick>();
1273  model::termlist tl;
1274  tl.push_back(model::term_description(varname, varname, true));
1275  tl.push_back(model::term_description(varname, multname, true));
1276  model::varnamelist vl(1, varname);
1277  vl.push_back(multname);
1278  model::varnamelist dl;
1279  return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
1280  }
1281 
1282 
1283 
1284 
1285 
1286  // ----------------------------------------------------------------------
1287  //
1288  // High-level Generic assembly operators
1289  //
1290  // ----------------------------------------------------------------------
1291 
1292 
1293  static void ga_init_scalar_(bgeot::multi_index &mi) { mi.resize(0); }
1294  static void ga_init_square_matrix_(bgeot::multi_index &mi, size_type N)
1295  { mi.resize(2); mi[0] = mi[1] = N; }
1296 
1297 
1298  // Matrix_i2 Operator (second invariant of square matrix of size >= 2)
1299  // (For 2x2 matrices, it is equivalent to det(M))
1300  struct matrix_i2_operator : public ga_nonlinear_operator {
1301  bool result_size(const arg_list &args, bgeot::multi_index &sizes) const {
1302  if (args.size() != 1 || args[0]->sizes().size() != 2
1303  || args[0]->sizes()[0] != args[0]->sizes()[1]) return false;
1304  ga_init_scalar_(sizes);
1305  return true;
1306  }
1307 
1308  // Value : (Trace(M))^2 - Trace(M^2))/2
1309  void value(const arg_list &args, base_tensor &result) const {
1310  size_type N = args[0]->sizes()[0];
1311  const base_tensor &t = *args[0];
1312  scalar_type tr = scalar_type(0);
1313  for (size_type i = 0; i < N; ++i) tr += t[i*N+i];
1314  scalar_type tr2 = scalar_type(0);
1315  for (size_type i = 0; i < N; ++i)
1316  for (size_type j = 0; j < N; ++j)
1317  tr2 += t[i+ j*N] * t[j + i*N];
1318  result[0] = (tr*tr-tr2)/scalar_type(2);
1319  }
1320 
1321  // Derivative : Trace(M)I - M^T
1322  void derivative(const arg_list &args, size_type,
1323  base_tensor &result) const {
1324  size_type N = args[0]->sizes()[0];
1325  const base_tensor &t = *args[0];
1326  scalar_type tr = scalar_type(0);
1327  for (size_type i = 0; i < N; ++i) tr += t[i*N+i];
1328  base_tensor::iterator it = result.begin();
1329  for (size_type j = 0; j < N; ++j)
1330  for (size_type i = 0; i < N; ++i, ++it)
1331  *it = ((i == j) ? tr : scalar_type(0)) - t[i*N+j];
1332  GMM_ASSERT1(it == result.end(), "Internal error");
1333  }
1334 
1335  // Second derivative : I@I - \delta_{il}\delta_{jk}
1336  void second_derivative(const arg_list &args, size_type, size_type,
1337  base_tensor &result) const {
1338  size_type N = args[0]->sizes()[0];
1339  gmm::clear(result.as_vector());
1340  for (size_type i = 0; i < N; ++i)
1341  for (size_type j = 0; j < N; ++j) {
1342  result[(N+1)*(i+N*N*j)] += scalar_type(1);
1343  result[(N+1)*N*j + i*(N*N*N + 1)] -= scalar_type(1);
1344  }
1345  }
1346  };
1347 
1348 
1349  // Matrix_j1 Operator
1350  struct matrix_j1_operator : public ga_nonlinear_operator {
1351  bool result_size(const arg_list &args, bgeot::multi_index &sizes) const {
1352  if (args.size() != 1 || args[0]->sizes().size() != 2
1353  || args[0]->sizes()[0] != args[0]->sizes()[1]) return false;
1354  ga_init_scalar_(sizes);
1355  return true;
1356  }
1357 
1358  // Value : Trace(M)/(det(M)^1/3)
1359  void value(const arg_list &args, base_tensor &result) const {
1360  size_type N = args[0]->sizes()[0];
1361  base_matrix M(N, N);
1362  gmm::copy(args[0]->as_vector(), M.as_vector());
1363  scalar_type det = bgeot::lu_det(&(*(M.begin())), N);
1364  scalar_type tr = scalar_type(0);
1365  for (size_type i = 0; i < N; ++i) tr += M(i,i);
1366  if (det > 0)
1367  result[0] = tr / pow(det, scalar_type(1)/scalar_type(3));
1368  else
1369  result[0] = 1.E200;
1370  }
1371 
1372  // Derivative : (I-Trace(M)*M^{-T}/3)/(det(M)^1/3)
1373  void derivative(const arg_list &args, size_type,
1374  base_tensor &result) const {
1375  size_type N = args[0]->sizes()[0];
1376  base_matrix M(N, N);
1377  gmm::copy(args[0]->as_vector(), M.as_vector());
1378  scalar_type tr = scalar_type(0);
1379  for (size_type i = 0; i < N; ++i) tr += M(i,i);
1380  scalar_type det = bgeot::lu_inverse(&(*(M.begin())), N);
1381  if (det > 0) {
1382  base_tensor::iterator it = result.begin();
1383  for (size_type j = 0; j < N; ++j)
1384  for (size_type i = 0; i < N; ++i, ++it)
1385  *it = (((i == j) ? scalar_type(1) : scalar_type(0))
1386  - tr*M(j,i)/scalar_type(3))
1387  / pow(det, scalar_type(1)/scalar_type(3));
1388  GMM_ASSERT1(it == result.end(), "Internal error");
1389  } else
1390  std::fill(result.begin(), result.end(), 1.E200);
1391  }
1392 
1393  // Second derivative : (-M^{-T}@I + Trace(M)*M^{-T}_{ik}M^{-T}_{lj}
1394  // -I@M^{-T} + Trace(M)*M^{-T}@M^{-T}/3)/(3det(M)^1/3)
1395  void second_derivative(const arg_list &args, size_type, size_type,
1396  base_tensor &result) const {
1397  size_type N = args[0]->sizes()[0];
1398  base_matrix M(N, N);
1399  gmm::copy(args[0]->as_vector(), M.as_vector());
1400  scalar_type tr = scalar_type(0);
1401  for (size_type i = 0; i < N; ++i) tr += M(i,i);
1402  scalar_type det = bgeot::lu_inverse(&(*(M.begin())), N);
1403  if (det > 0) {
1404  base_tensor::iterator it = result.begin();
1405  for (size_type l = 0; l < N; ++l)
1406  for (size_type k = 0; k < N; ++k)
1407  for (size_type j = 0; j < N; ++j)
1408  for (size_type i = 0; i < N; ++i, ++it)
1409  *it = (- ((k == l) ? M(j, i) : scalar_type(0))
1410  + tr*M(i,k)*M(l,j)
1411  - ((i == j) ? M(l, k) : scalar_type(0))
1412  + tr*M(j,i)*M(k,l)/ scalar_type(3))
1413  / (scalar_type(3)*pow(det, scalar_type(1)/scalar_type(3)));
1414  GMM_ASSERT1(it == result.end(), "Internal error");
1415  } else
1416  std::fill(result.begin(), result.end(), 1.E200);
1417  }
1418  };
1419 
1420 
1421  // Matrix_j2 Operator
1422  struct matrix_j2_operator : public ga_nonlinear_operator {
1423  bool result_size(const arg_list &args, bgeot::multi_index &sizes) const {
1424  if (args.size() != 1 || args[0]->sizes().size() != 2
1425  || args[0]->sizes()[0] != args[0]->sizes()[1]) return false;
1426  ga_init_scalar_(sizes);
1427  return true;
1428  }
1429 
1430  // Value : i2(M)/(det(M)^2/3)
1431  void value(const arg_list &args, base_tensor &result) const {
1432  size_type N = args[0]->sizes()[0];
1433  base_matrix M(N, N);
1434  gmm::copy(args[0]->as_vector(), M.as_vector());
1435  scalar_type tr = scalar_type(0);
1436  for (size_type i = 0; i < N; ++i) tr += M(i,i);
1437  scalar_type tr2 = scalar_type(0);
1438  for (size_type i = 0; i < N; ++i)
1439  for (size_type j = 0; j < N; ++j)
1440  tr2 += M(i,j)*M(j,i);
1441  scalar_type i2 = (tr*tr-tr2)/scalar_type(2);
1442  scalar_type det = bgeot::lu_det(&(*(M.begin())), N);
1443  if (det > 0)
1444  result[0] = i2 / pow(det, scalar_type(2)/scalar_type(3));
1445  else
1446  result[0] = 1.E200;
1447  }
1448 
1449  // Derivative : (Trace(M)*I-M^T-2i2(M)M^{-T}/3)/(det(M)^2/3)
1450  void derivative(const arg_list &args, size_type,
1451  base_tensor &result) const {
1452  size_type N = args[0]->sizes()[0];
1453  const base_tensor &t = *args[0];
1454  base_matrix M(N, N);
1455  gmm::copy(t.as_vector(), M.as_vector());
1456  scalar_type tr = scalar_type(0);
1457  for (size_type i = 0; i < N; ++i) tr += M(i,i);
1458  scalar_type tr2 = scalar_type(0);
1459  for (size_type i = 0; i < N; ++i)
1460  for (size_type j = 0; j < N; ++j)
1461  tr2 += M(i,j)*M(j,i);
1462  scalar_type i2 = (tr*tr-tr2)/scalar_type(2);
1463  scalar_type det = bgeot::lu_inverse(&(*(M.begin())), N);
1464  base_tensor::iterator it = result.begin();
1465  for (size_type j = 0; j < N; ++j)
1466  for (size_type i = 0; i < N; ++i, ++it)
1467  *it = (((i == j) ? tr : scalar_type(0)) - t[j+N*i]
1468  - scalar_type(2)*i2*M(j,i)/(scalar_type(3)))
1469  / pow(det, scalar_type(2)/scalar_type(3));
1470  GMM_ASSERT1(it == result.end(), "Internal error");
1471  }
1472 
1473  // Second derivative
1474  void second_derivative(const arg_list &args, size_type, size_type,
1475  base_tensor &result) const {
1476  size_type N = args[0]->sizes()[0];
1477  const base_tensor &t = *args[0];
1478  base_matrix M(N, N);
1479  gmm::copy(t.as_vector(), M.as_vector());
1480  scalar_type tr = scalar_type(0);
1481  for (size_type i = 0; i < N; ++i) tr += M(i,i);
1482  scalar_type tr2 = scalar_type(0);
1483  for (size_type i = 0; i < N; ++i)
1484  for (size_type j = 0; j < N; ++j)
1485  tr2 += M(i,j)*M(j,i);
1486  scalar_type i2 = (tr*tr-tr2)/scalar_type(2);
1487  scalar_type det = bgeot::lu_inverse(&(*(M.begin())), N);
1488  base_tensor::iterator it = result.begin();
1489  for (size_type l = 0; l < N; ++l)
1490  for (size_type k = 0; k < N; ++k)
1491  for (size_type j = 0; j < N; ++j)
1492  for (size_type i = 0; i < N; ++i, ++it)
1493  *it = ( + (((i==j) && (k==l)) ? 1. : 0.)
1494  - (((i==l) && (k==j)) ? 1. : 0.)
1495  + 10.*i2*M(j,i)*M(l,k)/(9.)
1496  - 2.*(M(j,i)*(tr*((k==l) ? 1.:0.)-t[l+N*k]))/(3.)
1497  - 2.*(M(l,k)*(tr*((i==j) ? 1.:0.)-t[j+N*i]))/(3.)
1498  - 2.*i2*(M(j,i)*M(l,k)-M(j,k)*M(l,i))/(3.))
1499  / pow(det, scalar_type(2)/scalar_type(3));
1500  }
1501  };
1502 
1503  // Right-Cauchy-Green operator (F^{T}F)
1504  struct Right_Cauchy_Green_operator : public ga_nonlinear_operator {
1505  bool result_size(const arg_list &args, bgeot::multi_index &sizes) const {
1506  if (args.size() != 1 || args[0]->sizes().size() != 2) return false;
1507  ga_init_square_matrix_(sizes, args[0]->sizes()[1]);
1508  return true;
1509  }
1510 
1511  // Value : F^{T}F
1512  void value(const arg_list &args, base_tensor &result) const {
1513  // to be verified
1514  size_type m = args[0]->sizes()[0], n = args[0]->sizes()[1];
1515  base_tensor::iterator it = result.begin();
1516  for (size_type j = 0; j < n; ++j)
1517  for (size_type i = 0; i < n; ++i, ++it) {
1518  *it = scalar_type(0);
1519  for (size_type k = 0; k < m; ++k)
1520  *it += (*(args[0]))[i*m+k] * (*(args[0]))[j*m+k];
1521  }
1522  }
1523 
1524  // Derivative : F{kj}delta{li}+F{ki}delta{lj}
1525  // (comes from H -> H^{T}F + F^{T}H)
1526  void derivative(const arg_list &args, size_type,
1527  base_tensor &result) const { // to be verified
1528  size_type m = args[0]->sizes()[0], n = args[0]->sizes()[1];
1529  base_tensor::iterator it = result.begin();
1530  for (size_type l = 0; l < n; ++l)
1531  for (size_type k = 0; k < m; ++k)
1532  for (size_type j = 0; j < n; ++j)
1533  for (size_type i = 0; i < n; ++i, ++it) {
1534  *it = scalar_type(0);
1535  if (l == i) *it += (*(args[0]))[j*m+k];
1536  if (l == j) *it += (*(args[0]))[i*m+k];
1537  }
1538  GMM_ASSERT1(it == result.end(), "Internal error");
1539  }
1540 
1541  // Second derivative :
1542  // A{ijklop}=delta{ok}delta{li}delta{pj} + delta{ok}delta{pi}delta{lj}
1543  // comes from (H,K) -> H^{T}K + K^{T}H
1544  void second_derivative(const arg_list &args, size_type, size_type,
1545  base_tensor &result) const { // to be verified
1546  size_type m = args[0]->sizes()[0], n = args[0]->sizes()[1];
1547  base_tensor::iterator it = result.begin();
1548  for (size_type p = 0; p < n; ++p)
1549  for (size_type o = 0; o < m; ++o)
1550  for (size_type l = 0; l < n; ++l)
1551  for (size_type k = 0; k < m; ++k)
1552  for (size_type j = 0; j < n; ++j)
1553  for (size_type i = 0; i < n; ++i, ++it) {
1554  *it = scalar_type(0);
1555  if (o == k) {
1556  if (l == i && p == j) *it += scalar_type(1);
1557  if (p == i && l == j) *it += scalar_type(1);
1558  }
1559  }
1560  GMM_ASSERT1(it == result.end(), "Internal error");
1561  }
1562  };
1563 
1564  // Left-Cauchy-Green operator (FF^{T})
1565  struct Left_Cauchy_Green_operator : public ga_nonlinear_operator {
1566  bool result_size(const arg_list &args, bgeot::multi_index &sizes) const {
1567  if (args.size() != 1 || args[0]->sizes().size() != 2) return false;
1568  ga_init_square_matrix_(sizes, args[0]->sizes()[0]);
1569  return true;
1570  }
1571 
1572  // Value : FF^{T}
1573  void value(const arg_list &args, base_tensor &result) const {
1574  // to be verified
1575  size_type m = args[0]->sizes()[0], n = args[0]->sizes()[1];
1576  base_tensor::iterator it = result.begin();
1577  for (size_type j = 0; j < m; ++j)
1578  for (size_type i = 0; i < m; ++i, ++it) {
1579  *it = scalar_type(0);
1580  for (size_type k = 0; k < n; ++k)
1581  *it += (*(args[0]))[k*m+i] * (*(args[0]))[k*m+j];
1582  }
1583  }
1584 
1585  // Derivative : F{jl}delta{ik}+F{il}delta{kj}
1586  // (comes from H -> HF^{T} + FH^{T})
1587  void derivative(const arg_list &args, size_type,
1588  base_tensor &result) const {
1589  size_type m = args[0]->sizes()[0], n = args[0]->sizes()[1];
1590  base_tensor::iterator it = result.begin();
1591  for (size_type l = 0; l < n; ++l)
1592  for (size_type k = 0; k < m; ++k)
1593  for (size_type j = 0; j < m; ++j)
1594  for (size_type i = 0; i < m; ++i, ++it) {
1595  *it = scalar_type(0);
1596  if (k == i) *it += (*(args[0]))[l*m+j];
1597  if (k == j) *it += (*(args[0]))[l*m+i];
1598  }
1599  GMM_ASSERT1(it == result.end(), "Internal error");
1600  }
1601 
1602  // Second derivative :
1603  // A{ijklop}=delta{ki}delta{lp}delta{oj} + delta{oi}delta{pl}delta{kj}
1604  // comes from (H,K) -> HK^{T} + KH^{T}
1605  void second_derivative(const arg_list &args, size_type, size_type,
1606  base_tensor &result) const {
1607  size_type m = args[0]->sizes()[0], n = args[0]->sizes()[1];
1608  base_tensor::iterator it = result.begin();
1609  for (size_type p = 0; p < n; ++p)
1610  for (size_type o = 0; o < m; ++o)
1611  for (size_type l = 0; l < n; ++l)
1612  for (size_type k = 0; k < m; ++k)
1613  for (size_type j = 0; j < m; ++j)
1614  for (size_type i = 0; i < m; ++i, ++it) {
1615  *it = scalar_type(0);
1616  if (p == l) {
1617  if (k == i && o == j) *it += scalar_type(1);
1618  if (o == i && k == j) *it += scalar_type(1);
1619  }
1620  }
1621  GMM_ASSERT1(it == result.end(), "Internal error");
1622  }
1623  };
1624 
1625 
1626  // Green-Lagrangian operator (F^{T}F-I)/2
1627  struct Green_Lagrangian_operator : public ga_nonlinear_operator {
1628  bool result_size(const arg_list &args, bgeot::multi_index &sizes) const {
1629  if (args.size() != 1 || args[0]->sizes().size() != 2) return false;
1630  ga_init_square_matrix_(sizes, args[0]->sizes()[1]);
1631  return true;
1632  }
1633 
1634  // Value : F^{T}F
1635  void value(const arg_list &args, base_tensor &result) const {
1636  // to be verified
1637  size_type m = args[0]->sizes()[0], n = args[0]->sizes()[1];
1638  base_tensor::iterator it = result.begin();
1639  for (size_type j = 0; j < n; ++j)
1640  for (size_type i = 0; i < n; ++i, ++it) {
1641  *it = scalar_type(0);
1642  for (size_type k = 0; k < m; ++k)
1643  *it += (*(args[0]))[i*m+k]*(*(args[0]))[j*m+k]*scalar_type(0.5);
1644  if (i == j) *it -= scalar_type(0.5);
1645  }
1646  }
1647 
1648  // Derivative : (F{kj}delta{li}+F{ki}delta{lj})/2
1649  // (comes from H -> (H^{T}F + F^{T}H)/2)
1650  void derivative(const arg_list &args, size_type,
1651  base_tensor &result) const { // to be verified
1652  size_type m = args[0]->sizes()[0], n = args[0]->sizes()[1];
1653  base_tensor::iterator it = result.begin();
1654  for (size_type l = 0; l < n; ++l)
1655  for (size_type k = 0; k < m; ++k)
1656  for (size_type j = 0; j < n; ++j)
1657  for (size_type i = 0; i < n; ++i, ++it) {
1658  *it = scalar_type(0);
1659  if (l == i) *it += (*(args[0]))[j*m+k]*scalar_type(0.5);
1660  if (l == j) *it += (*(args[0]))[i*m+k]*scalar_type(0.5);
1661  }
1662  GMM_ASSERT1(it == result.end(), "Internal error");
1663  }
1664 
1665  // Second derivative :
1666  // A{ijklop}=(delta{ok}delta{li}delta{pj} + delta{ok}delta{pi}delta{lj})/2
1667  // comes from (H,K) -> (H^{T}K + K^{T}H)/2
1668  void second_derivative(const arg_list &args, size_type, size_type,
1669  base_tensor &result) const { // to be verified
1670  size_type m = args[0]->sizes()[0], n = args[0]->sizes()[1];
1671  base_tensor::iterator it = result.begin();
1672  for (size_type p = 0; p < n; ++p)
1673  for (size_type o = 0; o < m; ++o)
1674  for (size_type l = 0; l < n; ++l)
1675  for (size_type k = 0; k < m; ++k)
1676  for (size_type j = 0; j < n; ++j)
1677  for (size_type i = 0; i < n; ++i, ++it) {
1678  *it = scalar_type(0);
1679  if (o == k) {
1680  if (l == i && p == j) *it += scalar_type(0.5);
1681  if (p == i && l == j) *it += scalar_type(0.5);
1682  }
1683  }
1684  GMM_ASSERT1(it == result.end(), "Internal error");
1685  }
1686  };
1687 
1688  // Cauchy stress tensor from second Piola-Kirchhoff stress tensor
1689  // (I+Grad_u)\sigma(I+Grad_u')/det(I+Grad_u)
1690  struct Cauchy_stress_from_PK2 : public ga_nonlinear_operator {
1691  bool result_size(const arg_list &args, bgeot::multi_index &sizes) const {
1692  if (args.size() != 2 || args[0]->sizes().size() != 2
1693  || args[1]->sizes().size() != 2
1694  || args[0]->sizes()[0] != args[0]->sizes()[1]
1695  || args[1]->sizes()[0] != args[0]->sizes()[1]
1696  || args[1]->sizes()[1] != args[0]->sizes()[1]) return false;
1697  ga_init_square_matrix_(sizes, args[0]->sizes()[1]);
1698  return true;
1699  }
1700 
1701  // Value : (I+Grad_u)\sigma(I+Grad_u')/det(I+Grad_u)
1702  void value(const arg_list &args, base_tensor &result) const {
1703  size_type N = args[0]->sizes()[0];
1704  base_matrix F(N, N), sigma(N,N), aux(N, N);
1705  gmm::copy(args[0]->as_vector(), sigma.as_vector());
1706  gmm::copy(args[1]->as_vector(), F.as_vector());
1707  gmm::add(gmm::identity_matrix(), F);
1708  gmm::mult(F, sigma, aux);
1709  gmm::mult(aux, gmm::transposed(F), sigma);
1710  scalar_type det = bgeot::lu_det(&(*(F.begin())), N);
1711  gmm::scale(sigma, scalar_type(1)/det);
1712  gmm::copy(sigma.as_vector(), result.as_vector());
1713  }
1714 
1715  // Derivative :
1716  void derivative(const arg_list &args, size_type nder,
1717  base_tensor &result) const { // to be verified
1718  size_type N = args[0]->sizes()[0];
1719  base_matrix F(N, N);
1720  gmm::copy(args[1]->as_vector(), F.as_vector());
1721  gmm::add(gmm::identity_matrix(), F);
1722  scalar_type det = bgeot::lu_det(&(*(F.begin())), N);
1723 
1724  base_tensor::iterator it = result.begin();
1725 
1726  switch (nder) {
1727  case 1: // Derivative with respect to sigma : F_ik F_jl / det(F)
1728  for (size_type l = 0; l < N; ++l)
1729  for (size_type k = 0; k < N; ++k)
1730  for (size_type j = 0; j < N; ++j)
1731  for (size_type i = 0; i < N; ++i, ++it)
1732  *it = F(i,k) * F(j,l) / det;
1733  break;
1734 
1735  case 2:
1736  {
1737  // Derivative with respect to Grad_u :
1738  // (delta_ik sigma_lm F_jm + F_im sigma_mk delta_lj
1739  // - F_im sigma_mn F_jn F^{-1}_lk) / det(F)
1740  base_matrix sigma(N,N), aux(N,N), aux2(N,N);
1741  gmm::copy(args[0]->as_vector(), sigma.as_vector());
1742  gmm::mult(sigma, gmm::transposed(F), aux);
1743  gmm::mult(F, aux, aux2);
1744  bgeot::lu_inverse(&(*(F.begin())), N);
1745  for (size_type l = 0; l < N; ++l)
1746  for (size_type k = 0; k < N; ++k)
1747  for (size_type j = 0; j < N; ++j)
1748  for (size_type i = 0; i < N; ++i, ++it) {
1749  *it = scalar_type(0);
1750  if (i == k) *it += aux(l, j) / det;
1751  if (l == j) *it += aux(k, i) / det;
1752  *it -= aux2(i,j) * F(l,k) / det;
1753  }
1754  }
1755  break;
1756  }
1757  GMM_ASSERT1(it == result.end(), "Internal error");
1758  }
1759 
1760  // Second derivative : to be implemented
1761  void second_derivative(const arg_list &, size_type, size_type,
1762  base_tensor &) const {
1763  GMM_ASSERT1(false, "Sorry, not implemented");
1764  }
1765  };
1766 
1767 
1768  struct AHL_wrapper_sigma : public ga_nonlinear_operator {
1769  phyperelastic_law AHL;
1770  bool result_size(const arg_list &args, bgeot::multi_index &sizes) const {
1771  if (args.size() != 2 || args[0]->sizes().size() != 2
1772  || args[1]->size() != AHL->nb_params()
1773  || args[0]->sizes()[0] != args[0]->sizes()[1]) return false;
1774  ga_init_square_matrix_(sizes, args[0]->sizes()[0]);
1775  return true;
1776  }
1777 
1778  // Value :
1779  void value(const arg_list &args, base_tensor &result) const {
1780  size_type N = args[0]->sizes()[0];
1781  base_vector params(AHL->nb_params());
1782  gmm::copy(args[1]->as_vector(), params);
1783  base_matrix Gu(N, N), E(N,N), sigma(N,N);
1784  gmm::copy(args[0]->as_vector(), Gu.as_vector());
1785  gmm::mult(gmm::transposed(Gu), Gu, E);
1786  gmm::add(Gu, E); gmm::add(gmm::transposed(Gu), E);
1787  gmm::scale(E, scalar_type(0.5));
1788  gmm::add(gmm::identity_matrix(), Gu);
1789  scalar_type det = bgeot::lu_det(&(*(Gu.begin())), N);
1790 
1791  AHL->sigma(E, sigma, params, det);
1792  gmm::copy(sigma.as_vector(), result.as_vector());
1793  }
1794 
1795  // Derivative :
1796  void derivative(const arg_list &args, size_type nder,
1797  base_tensor &result) const {
1798  size_type N = args[0]->sizes()[0];
1799  base_vector params(AHL->nb_params());
1800  gmm::copy(args[1]->as_vector(), params);
1801  base_tensor grad_sigma(N, N, N, N);
1802  base_matrix Gu(N, N), E(N,N);
1803  gmm::copy(args[0]->as_vector(), Gu.as_vector());
1804  gmm::mult(gmm::transposed(Gu), Gu, E);
1805  gmm::add(Gu, E); gmm::add(gmm::transposed(Gu), E);
1806  gmm::scale(E, scalar_type(0.5));
1807  gmm::add(gmm::identity_matrix(), Gu);
1808  scalar_type det = bgeot::lu_det(&(*(Gu.begin())), N);
1809 
1810  GMM_ASSERT1(nder == 1, "Sorry, the derivative of this hyperelastic "
1811  "law with respect to its parameters is not available.");
1812 
1813  AHL->grad_sigma(E, grad_sigma, params, det);
1814 
1815  base_tensor::iterator it = result.begin();
1816  for (size_type l = 0; l < N; ++l)
1817  for (size_type k = 0; k < N; ++k)
1818  for (size_type j = 0; j < N; ++j)
1819  for (size_type i = 0; i < N; ++i, ++it) {
1820  *it = scalar_type(0);
1821  for (size_type m = 0; m < N; ++m)
1822  *it += grad_sigma(i,j,m,l) * Gu(k, m); // to be optimized
1823  }
1824  GMM_ASSERT1(it == result.end(), "Internal error");
1825  }
1826 
1827 
1828  // Second derivative : not implemented (not necessary)
1829  void second_derivative(const arg_list &, size_type, size_type,
1830  base_tensor &) const {
1831  GMM_ASSERT1(false, "Sorry, second derivative not implemented");
1832  }
1833 
1834  AHL_wrapper_sigma(const phyperelastic_law &A) : AHL(A) {}
1835 
1836  };
1837 
1838 
1839  struct AHL_wrapper_potential : public ga_nonlinear_operator {
1840  phyperelastic_law AHL;
1841  bool result_size(const arg_list &args, bgeot::multi_index &sizes) const {
1842  if (args.size() != 2 || args[0]->sizes().size() != 2
1843  || args[1]->size() != AHL->nb_params()
1844  || args[0]->sizes()[0] != args[0]->sizes()[1]) return false;
1845  ga_init_scalar_(sizes);
1846  return true;
1847  }
1848 
1849  // Value :
1850  void value(const arg_list &args, base_tensor &result) const {
1851  size_type N = args[0]->sizes()[0];
1852  base_vector params(AHL->nb_params());
1853  gmm::copy(args[1]->as_vector(), params);
1854  base_matrix Gu(N, N), E(N,N);
1855  gmm::copy(args[0]->as_vector(), Gu.as_vector());
1856  gmm::mult(gmm::transposed(Gu), Gu, E);
1857  gmm::add(Gu, E); gmm::add(gmm::transposed(Gu), E);
1858  gmm::scale(E, scalar_type(0.5));
1859  gmm::add(gmm::identity_matrix(), Gu);
1860  scalar_type det = bgeot::lu_det(&(*(Gu.begin())), N); // not necessary here
1861 
1862  result[0] = AHL->strain_energy(E, params, det);
1863  }
1864 
1865  // Derivative :
1866  void derivative(const arg_list &args, size_type nder,
1867  base_tensor &result) const {
1868  size_type N = args[0]->sizes()[0];
1869  base_vector params(AHL->nb_params());
1870  gmm::copy(args[1]->as_vector(), params);
1871  base_matrix Gu(N, N), E(N,N), sigma(N,N);
1872  gmm::copy(args[0]->as_vector(), Gu.as_vector());
1873  gmm::mult(gmm::transposed(Gu), Gu, E);
1874  gmm::add(Gu, E); gmm::add(gmm::transposed(Gu), E);
1875  gmm::scale(E, scalar_type(0.5));
1876  gmm::add(gmm::identity_matrix(), Gu);
1877  scalar_type det = bgeot::lu_det(&(*(Gu.begin())), N); // not necessary here
1878 
1879  GMM_ASSERT1(nder == 1, "Sorry, Cannot derive the potential with "
1880  "respect to law parameters.");
1881 
1882  AHL->sigma(E, sigma, params, det);
1883  gmm::mult(Gu, sigma, E);
1884  gmm::copy(E.as_vector(), result.as_vector());
1885  }
1886 
1887 
1888  // Second derivative :
1889  void second_derivative(const arg_list &args, size_type nder1,
1890  size_type nder2, base_tensor &result) const {
1891 
1892  size_type N = args[0]->sizes()[0];
1893  base_vector params(AHL->nb_params());
1894  gmm::copy(args[1]->as_vector(), params);
1895  base_tensor grad_sigma(N, N, N, N);
1896  base_matrix Gu(N, N), E(N,N), sigma(N,N);
1897  gmm::copy(args[0]->as_vector(), Gu.as_vector());
1898  gmm::mult(gmm::transposed(Gu), Gu, E);
1899  gmm::add(Gu, E); gmm::add(gmm::transposed(Gu), E);
1900  gmm::scale(E, scalar_type(0.5));
1901  gmm::add(gmm::identity_matrix(), Gu);
1902  scalar_type det = bgeot::lu_det(&(*(Gu.begin())), N);
1903 
1904  GMM_ASSERT1(nder1 == 1 && nder2 == 1, "Sorry, Cannot derive the "
1905  "potential with respect to law parameters.");
1906 
1907  AHL->sigma(E, sigma, params, det);
1908  AHL->grad_sigma(E, grad_sigma, params, det);
1909 
1910  base_tensor::iterator it = result.begin();
1911  for (size_type l = 0; l < N; ++l)
1912  for (size_type k = 0; k < N; ++k)
1913  for (size_type j = 0; j < N; ++j)
1914  for (size_type i = 0; i < N; ++i, ++it) {
1915  *it = scalar_type(0);
1916  if (i == k) *it += sigma(l,j);
1917 
1918  for (size_type m = 0; m < N; ++m) // to be optimized
1919  for (size_type n = 0; n < N; ++n)
1920  *it += grad_sigma(n,j,m,l) * Gu(k, m) * Gu(i, n);
1921  }
1922  GMM_ASSERT1(it == result.end(), "Internal error");
1923 
1924  }
1925 
1926  AHL_wrapper_potential(const phyperelastic_law &A) : AHL(A) {}
1927 
1928  };
1929 
1930 
1931  // Saint-Venant Kirchhoff tensor:
1932  // lambda Tr(E)Id + 2*mu*E with E = (Grad_u+Grad_u'+Grad_u'Grad_u)/2
1933  struct Saint_Venant_Kirchhoff_sigma : public ga_nonlinear_operator {
1934  bool result_size(const arg_list &args, bgeot::multi_index &sizes) const {
1935  if (args.size() != 2 || args[0]->sizes().size() != 2
1936  || args[1]->size() != 2
1937  || args[0]->sizes()[0] != args[0]->sizes()[1]) return false;
1938  ga_init_square_matrix_(sizes, args[0]->sizes()[0]);
1939  return true;
1940  }
1941 
1942  // Value : Tr(E) + 2*mu*E
1943  void value(const arg_list &args, base_tensor &result) const {
1944  size_type N = args[0]->sizes()[0];
1945  scalar_type lambda = (*(args[1]))[0], mu = (*(args[1]))[1];
1946  base_matrix Gu(N, N), E(N,N);
1947  gmm::copy(args[0]->as_vector(), Gu.as_vector());
1948  gmm::mult(gmm::transposed(Gu), Gu, E);
1949  gmm::add(Gu, E); gmm::add(gmm::transposed(Gu), E);
1950  gmm::scale(E, scalar_type(0.5));
1951  scalar_type trE = gmm::mat_trace(E);
1952 
1953  base_tensor::iterator it = result.begin();
1954  for (size_type j = 0; j < N; ++j)
1955  for (size_type i = 0; i < N; ++i, ++it) {
1956  *it = 2*mu*E(i,j);
1957  if (i == j) *it += lambda*trE;
1958  }
1959  }
1960 
1961  // Derivative / u: lambda Tr(H + H'*U) Id + mu(H + H' + H'*U + U'*H)
1962  // with U = Grad_u, H = Grad_Test_u
1963  // Implementation: A{ijkl} = lambda(delta{kl}delta{ij} + U{kl}delta{ij})
1964  // + mu(delta{ik}delta{jl} + delta{il}delta{jk} + U{kj}delta{il} +
1965  // U{ki}delta{lj})
1966  void derivative(const arg_list &args, size_type nder,
1967  base_tensor &result) const {
1968  size_type N = args[0]->sizes()[0];
1969  scalar_type lambda = (*(args[1]))[0], mu = (*(args[1]))[1], trE;
1970  base_matrix Gu(N, N), E(N,N);
1971  gmm::copy(args[0]->as_vector(), Gu.as_vector());
1972  if (nder > 1) {
1973  gmm::mult(gmm::transposed(Gu), Gu, E);
1974  gmm::add(Gu, E); gmm::add(gmm::transposed(Gu), E);
1975  gmm::scale(E, scalar_type(0.5));
1976  }
1977  base_tensor::iterator it = result.begin();
1978  switch (nder) {
1979  case 1: // Derivative with respect to u
1980  for (size_type l = 0; l < N; ++l)
1981  for (size_type k = 0; k < N; ++k)
1982  for (size_type j = 0; j < N; ++j)
1983  for (size_type i = 0; i < N; ++i, ++it) {
1984  *it = scalar_type(0);
1985  if (i == j && k == l) *it += lambda;
1986  if (i == j) *it += lambda*Gu(k,l);
1987  if (i == k && j == l) *it += mu;
1988  if (i == l && j == k) *it += mu;
1989  if (i == l) *it += mu*Gu(k,j);
1990  if (l == j) *it += mu*Gu(k,i);
1991  }
1992  break;
1993  case 2:
1994  // Derivative with respect to lambda: Tr(E)Id
1995  trE = gmm::mat_trace(E);
1996  for (size_type j = 0; j < N; ++j)
1997  for (size_type i = 0; i < N; ++i, ++it) {
1998  *it = scalar_type(0); if (i == j) *it += trE;
1999  }
2000  // Derivative with respect to mu: 2E
2001  for (size_type j = 0; j < N; ++j)
2002  for (size_type i = 0; i < N; ++i, ++it) {
2003  *it += 2*E(i,j);
2004  }
2005  break;
2006  default: GMM_ASSERT1(false, "Internal error");
2007  }
2008  GMM_ASSERT1(it == result.end(), "Internal error");
2009  }
2010 
2011  // Second derivative : not implemented
2012  void second_derivative(const arg_list &, size_type, size_type,
2013  base_tensor &) const {
2014  GMM_ASSERT1(false, "Sorry, second derivative not implemented");
2015  }
2016  };
2017 
2018 
2019 
2020  static bool init_predef_operators() {
2021 
2022  ga_predef_operator_tab &PREDEF_OPERATORS
2023  = dal::singleton<ga_predef_operator_tab>::instance();
2024 
2025  PREDEF_OPERATORS.add_method
2026  ("Matrix_i2", std::make_shared<matrix_i2_operator>());
2027  PREDEF_OPERATORS.add_method
2028  ("Matrix_j1", std::make_shared<matrix_j1_operator>());
2029  PREDEF_OPERATORS.add_method
2030  ("Matrix_j2", std::make_shared<matrix_j2_operator>());
2031  PREDEF_OPERATORS.add_method
2032  ("Right_Cauchy_Green", std::make_shared<Right_Cauchy_Green_operator>());
2033  PREDEF_OPERATORS.add_method
2034  ("Left_Cauchy_Green", std::make_shared<Left_Cauchy_Green_operator>());
2035  PREDEF_OPERATORS.add_method
2036  ("Green_Lagrangian", std::make_shared<Green_Lagrangian_operator>());
2037 
2038  PREDEF_OPERATORS.add_method
2039  ("Cauchy_stress_from_PK2", std::make_shared<Cauchy_stress_from_PK2>());
2040 
2041  PREDEF_OPERATORS.add_method
2042  ("Saint_Venant_Kirchhoff_sigma",
2043  std::make_shared<Saint_Venant_Kirchhoff_sigma>());
2044  PREDEF_OPERATORS.add_method
2045  ("Saint_Venant_Kirchhoff_PK2",
2046  std::make_shared<Saint_Venant_Kirchhoff_sigma>());
2047  PREDEF_OPERATORS.add_method
2048  ("Saint_Venant_Kirchhoff_potential",
2049  std::make_shared<AHL_wrapper_potential>
2050  (std::make_shared<SaintVenant_Kirchhoff_hyperelastic_law>()));
2051  PREDEF_OPERATORS.add_method
2052  ("Plane_Strain_Saint_Venant_Kirchhoff_sigma",
2053  std::make_shared<Saint_Venant_Kirchhoff_sigma>());
2054  PREDEF_OPERATORS.add_method
2055  ("Plane_Strain_Saint_Venant_Kirchhoff_PK2",
2056  std::make_shared<Saint_Venant_Kirchhoff_sigma>());
2057  PREDEF_OPERATORS.add_method
2058  ("Plane_Strain_Saint_Venant_Kirchhoff_potential",
2059  std::make_shared<AHL_wrapper_potential>
2060  (std::make_shared<SaintVenant_Kirchhoff_hyperelastic_law>()));
2061 
2062  phyperelastic_law gbklaw
2063  = std::make_shared<generalized_Blatz_Ko_hyperelastic_law>();
2064  PREDEF_OPERATORS.add_method
2065  ("Generalized_Blatz_Ko_sigma",
2066  std::make_shared<AHL_wrapper_sigma>(gbklaw));
2067  PREDEF_OPERATORS.add_method
2068  ("Generalized_Blatz_Ko_PK2",
2069  std::make_shared<AHL_wrapper_sigma>(gbklaw));
2070  PREDEF_OPERATORS.add_method
2071  ("Generalized_Blatz_Ko_potential",
2072  std::make_shared<AHL_wrapper_potential>
2073  (std::make_shared<generalized_Blatz_Ko_hyperelastic_law>()));
2074  PREDEF_OPERATORS.add_method
2075  ("Plane_Strain_Generalized_Blatz_Ko_sigma",
2076  std::make_shared<AHL_wrapper_sigma>
2077  (std::make_shared<plane_strain_hyperelastic_law>(gbklaw)));
2078  PREDEF_OPERATORS.add_method
2079  ("Plane_Strain_Generalized_Blatz_Ko_PK2",
2080  std::make_shared<AHL_wrapper_sigma>
2081  (std::make_shared<plane_strain_hyperelastic_law>(gbklaw)));
2082  PREDEF_OPERATORS.add_method
2083  ("Plane_Strain_Generalized_Blatz_Ko_potential",
2084  std::make_shared<AHL_wrapper_potential>
2085  (std::make_shared<plane_strain_hyperelastic_law>(gbklaw)));
2086 
2087  phyperelastic_law cigelaw
2088  = std::make_shared<Ciarlet_Geymonat_hyperelastic_law>();
2089  PREDEF_OPERATORS.add_method
2090  ("Ciarlet_Geymonat_PK2", std::make_shared<AHL_wrapper_sigma>(cigelaw));
2091  PREDEF_OPERATORS.add_method
2092  ("Ciarlet_Geymonat_sigma", std::make_shared<AHL_wrapper_sigma>(cigelaw));
2093  PREDEF_OPERATORS.add_method
2094  ("Ciarlet_Geymonat_potential",
2095  std::make_shared<AHL_wrapper_potential>
2096  (std::make_shared<Ciarlet_Geymonat_hyperelastic_law>()));
2097  PREDEF_OPERATORS.add_method
2098  ("Plane_Strain_Ciarlet_Geymonat_sigma",
2099  std::make_shared<AHL_wrapper_sigma>
2100  (std::make_shared<plane_strain_hyperelastic_law>(cigelaw)));
2101  PREDEF_OPERATORS.add_method
2102  ("Plane_Strain_Ciarlet_Geymonat_PK2",
2103  std::make_shared<AHL_wrapper_sigma>
2104  (std::make_shared<plane_strain_hyperelastic_law>(cigelaw)));
2105  PREDEF_OPERATORS.add_method
2106  ("Plane_Strain_Ciarlet_Geymonat_potential",
2107  std::make_shared<AHL_wrapper_potential>
2108  (std::make_shared<plane_strain_hyperelastic_law>(cigelaw)));
2109 
2110  phyperelastic_law morilaw
2111  = std::make_shared<Mooney_Rivlin_hyperelastic_law>();
2112  PREDEF_OPERATORS.add_method
2113  ("Incompressible_Mooney_Rivlin_sigma",
2114  std::make_shared<AHL_wrapper_sigma>(morilaw));
2115  PREDEF_OPERATORS.add_method
2116  ("Incompressible_Mooney_Rivlin_PK2",
2117  std::make_shared<AHL_wrapper_sigma>(morilaw));
2118  PREDEF_OPERATORS.add_method
2119  ("Incompressible_Mooney_Rivlin_potential",
2120  std::make_shared<AHL_wrapper_potential>
2121  (std::make_shared<Mooney_Rivlin_hyperelastic_law>()));
2122  PREDEF_OPERATORS.add_method
2123  ("Plane_Strain_Incompressible_Mooney_Rivlin_PK2",
2124  std::make_shared<AHL_wrapper_sigma>
2125  (std::make_shared<plane_strain_hyperelastic_law>(morilaw)));
2126  PREDEF_OPERATORS.add_method
2127  ("Plane_Strain_Incompressible_Mooney_Rivlin_sigma",
2128  std::make_shared<AHL_wrapper_sigma>
2129  (std::make_shared<plane_strain_hyperelastic_law>(morilaw)));
2130  PREDEF_OPERATORS.add_method
2131  ("Plane_Strain_Incompressible_Mooney_Rivlin_potential",
2132  std::make_shared<AHL_wrapper_potential>
2133  (std::make_shared<plane_strain_hyperelastic_law>(morilaw)));
2134 
2135  phyperelastic_law cmorilaw
2136  = std::make_shared<Mooney_Rivlin_hyperelastic_law>(true);
2137  PREDEF_OPERATORS.add_method
2138  ("Compressible_Mooney_Rivlin_sigma",
2139  std::make_shared<AHL_wrapper_sigma>(cmorilaw));
2140  PREDEF_OPERATORS.add_method
2141  ("Compressible_Mooney_Rivlin_PK2",
2142  std::make_shared<AHL_wrapper_sigma>(cmorilaw));
2143  PREDEF_OPERATORS.add_method
2144  ("Compressible_Mooney_Rivlin_potential",
2145  std::make_shared<AHL_wrapper_potential>
2146  (std::make_shared<Mooney_Rivlin_hyperelastic_law>(true)));
2147  PREDEF_OPERATORS.add_method
2148  ("Plane_Strain_Compressible_Mooney_Rivlin_PK2",
2149  std::make_shared<AHL_wrapper_sigma>
2150  (std::make_shared<plane_strain_hyperelastic_law>(cmorilaw)));
2151  PREDEF_OPERATORS.add_method
2152  ("Plane_Strain_Compressible_Mooney_Rivlin_sigma",
2153  std::make_shared<AHL_wrapper_sigma>
2154  (std::make_shared<plane_strain_hyperelastic_law>(cmorilaw)));
2155  PREDEF_OPERATORS.add_method
2156  ("Plane_Strain_Compressible_Mooney_Rivlin_potential",
2157  std::make_shared<AHL_wrapper_potential>
2158  (std::make_shared<plane_strain_hyperelastic_law>(cmorilaw)));
2159 
2160  phyperelastic_law ineolaw
2161  = std::make_shared<Mooney_Rivlin_hyperelastic_law>(false, true);
2162  PREDEF_OPERATORS.add_method
2163  ("Incompressible_Neo_Hookean_sigma",
2164  std::make_shared<AHL_wrapper_sigma>(ineolaw));
2165  PREDEF_OPERATORS.add_method
2166  ("Incompressible_Neo_Hookean_PK2",
2167  std::make_shared<AHL_wrapper_sigma>(ineolaw));
2168  PREDEF_OPERATORS.add_method
2169  ("Incompressible_Neo_Hookean_potential",
2170  std::make_shared<AHL_wrapper_potential>
2171  (std::make_shared<Mooney_Rivlin_hyperelastic_law>(false,true)));
2172  PREDEF_OPERATORS.add_method
2173  ("Plane_Strain_Incompressible_Neo_Hookean_sigma",
2174  std::make_shared<AHL_wrapper_sigma>
2175  (std::make_shared<plane_strain_hyperelastic_law>(ineolaw)));
2176  PREDEF_OPERATORS.add_method
2177  ("Plane_Strain_Incompressible_Neo_Hookean_PK2",
2178  std::make_shared<AHL_wrapper_sigma>
2179  (std::make_shared<plane_strain_hyperelastic_law>(ineolaw)));
2180  PREDEF_OPERATORS.add_method
2181  ("Plane_Strain_Incompressible_Neo_Hookean_potential",
2182  std::make_shared<AHL_wrapper_potential>
2183  (std::make_shared<plane_strain_hyperelastic_law>(ineolaw)));
2184 
2185  phyperelastic_law cneolaw
2186  = std::make_shared<Mooney_Rivlin_hyperelastic_law>(true, true);
2187  PREDEF_OPERATORS.add_method
2188  ("Compressible_Neo_Hookean_sigma",
2189  std::make_shared<AHL_wrapper_sigma>(cneolaw));
2190  PREDEF_OPERATORS.add_method
2191  ("Compressible_Neo_Hookean_PK2",
2192  std::make_shared<AHL_wrapper_sigma>(cneolaw));
2193  PREDEF_OPERATORS.add_method
2194  ("Compressible_Neo_Hookean_potential",
2195  std::make_shared<AHL_wrapper_potential>
2196  (std::make_shared<Mooney_Rivlin_hyperelastic_law>(true,true)));
2197  PREDEF_OPERATORS.add_method
2198  ("Plane_Strain_Compressible_Neo_Hookean_sigma",
2199  std::make_shared<AHL_wrapper_sigma>
2200  (std::make_shared<plane_strain_hyperelastic_law>(cneolaw)));
2201  PREDEF_OPERATORS.add_method
2202  ("Plane_Strain_Compressible_Neo_Hookean_PK2",
2203  std::make_shared<AHL_wrapper_sigma>
2204  (std::make_shared<plane_strain_hyperelastic_law>(cneolaw)));
2205  PREDEF_OPERATORS.add_method
2206  ("Plane_Strain_Compressible_Neo_Hookean_potential",
2207  std::make_shared<AHL_wrapper_potential>
2208  (std::make_shared<plane_strain_hyperelastic_law>(cneolaw)));
2209 
2210  phyperelastic_law cneobolaw
2211  = std::make_shared<Neo_Hookean_hyperelastic_law>(true);
2212  PREDEF_OPERATORS.add_method
2213  ("Compressible_Neo_Hookean_Bonet_sigma",
2214  std::make_shared<AHL_wrapper_sigma>(cneobolaw));
2215  PREDEF_OPERATORS.add_method
2216  ("Compressible_Neo_Hookean_Bonet_PK2",
2217  std::make_shared<AHL_wrapper_sigma>(cneobolaw));
2218  PREDEF_OPERATORS.add_method
2219  ("Compressible_Neo_Hookean_Bonet_potential",
2220  std::make_shared<AHL_wrapper_potential>
2221  (std::make_shared<Neo_Hookean_hyperelastic_law>(true)));
2222  PREDEF_OPERATORS.add_method
2223  ("Plane_Strain_Compressible_Neo_Hookean_Bonet_sigma",
2224  std::make_shared<AHL_wrapper_sigma>
2225  (std::make_shared<plane_strain_hyperelastic_law>(cneobolaw)));
2226  PREDEF_OPERATORS.add_method
2227  ("Plane_Strain_Compressible_Neo_Hookean_Bonet_PK2",
2228  std::make_shared<AHL_wrapper_sigma>
2229  (std::make_shared<plane_strain_hyperelastic_law>(cneobolaw)));
2230  PREDEF_OPERATORS.add_method
2231  ("Plane_Strain_Compressible_Neo_Hookean_Bonet_potential",
2232  std::make_shared<AHL_wrapper_potential>
2233  (std::make_shared<plane_strain_hyperelastic_law>(cneobolaw)));
2234 
2235  phyperelastic_law cneocilaw
2236  = std::make_shared<Neo_Hookean_hyperelastic_law>(false);
2237  PREDEF_OPERATORS.add_method
2238  ("Compressible_Neo_Hookean_Ciarlet_sigma",
2239  std::make_shared<AHL_wrapper_sigma>(cneocilaw));
2240  PREDEF_OPERATORS.add_method
2241  ("Compressible_Neo_Hookean_Ciarlet_PK2",
2242  std::make_shared<AHL_wrapper_sigma>(cneocilaw));
2243  PREDEF_OPERATORS.add_method
2244  ("Compressible_Neo_Hookean_Ciarlet_potential",
2245  std::make_shared<AHL_wrapper_potential>
2246  (std::make_shared<Neo_Hookean_hyperelastic_law>(false)));
2247  PREDEF_OPERATORS.add_method
2248  ("Plane_Strain_Compressible_Neo_Hookean_Ciarlet_sigma",
2249  std::make_shared<AHL_wrapper_sigma>
2250  (std::make_shared<plane_strain_hyperelastic_law>(cneocilaw)));
2251  PREDEF_OPERATORS.add_method
2252  ("Plane_Strain_Compressible_Neo_Hookean_Ciarlet_PK2",
2253  std::make_shared<AHL_wrapper_sigma>
2254  (std::make_shared<plane_strain_hyperelastic_law>(cneocilaw)));
2255  PREDEF_OPERATORS.add_method
2256  ("Plane_Strain_Compressible_Neo_Hookean_Ciarlet_potential",
2257  std::make_shared<AHL_wrapper_potential>
2258  (std::make_shared<plane_strain_hyperelastic_law>(cneocilaw)));
2259 
2260  return true;
2261  }
2262 
2263  // declared in getfem_generic_assembly.cc
2264  bool predef_operators_nonlinear_elasticity_initialized
2265  = init_predef_operators();
2266 
2267 
2268  std::string adapt_law_name(const std::string &lawname, size_type N) {
2269  std::string adapted_lawname = lawname;
2270 
2271  for (size_type i = 0; i < lawname.size(); ++i)
2272  if (adapted_lawname[i] == ' ') adapted_lawname[i] = '_';
2273 
2274  if (adapted_lawname.compare("SaintVenant_Kirchhoff") == 0) {
2275  adapted_lawname = "Saint_Venant_Kirchhoff";
2276  } else if (adapted_lawname.compare("Saint_Venant_Kirchhoff") == 0) {
2277 
2278  } else if (adapted_lawname.compare("Generalized_Blatz_Ko") == 0) {
2279  if (N == 2) adapted_lawname = "Plane_Strain_" + adapted_lawname;
2280  } else if (adapted_lawname.compare("Ciarlet_Geymonat") == 0) {
2281  if (N == 2) adapted_lawname = "Plane_Strain_" + adapted_lawname;
2282  } else if (adapted_lawname.compare("Incompressible_Mooney_Rivlin") == 0) {
2283  if (N == 2) adapted_lawname = "Plane_Strain_" + adapted_lawname;
2284  } else if (adapted_lawname.compare("Compressible_Mooney_Rivlin") == 0) {
2285  if (N == 2) adapted_lawname = "Plane_Strain_" + adapted_lawname;
2286  } else if (adapted_lawname.compare("Incompressible_Neo_Hookean") == 0) {
2287  if (N == 2) adapted_lawname = "Plane_Strain_" + adapted_lawname;
2288  } else if (adapted_lawname.compare("Compressible_Neo_Hookean") == 0 ||
2289  adapted_lawname.compare("Compressible_Neo_Hookean_Bonet") == 0 ||
2290  adapted_lawname.compare("Compressible_Neo_Hookean_Ciarlet") == 0 ) {
2291  if (N == 2) adapted_lawname = "Plane_Strain_" + adapted_lawname;
2292  } else
2293  GMM_ASSERT1(false, lawname << " is not a known hyperelastic law");
2294 
2295  return adapted_lawname;
2296  }
2297 
2298 
2299  size_type add_finite_strain_elasticity_brick
2300  (model &md, const mesh_im &mim, const std::string &lawname,
2301  const std::string &varname, const std::string &params,
2302  size_type region) {
2303  std::string test_varname = "Test_" + sup_previous_and_dot_to_varname(varname);
2304  size_type N = mim.linked_mesh().dim();
2305  GMM_ASSERT1(N >= 2 && N <= 3,
2306  "Finite strain elasticity brick works only in 2D or 3D");
2307 
2308  const mesh_fem *mf = md.pmesh_fem_of_variable(varname);
2309  GMM_ASSERT1(mf, "Finite strain elasticity brick can only be applied on "
2310  "fem variables");
2311  size_type Q = mf->get_qdim();
2312  GMM_ASSERT1(Q == N, "Finite strain elasticity brick can only be applied "
2313  "on a fem variable having the same dimension as the mesh");
2314 
2315  std::string adapted_lawname = adapt_law_name(lawname, N);
2316 
2317  std::string expr = "((Id(meshdim)+Grad_"+varname+")*(" + adapted_lawname
2318  + "_PK2(Grad_"+varname+","+params+"))):Grad_" + test_varname;
2319 
2320  return add_nonlinear_generic_assembly_brick
2321  (md, mim, expr, region, true, false,
2322  "Finite strain elasticity brick for " + adapted_lawname + " law");
2323  }
2324 
2325  size_type add_finite_strain_incompressibility_brick
2326  (model &md, const mesh_im &mim, const std::string &varname,
2327  const std::string &multname, size_type region) {
2328  std::string test_varname = "Test_" + sup_previous_and_dot_to_varname(varname);
2329  std::string test_multname = "Test_" + sup_previous_and_dot_to_varname(multname);
2330 
2331  std::string expr
2332  = "(" + test_multname+ ")*(1-Det(Id(meshdim)+Grad_" + varname + "))"
2333  + "-(" + multname + ")*(Det(Id(meshdim)+Grad_" + varname + ")"
2334  + "*((Inv(Id(meshdim)+Grad_" + varname + "))':Grad_"
2335  + test_varname + "))" ;
2336  return add_nonlinear_generic_assembly_brick
2337  (md, mim, expr, region, true, false,
2338  "Finite strain incompressibility brick");
2339  }
2340 
2341  void compute_finite_strain_elasticity_Von_Mises
2342  (model &md, const std::string &lawname, const std::string &varname,
2343  const std::string &params, const mesh_fem &mf_vm,
2344  model_real_plain_vector &VM, const mesh_region &rg) {
2345  size_type N = mf_vm.linked_mesh().dim();
2346 
2347  std::string adapted_lawname = adapt_law_name(lawname, N);
2348 
2349  std::string expr = "sqrt(3/2)*Norm(Deviator(Cauchy_stress_from_PK2("
2350  + adapted_lawname + "_PK2(Grad_" + varname + "," + params + "),Grad_"
2351  + varname + ")))";
2352  ga_interpolation_Lagrange_fem(md, expr, mf_vm, VM, rg);
2353  }
2354 
2355 
2356 
2357 } /* end of namespace getfem. */
2358 
dal::singleton::instance
static T & instance()
Instance from the current thread.
Definition: dal_singleton.h:165
getfem::SaintVenant_Kirchhoff_hyperelastic_law::grad_sigma_updated_lagrangian
virtual void grad_sigma_updated_lagrangian(const base_matrix &F, const base_matrix &E, const base_vector &params, scalar_type det_trans, base_tensor &grad_sigma_ul) const
cauchy-truesdel tangent moduli, used in updated lagrangian
Definition: getfem_nonlinear_elasticity.cc:423
getfem_nonlinear_elasticity.h
Non-linear elasticty and incompressibility bricks.
getfem::mesh_fem::get_qdim
virtual dim_type get_qdim() const
Return the Q dimension.
Definition: getfem_mesh_fem.h:312
getfem::add_nonlinear_incompressibility_brick
size_type add_nonlinear_incompressibility_brick(model &md, const mesh_im &mim, const std::string &varname, const std::string &multname, size_type region=size_type(-1))
Add a nonlinear incompressibility term (for large strain elasticity) to the model with respect to the...
Definition: getfem_nonlinear_elasticity.cc:1270
getfem::abstract_hyperelastic_law::grad_sigma_updated_lagrangian
virtual void grad_sigma_updated_lagrangian(const base_matrix &F, const base_matrix &E, const base_vector &params, scalar_type det_trans, base_tensor &grad_sigma_ul) const
cauchy-truesdel tangent moduli, used in updated lagrangian
Definition: getfem_nonlinear_elasticity.cc:363
gmm::resize
void resize(M &v, size_type m, size_type n)
*‍/
Definition: gmm_blas.h:231
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
getfem::compute_gradient
void compute_gradient(const mesh_fem &mf, const mesh_fem &mf_target, const VECT1 &UU, VECT2 &VV)
Compute the gradient of a field on a getfem::mesh_fem.
Definition: getfem_derivatives.h:62
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:59
getfem::abstract_hyperelastic_law::cauchy_updated_lagrangian
virtual void cauchy_updated_lagrangian(const base_matrix &F, const base_matrix &E, base_matrix &cauchy_stress, const base_vector &params, scalar_type det_trans) const
True Cauchy stress (for Updated Lagrangian formulation)
Definition: getfem_nonlinear_elasticity.cc:348
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:2897
getfem::mesh_im
Describe an integration method linked to a mesh.
Definition: getfem_mesh_im.h:47
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::add_finite_strain_incompressibility_brick
size_type add_finite_strain_incompressibility_brick(model &md, const mesh_im &mim, const std::string &varname, const std::string &multname, size_type region=size_type(-1))
Add a finite strain incompressibility term (for large strain elasticity) to the model with respect to...
Definition: getfem_nonlinear_elasticity.cc:2326
getfem::mesh_fem
Describe a finite element method linked to a mesh.
Definition: getfem_mesh_fem.h:148
gmm::mat_euclidean_norm_sqr
number_traits< typename linalg_traits< M >::value_type >::magnitude_type mat_euclidean_norm_sqr(const M &m)
*‍/
Definition: gmm_blas.h:626
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::model
`‘Model’' variables store the variables, the data and the description of a model.
Definition: getfem_models.h:114
getfem::compute_finite_strain_elasticity_Von_Mises
void compute_finite_strain_elasticity_Von_Mises(model &md, const std::string &lawname, const std::string &varname, const std::string &params, const mesh_fem &mf_vm, model_real_plain_vector &VM, const mesh_region &rg=mesh_region::all_convexes())
Interpolate the Von-Mises stress of a field varname with respect to the nonlinear elasticity constitu...
Definition: getfem_nonlinear_elasticity.cc:2342
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
getfem_models.h
Model representation in Getfem.
getfem::asm_nonlinear_elasticity_tangent_matrix
void asm_nonlinear_elasticity_tangent_matrix(const MAT &K_, const mesh_im &mim, const getfem::mesh_fem &mf, const VECT1 &U, const getfem::mesh_fem *mf_data, const VECT2 &PARAMS, const abstract_hyperelastic_law &AHL, const mesh_region &rg=mesh_region::all_convexes())
Tangent matrix for the non-linear elasticity.
Definition: getfem_nonlinear_elasticity.h:341
getfem_generic_assembly.h
A language for generic assembly of pde boundary value problems.
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::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:1033
getfem::mesh_fem::linked_mesh
const mesh & linked_mesh() const
Return a reference to the underlying mesh.
Definition: getfem_mesh_fem.h:279
getfem::add_finite_strain_elasticity_brick
size_type add_finite_strain_elasticity_brick(model &md, const mesh_im &mim, const std::string &lawname, const std::string &varname, const std::string &params, size_type region=size_type(-1))
Add a finite strain elasticity brick to the model with respect to the variable varname (the displacem...
Definition: getfem_nonlinear_elasticity.cc:2300
gmm::nnz
size_type nnz(const L &l)
count the number of non-zero entries of a vector or matrix.
Definition: gmm_blas.h:68
gmm::fill_random
void fill_random(L &l, double cfill)
*‍/
Definition: gmm_blas.h:154
bgeot::alpha
size_type alpha(short_type n, short_type d)
Return the value of which is the number of monomials of a polynomial of variables and degree .
Definition: bgeot_poly.cc:47
gmm::mat_trace
linalg_traits< M >::value_type mat_trace(const M &m)
Trace of a matrix.
Definition: gmm_blas.h:528
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::mesh_region
structure used to hold a set of convexes and/or convex faces.
Definition: getfem_mesh_region.h:55
getfem::pbrick
std::shared_ptr< const virtual_brick > pbrick
type of pointer on a brick
Definition: getfem_models.h:49
getfem::add_nonlinear_elasticity_brick
size_type add_nonlinear_elasticity_brick(model &md, const mesh_im &mim, const std::string &varname, const phyperelastic_law &AHL, const std::string &dataname, size_type region=size_type(-1))
Add a nonlinear (large strain) elasticity term to the model with respect to the variable varname (dep...
Definition: getfem_nonlinear_elasticity.cc:1020