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"
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
82  { if (!di1_c) compute_di1(); return di1; }
83
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());
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
121  { if (!di2_c) compute_di2(); return di2; }
122
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
159  { if (!di3_c) compute_di3(); return di3; }
160
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() {
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();
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
202  { if (!dj1_c) compute_dj1(); return dj1; }
203
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() {
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);
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
247  { if (!dj2_c) compute_dj2(); return dj2; }
248
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);
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);
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);
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
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
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
363  (const base_matrix& F, const base_matrix& E,
364  const base_vector &params, scalar_type det_trans,
366  {
367  size_type N = E.ncols();
368  base_tensor Cse(N,N,N,N);
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  {
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)
384  F(i,m)*F(j,n)*F(k,p)*F(l,q)*Cse(m,n,p,q);
385  }
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)) {
407  }
408  }
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
424  const base_matrix& E,
425  const base_vector &params,
426  scalar_type det_trans,
428  {
429  size_type N = E.ncols();
430  base_tensor Cse(N,N,N,N);
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);
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
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));
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));
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
544  }
545  if (compressible) {
546  scalar_type D1 = params[i++];
547  scalar_type di3 = D1 - D1 / sqrt(gmm::abs(ci.i3()));
549
550 // shouldn't negative det_trans be handled here???
551  if (det_trans <= scalar_type(0))
553  }
554  }
555
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));
565  compute_invariants ci(C);
566
567  size_type i=0;
568  scalar_type C1 = params[i++]; // C10
570  scalar_type(4)*C1), result.as_vector());
571  if (!neohookean) {
572  scalar_type C2 = params[i++]; // C01
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()));
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));
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));
642  compute_invariants ci(C);
643
644  scalar_type lambda = params[0];
645  scalar_type mu = params[1];
647  if (bonet)
649  (lambda/2 * log(ci.i3()) - mu) / ci.i3()), result);
650  else
652  lambda/2 - lambda/(2*ci.i3()) - mu / ci.i3()), result);
653  if (det_trans <= scalar_type(0))
655  }
656
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));
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());
675  (lambda * logi3 - 2*mu) / ci.i3()),
676  result.as_vector());
677  coeff = (lambda + 2 * mu - lambda * logi3) / gmm::sqr(ci.i3());
678  } else {
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));
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));
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);
745  if (det_trans <= scalar_type(0))
747
748  }
749
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));
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
773  scalar_type(4)*di1), result.as_vector());
775  scalar_type(4)*di2), result.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];
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);
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);
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))
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
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);
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
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);
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);
988
989  if (version & model::BUILD_MATRIX) {
990  gmm::clear(matl[0]);
991  GMM_TRACE2("Nonlinear elasticity stiffness matrix assembly");
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
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
1056  unsigned NP = unsigned(AHL->nb_params()), NFem = mf_u.get_qdim();
1058  model_real_plain_vector PARAMS(mf_vm.nb_dof()*NP);
1059  if (mf_params) interpolation(*mf_params, mf_vm, params, PARAMS);
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) {
1076  for (unsigned int alpha = 0; alpha <N; ++alpha)
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){
1091  }
1092  normEz= gmm::vect_norm2(ez);
1093  }
1094  //jyh : end compute ez
1095  }
1098
1099  /* jyh : complete gradphi for virtual 3rd dim (perpendicular to
1100  deformed surface, same thickness) */
1101  if (NFem == 3 && N == 2) {
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)
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
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
1157  model_real_plain_vector PARAMS(mf_sigma.nb_dof()*NP);
1158
1159
1161  if (mf_params) {
1162  for (size_type i = 0; i < mf_sigma.nb_dof(); ++i)
1164  interpolation(*mf_params, mti, params, PARAMS);
1165  }
1166
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
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
1193  for (unsigned int alpha = 0; alpha <N; ++alpha)
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);
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
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
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
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());
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());
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);
1787  gmm::scale(E, scalar_type(0.5));
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);
1806  gmm::scale(E, scalar_type(0.5));
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
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);
1858  gmm::scale(E, scalar_type(0.5));
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);
1875  gmm::scale(E, scalar_type(0.5));
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);
1900  gmm::scale(E, scalar_type(0.5));
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);
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:
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);
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)
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);
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
2024
2026  ("Matrix_i2", std::make_shared<matrix_i2_operator>());
2028  ("Matrix_j1", std::make_shared<matrix_j1_operator>());
2030  ("Matrix_j2", std::make_shared<matrix_j2_operator>());
2032  ("Right_Cauchy_Green", std::make_shared<Right_Cauchy_Green_operator>());
2034  ("Left_Cauchy_Green", std::make_shared<Left_Cauchy_Green_operator>());
2036  ("Green_Lagrangian", std::make_shared<Green_Lagrangian_operator>());
2037
2039  ("Cauchy_stress_from_PK2", std::make_shared<Cauchy_stress_from_PK2>());
2040
2042  ("Saint_Venant_Kirchhoff_sigma",
2043  std::make_shared<Saint_Venant_Kirchhoff_sigma>());
2045  ("Saint_Venant_Kirchhoff_PK2",
2046  std::make_shared<Saint_Venant_Kirchhoff_sigma>());
2048  ("Saint_Venant_Kirchhoff_potential",
2049  std::make_shared<AHL_wrapper_potential>
2050  (std::make_shared<SaintVenant_Kirchhoff_hyperelastic_law>()));
2052  ("Plane_Strain_Saint_Venant_Kirchhoff_sigma",
2053  std::make_shared<Saint_Venant_Kirchhoff_sigma>());
2055  ("Plane_Strain_Saint_Venant_Kirchhoff_PK2",
2056  std::make_shared<Saint_Venant_Kirchhoff_sigma>());
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>();
2065  ("Generalized_Blatz_Ko_sigma",
2066  std::make_shared<AHL_wrapper_sigma>(gbklaw));
2068  ("Generalized_Blatz_Ko_PK2",
2069  std::make_shared<AHL_wrapper_sigma>(gbklaw));
2071  ("Generalized_Blatz_Ko_potential",
2072  std::make_shared<AHL_wrapper_potential>
2073  (std::make_shared<generalized_Blatz_Ko_hyperelastic_law>()));
2075  ("Plane_Strain_Generalized_Blatz_Ko_sigma",
2076  std::make_shared<AHL_wrapper_sigma>
2077  (std::make_shared<plane_strain_hyperelastic_law>(gbklaw)));
2079  ("Plane_Strain_Generalized_Blatz_Ko_PK2",
2080  std::make_shared<AHL_wrapper_sigma>
2081  (std::make_shared<plane_strain_hyperelastic_law>(gbklaw)));
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>();
2090  ("Ciarlet_Geymonat_PK2", std::make_shared<AHL_wrapper_sigma>(cigelaw));
2092  ("Ciarlet_Geymonat_sigma", std::make_shared<AHL_wrapper_sigma>(cigelaw));
2094  ("Ciarlet_Geymonat_potential",
2095  std::make_shared<AHL_wrapper_potential>
2096  (std::make_shared<Ciarlet_Geymonat_hyperelastic_law>()));
2098  ("Plane_Strain_Ciarlet_Geymonat_sigma",
2099  std::make_shared<AHL_wrapper_sigma>
2100  (std::make_shared<plane_strain_hyperelastic_law>(cigelaw)));
2102  ("Plane_Strain_Ciarlet_Geymonat_PK2",
2103  std::make_shared<AHL_wrapper_sigma>
2104  (std::make_shared<plane_strain_hyperelastic_law>(cigelaw)));
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>();
2113  ("Incompressible_Mooney_Rivlin_sigma",
2114  std::make_shared<AHL_wrapper_sigma>(morilaw));
2116  ("Incompressible_Mooney_Rivlin_PK2",
2117  std::make_shared<AHL_wrapper_sigma>(morilaw));
2119  ("Incompressible_Mooney_Rivlin_potential",
2120  std::make_shared<AHL_wrapper_potential>
2121  (std::make_shared<Mooney_Rivlin_hyperelastic_law>()));
2123  ("Plane_Strain_Incompressible_Mooney_Rivlin_PK2",
2124  std::make_shared<AHL_wrapper_sigma>
2125  (std::make_shared<plane_strain_hyperelastic_law>(morilaw)));
2127  ("Plane_Strain_Incompressible_Mooney_Rivlin_sigma",
2128  std::make_shared<AHL_wrapper_sigma>
2129  (std::make_shared<plane_strain_hyperelastic_law>(morilaw)));
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);
2138  ("Compressible_Mooney_Rivlin_sigma",
2139  std::make_shared<AHL_wrapper_sigma>(cmorilaw));
2141  ("Compressible_Mooney_Rivlin_PK2",
2142  std::make_shared<AHL_wrapper_sigma>(cmorilaw));
2144  ("Compressible_Mooney_Rivlin_potential",
2145  std::make_shared<AHL_wrapper_potential>
2146  (std::make_shared<Mooney_Rivlin_hyperelastic_law>(true)));
2148  ("Plane_Strain_Compressible_Mooney_Rivlin_PK2",
2149  std::make_shared<AHL_wrapper_sigma>
2150  (std::make_shared<plane_strain_hyperelastic_law>(cmorilaw)));
2152  ("Plane_Strain_Compressible_Mooney_Rivlin_sigma",
2153  std::make_shared<AHL_wrapper_sigma>
2154  (std::make_shared<plane_strain_hyperelastic_law>(cmorilaw)));
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);
2163  ("Incompressible_Neo_Hookean_sigma",
2164  std::make_shared<AHL_wrapper_sigma>(ineolaw));
2166  ("Incompressible_Neo_Hookean_PK2",
2167  std::make_shared<AHL_wrapper_sigma>(ineolaw));
2169  ("Incompressible_Neo_Hookean_potential",
2170  std::make_shared<AHL_wrapper_potential>
2171  (std::make_shared<Mooney_Rivlin_hyperelastic_law>(false,true)));
2173  ("Plane_Strain_Incompressible_Neo_Hookean_sigma",
2174  std::make_shared<AHL_wrapper_sigma>
2175  (std::make_shared<plane_strain_hyperelastic_law>(ineolaw)));
2177  ("Plane_Strain_Incompressible_Neo_Hookean_PK2",
2178  std::make_shared<AHL_wrapper_sigma>
2179  (std::make_shared<plane_strain_hyperelastic_law>(ineolaw)));
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);
2188  ("Compressible_Neo_Hookean_sigma",
2189  std::make_shared<AHL_wrapper_sigma>(cneolaw));
2191  ("Compressible_Neo_Hookean_PK2",
2192  std::make_shared<AHL_wrapper_sigma>(cneolaw));
2194  ("Compressible_Neo_Hookean_potential",
2195  std::make_shared<AHL_wrapper_potential>
2196  (std::make_shared<Mooney_Rivlin_hyperelastic_law>(true,true)));
2198  ("Plane_Strain_Compressible_Neo_Hookean_sigma",
2199  std::make_shared<AHL_wrapper_sigma>
2200  (std::make_shared<plane_strain_hyperelastic_law>(cneolaw)));
2202  ("Plane_Strain_Compressible_Neo_Hookean_PK2",
2203  std::make_shared<AHL_wrapper_sigma>
2204  (std::make_shared<plane_strain_hyperelastic_law>(cneolaw)));
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);
2213  ("Compressible_Neo_Hookean_Bonet_sigma",
2214  std::make_shared<AHL_wrapper_sigma>(cneobolaw));
2216  ("Compressible_Neo_Hookean_Bonet_PK2",
2217  std::make_shared<AHL_wrapper_sigma>(cneobolaw));
2219  ("Compressible_Neo_Hookean_Bonet_potential",
2220  std::make_shared<AHL_wrapper_potential>
2221  (std::make_shared<Neo_Hookean_hyperelastic_law>(true)));
2223  ("Plane_Strain_Compressible_Neo_Hookean_Bonet_sigma",
2224  std::make_shared<AHL_wrapper_sigma>
2225  (std::make_shared<plane_strain_hyperelastic_law>(cneobolaw)));
2227  ("Plane_Strain_Compressible_Neo_Hookean_Bonet_PK2",
2228  std::make_shared<AHL_wrapper_sigma>
2229  (std::make_shared<plane_strain_hyperelastic_law>(cneobolaw)));
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);
2238  ("Compressible_Neo_Hookean_Ciarlet_sigma",
2239  std::make_shared<AHL_wrapper_sigma>(cneocilaw));
2241  ("Compressible_Neo_Hookean_Ciarlet_PK2",
2242  std::make_shared<AHL_wrapper_sigma>(cneocilaw));
2244  ("Compressible_Neo_Hookean_Ciarlet_potential",
2245  std::make_shared<AHL_wrapper_potential>
2246  (std::make_shared<Neo_Hookean_hyperelastic_law>(false)));
2248  ("Plane_Strain_Compressible_Neo_Hookean_Ciarlet_sigma",
2249  std::make_shared<AHL_wrapper_sigma>
2250  (std::make_shared<plane_strain_hyperelastic_law>(cneocilaw)));
2252  ("Plane_Strain_Compressible_Neo_Hookean_Ciarlet_PK2",
2253  std::make_shared<AHL_wrapper_sigma>
2254  (std::make_shared<plane_strain_hyperelastic_law>(cneocilaw)));
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) {
2270
2271  for (size_type i = 0; i < lawname.size(); ++i)
2273
2274  if (adapted_lawname.compare("SaintVenant_Kirchhoff") == 0) {
2276  } else if (adapted_lawname.compare("Saint_Venant_Kirchhoff") == 0) {
2277
2278  } else if (adapted_lawname.compare("Generalized_Blatz_Ko") == 0) {
2280  } else if (adapted_lawname.compare("Ciarlet_Geymonat") == 0) {
2282  } else if (adapted_lawname.compare("Incompressible_Mooney_Rivlin") == 0) {
2284  } else if (adapted_lawname.compare("Compressible_Mooney_Rivlin") == 0) {
2286  } else if (adapted_lawname.compare("Incompressible_Neo_Hookean") == 0) {
2288  } else if (adapted_lawname.compare("Compressible_Neo_Hookean") == 0 ||
2290  adapted_lawname.compare("Compressible_Neo_Hookean_Ciarlet") == 0 ) {
2292  } else
2293  GMM_ASSERT1(false, lawname << " is not a known hyperelastic law");
2294
2296  }
2297
2298
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);
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
2316
2319
2321  (md, mim, expr, region, true, false,
2322  "Finite strain elasticity brick for " + adapted_lawname + " law");
2323  }
2324
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 + ")"
2335  + test_varname + "))" ;
2337  (md, mim, expr, region, true, false,
2338  "Finite strain incompressibility brick");
2339  }
2340
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) {
2346
2348
2349  std::string expr = "sqrt(3/2)*Norm(Deviator(Cauchy_stress_from_PK2("
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()
Definition: dal_singleton.h:165
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
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
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
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
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
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
Return a reference to the underlying mesh.
Definition: getfem_mesh_fem.h:279
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