GetFEM  5.4.2
getfem_nonlinear_elasticity.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2000-2020 Yves Renard, Julien Pommier
5 
6  This file is a part of GetFEM
7 
8  GetFEM is free software; you can redistribute it and/or modify it
9  under the terms of the GNU Lesser General Public License as published
10  by the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version along with the GCC Runtime Library
12  Exception either version 3.1 or (at your option) any later version.
13  This program is distributed in the hope that it will be useful, but
14  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  License and GCC Runtime Library Exception for more details.
17  You should have received a copy of the GNU Lesser General Public License
18  along with this program; if not, write to the Free Software Foundation,
19  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
20 
21  As a special exception, you may use this file as it is a part of a free
22  software library without restriction. Specifically, if other files
23  instantiate templates or use macros or inline functions from this file,
24  or you compile this file and link it with other files to produce an
25  executable, this file does not by itself cause the resulting executable
26  to be covered by the GNU Lesser General Public License. This exception
27  does not however invalidate any other reasons why the executable file
28  might be covered by the GNU Lesser General Public License.
29 
30 ===========================================================================*/
31 
32 /**@file getfem_nonlinear_elasticity.h
33  @author Yves Renard <Yves.Renard@insa-lyon.fr>,
34  @author Julien Pommier <Julien.Pommier@insa-toulouse.fr>
35  @date July 6, 2004.
36  @brief Non-linear elasticty and incompressibility bricks.
37 */
38 #ifndef GETFEM_NONLINEAR_ELASTICITY_H__
39 #define GETFEM_NONLINEAR_ELASTICITY_H__
40 
41 #include "getfem_models.h"
42 #include "getfem_assembling_tensors.h"
43 #include "getfem_derivatives.h"
44 #include "getfem_interpolation.h"
45 #include "getfem_generic_assembly.h"
46 #include "gmm/gmm_inoutput.h"
47 
48 namespace getfem {
49 
50 
51  int check_symmetry(const base_tensor &t);
52 
53  class abstract_hyperelastic_law;
54  typedef std::shared_ptr<const abstract_hyperelastic_law> phyperelastic_law;
55 
56  /** Base class for material law.
57  Inherit from this class to define a new law.
58  */
59  class abstract_hyperelastic_law {
60  public:
61  mutable int uvflag;
62  size_type nb_params_;
63  phyperelastic_law pl; /* optional reference */
64  void reset_unvalid_flag() const { uvflag = 0; }
65  void inc_unvalid_flag() const { uvflag++; }
66  int get_unvalid_flag() const { return uvflag; }
67  std::string adapted_tangent_term_assembly_fem_data; // should be filled
68  std::string adapted_tangent_term_assembly_cte_data; // to replace the
69  // default assembly
70 
71  virtual scalar_type strain_energy(const base_matrix &E,
72  const base_vector &params,
73  scalar_type det_trans) const = 0;
74  /**True Cauchy stress (for Updated Lagrangian formulation)*/
75  virtual void cauchy_updated_lagrangian(const base_matrix& F,
76  const base_matrix &E,
77  base_matrix &cauchy_stress,
78  const base_vector &params,
79  scalar_type det_trans) const;
80  virtual void sigma(const base_matrix &E, base_matrix &result,
81  const base_vector &params,
82  scalar_type det_trans) const = 0;
83  // the result of grad_sigma has to be completely symmetric.
84  virtual void grad_sigma(const base_matrix &E, base_tensor &result,
85  const base_vector &params, scalar_type det_trans) const = 0;
86 
87  /**cauchy-truesdel tangent moduli, used in updated lagrangian*/
88  virtual void grad_sigma_updated_lagrangian(const base_matrix& F,
89  const base_matrix& E,
90  const base_vector &params,
91  scalar_type det_trans,
92  base_tensor &grad_sigma_ul)const;
93 
94  size_type nb_params() const { return nb_params_; }
95  abstract_hyperelastic_law() { nb_params_ = 0; }
96  virtual ~abstract_hyperelastic_law() {}
97  static void random_E(base_matrix &E);
98  void test_derivatives(size_type N, scalar_type h,
99  const base_vector& param) const;
100  };
101 
102  /** Saint-Venant Kirchhoff hyperelastic law.
103 
104  This is the linear law used in linear elasticity, it is not well
105  suited to large strain. (the convexes may become flat)
106  */
107  struct SaintVenant_Kirchhoff_hyperelastic_law :
108  public abstract_hyperelastic_law {
109  /* W = lambda*0.5*trace(E)^2 + mu*tr(E^2) */
110  virtual scalar_type strain_energy(const base_matrix &E,
111  const base_vector &params,
112  scalar_type det_trans) const;
113  /* sigma = lambda*trace(E) + 2 mu * E */
114  virtual void sigma(const base_matrix &E, base_matrix &result,
115  const base_vector &params, scalar_type det_trans) const;
116  virtual void grad_sigma(const base_matrix &E, base_tensor &result,
117  const base_vector &params, scalar_type det_trans) const;
118  virtual void grad_sigma_updated_lagrangian(const base_matrix& F,
119  const base_matrix& E,
120  const base_vector &params,
121  scalar_type det_trans,
122  base_tensor &grad_sigma_ul)const;
123  SaintVenant_Kirchhoff_hyperelastic_law();
124  };
125 
126 
127  /** Linear law for a membrane (plane stress), orthotropic material
128  caracterized by Ex, Ey, vYX and G, with optional orthotropic prestresses.
129  due to Jean-Yves Heddebaut <jyhed@alpino.be>
130  */
131  struct membrane_elastic_law :
132  public abstract_hyperelastic_law {
133  virtual scalar_type strain_energy(const base_matrix &E,
134  const base_vector &params,
135  scalar_type det_trans) const;
136  virtual void sigma(const base_matrix &E, base_matrix &result,
137  const base_vector &params, scalar_type det_trans) const;
138  virtual void grad_sigma(const base_matrix &E, base_tensor &result,
139  const base_vector &params, scalar_type det_trans) const;
140  membrane_elastic_law() { nb_params_ = 6; }
141  };
142 
143 
144  /** Mooney-Rivlin hyperelastic law
145 
146  To be used for compressible and incompressible problems.
147  Following combinations are possible:
148  not compressible, not neohookean (default): 2 parameters (C1,C2)
149  not compressible, neohookean: 1 parameter (C1)
150  compressible, not neohookean: 3 parameters (C1,C2,D1)
151  compressible, neohookean: 2 parameters (C1,D1)
152  */
153  struct Mooney_Rivlin_hyperelastic_law : public abstract_hyperelastic_law {
154  const bool compressible, neohookean;
155  virtual scalar_type strain_energy(const base_matrix &E,
156  const base_vector &params, scalar_type det_trans) const;
157  virtual void sigma(const base_matrix &E, base_matrix &result,
158  const base_vector &params, scalar_type det_trans) const;
159  virtual void grad_sigma(const base_matrix &E, base_tensor &result,
160  const base_vector &params, scalar_type det_trans) const;
161  explicit Mooney_Rivlin_hyperelastic_law(bool compressible_=false,
162  bool neohookean_=false);
163  };
164 
165  /** Neo-Hookean hyperelastic law variants
166 
167  To be used for compressible problems.
168  For incompressible ones Mooney-Rivlin hyperelastic law can be used.
169  */
170  struct Neo_Hookean_hyperelastic_law : public abstract_hyperelastic_law {
171  const bool bonet;
172  virtual scalar_type strain_energy(const base_matrix &E,
173  const base_vector &params, scalar_type det_trans) const;
174  virtual void sigma(const base_matrix &E, base_matrix &result,
175  const base_vector &params, scalar_type det_trans) const;
176  virtual void grad_sigma(const base_matrix &E, base_tensor &result,
177  const base_vector &params, scalar_type det_trans) const;
178  explicit Neo_Hookean_hyperelastic_law(bool bonet_=true);
179  };
180 
181  /** Blatz_Ko hyperelastic law
182 
183  To be used for compressible problems.
184  */
185  struct generalized_Blatz_Ko_hyperelastic_law : public abstract_hyperelastic_law {
186  virtual scalar_type strain_energy(const base_matrix &E,
187  const base_vector &params, scalar_type det_trans) const;
188  virtual void sigma(const base_matrix &E, base_matrix &result,
189  const base_vector &params, scalar_type det_trans) const;
190  virtual void grad_sigma(const base_matrix &E, base_tensor &result,
191  const base_vector &params, scalar_type det_trans) const;
192  generalized_Blatz_Ko_hyperelastic_law();
193  };
194 
195 
196  /** Ciarlet-Geymonat hyperelastic law
197  */
198  struct Ciarlet_Geymonat_hyperelastic_law : public abstract_hyperelastic_law {
199  // parameters are lambda=params[0], mu=params[1], a=params[2]
200  // The parameter a has to verify a in ]0, mu/2[
201  virtual scalar_type strain_energy(const base_matrix &E,
202  const base_vector &params, scalar_type det_trans) const;
203  virtual void sigma(const base_matrix &E, base_matrix &result,
204  const base_vector &params, scalar_type det_trans) const;
205  virtual void grad_sigma(const base_matrix &E, base_tensor &result,
206  const base_vector &params, scalar_type det_trans) const;
207  Ciarlet_Geymonat_hyperelastic_law() { nb_params_ = 3; }
208  };
209 
210  /** Plane strain hyperelastic law (takes another law as a parameter)
211  */
212  struct plane_strain_hyperelastic_law : public abstract_hyperelastic_law {
213  virtual scalar_type strain_energy(const base_matrix &E,
214  const base_vector &params, scalar_type det_trans) const;
215  virtual void sigma(const base_matrix &E, base_matrix &result,
216  const base_vector &params, scalar_type det_trans) const;
217  virtual void grad_sigma(const base_matrix &E, base_tensor &result,
218  const base_vector &params, scalar_type det_trans) const;
219  plane_strain_hyperelastic_law(const phyperelastic_law &pl_)
220  { pl = pl_; nb_params_ = pl->nb_params(); }
221  };
222 
223 
224 
225 
226  template<typename VECT1, typename VECT2> class elasticity_nonlinear_term
227  : public getfem::nonlinear_elem_term {
228  const mesh_fem &mf;
229  std::vector<scalar_type> U;
230  const mesh_fem *mf_data;
231  const VECT2 &PARAMS;
232  size_type N;
233  size_type NFem;
234  const abstract_hyperelastic_law &AHL;
235  base_vector params, coeff;
236  base_matrix E, Sigma, gradU;
237  base_tensor tt;
238  bgeot::multi_index sizes_;
239  int version;
240  public:
241  elasticity_nonlinear_term(const mesh_fem &mf_, const VECT1 &U_,
242  const mesh_fem *mf_data_, const VECT2 &PARAMS_,
243  const abstract_hyperelastic_law &AHL_,
244  int version_)
245  : mf(mf_), U(mf_.nb_basic_dof()), mf_data(mf_data_), PARAMS(PARAMS_),
246  N(mf_.linked_mesh().dim()), NFem(mf_.get_qdim()), AHL(AHL_),
247  params(AHL_.nb_params()), E(N, N), Sigma(N, N), gradU(NFem, N),
248  tt(N, N, N, N), sizes_(NFem, N, NFem, N),
249  version(version_) {
250  switch (version) {
251  case 0 : break; // tangent term
252  case 1 : sizes_.resize(2); break; // rhs
253  case 2 : sizes_.resize(1); sizes_[0] = 1; break; // strain energy
254  case 3 : sizes_.resize(2); break; // Id + grad(u)
255  }
256 
257  mf.extend_vector(U_, U);
258  if (gmm::vect_size(PARAMS) == AHL_.nb_params())
259  gmm::copy(PARAMS, params);
260  }
261 
262  const bgeot::multi_index &sizes(size_type) const { return sizes_; }
263 
264  virtual void compute(getfem::fem_interpolation_context& ctx,
265  bgeot::base_tensor &t) {
266  size_type cv = ctx.convex_num();
267  slice_vector_on_basic_dof_of_element(mf, U, cv, coeff);
268  ctx.pf()->interpolation_grad(ctx, coeff, gradU, mf.get_qdim());
269 
270  for (unsigned int alpha = 0; alpha < N; ++alpha)
271  gradU(alpha, alpha) += scalar_type(1);
272 
273  if (version == 3) {
274  for (size_type n = 0; n < NFem; ++n)
275  for (size_type m = 0; m < N; ++m)
276  t(n,m) = gradU(n,m);
277  return;
278  }
279 
280  gmm::mult(gmm::transposed(gradU), gradU, E);
281  for (unsigned int alpha = 0; alpha < N; ++alpha)
282  E(alpha, alpha) -= scalar_type(1);
283  gmm::scale(E, scalar_type(0.5));
284 
285  scalar_type det_trans = gmm::lu_det(gradU);
286 
287  if (version == 2) {
288  t[0] = AHL.strain_energy(E, params, det_trans);
289  return;
290  }
291 
292  AHL.sigma(E, Sigma, params, det_trans);
293 
294  if (version == 0) {
295  AHL.grad_sigma(E, tt, params, det_trans);
296  for (size_type n = 0; n < NFem; ++n)
297  for (size_type m = 0; m < N; ++m)
298  for (size_type l = 0; l < N; ++l)
299  for (size_type k = 0; k < NFem; ++k) {
300  scalar_type aux = (k == n) ? Sigma(m,l) : 0.0;
301  for (size_type j = 0; j < N; ++j)
302  for (size_type i = 0; i < N; ++i) {
303  aux += gradU(n ,j) * gradU(k, i) * tt(j, m, i, l);
304  }
305  t(n, m, k, l) = aux;
306  }
307  } else {
308  if (det_trans < scalar_type(0)) AHL.inc_unvalid_flag();
309  for (size_type i = 0; i < NFem; ++i)
310  for (size_type j = 0; j < N; ++j) {
311  scalar_type aux(0);
312  for (size_type k = 0; k < N; ++k)
313  aux += gradU(i, k) * Sigma(k, j);
314  t(i,j) = aux;
315  }
316  }
317  }
318 
319  virtual void prepare(fem_interpolation_context& ctx, size_type ) {
320  if (mf_data) {
321  size_type cv = ctx.convex_num();
322  size_type nb = AHL.nb_params();
323  coeff.resize(mf_data->nb_basic_dof_of_element(cv)*nb);
324  for (size_type i = 0; i < mf_data->nb_basic_dof_of_element(cv); ++i)
325  for (size_type k = 0; k < nb; ++k)
326  coeff[i * nb + k]
327  = PARAMS[mf_data->ind_basic_dof_of_element(cv)[i]*nb+k];
328  ctx.pf()->interpolation(ctx, coeff, params, dim_type(nb));
329  }
330  }
331 
332  };
333 
334 
335  /**
336  Tangent matrix for the non-linear elasticity
337  @ingroup asm
338  */
339  template<typename MAT, typename VECT1, typename VECT2>
340  void asm_nonlinear_elasticity_tangent_matrix
341  (const MAT &K_, const mesh_im &mim, const getfem::mesh_fem &mf,
342  const VECT1 &U, const getfem::mesh_fem *mf_data, const VECT2 &PARAMS,
343  const abstract_hyperelastic_law &AHL,
344  const mesh_region &rg = mesh_region::all_convexes()) {
345  MAT &K = const_cast<MAT &>(K_);
346  GMM_ASSERT1(mf.get_qdim() >= mf.linked_mesh().dim(),
347  "wrong qdim for the mesh_fem");
348 
349  elasticity_nonlinear_term<VECT1, VECT2>
350  nterm(mf, U, mf_data, PARAMS, AHL, 0);
351  elasticity_nonlinear_term<VECT1, VECT2>
352  nterm2(mf, U, mf_data, PARAMS, AHL, 3);
353 
354  getfem::generic_assembly assem;
355  if (mf_data)
356  if (AHL.adapted_tangent_term_assembly_fem_data.size() > 0)
357  assem.set(AHL.adapted_tangent_term_assembly_cte_data);
358  else
359  assem.set("M(#1,#1)+=sym(comp(NonLin$1(#1,#2)(i,j,k,l).vGrad(#1)(:,i,j).vGrad(#1)(:,k,l)))");
360  else
361  if (AHL.adapted_tangent_term_assembly_cte_data.size() > 0)
362  assem.set(AHL.adapted_tangent_term_assembly_cte_data);
363  else
364  assem.set("M(#1,#1)+=sym(comp(NonLin$1(#1)(i,j,k,l).vGrad(#1)(:,i,j).vGrad(#1)(:,k,l)))");
365  assem.push_mi(mim);
366  assem.push_mf(mf);
367  if (mf_data) assem.push_mf(*mf_data);
368  assem.push_data(PARAMS);
369  assem.push_nonlinear_term(&nterm);
370  assem.push_nonlinear_term(&nterm2);
371  assem.push_mat(K);
372  assem.assembly(rg);
373  }
374 
375 
376  /**
377  @ingroup asm
378  */
379  template<typename VECT1, typename VECT2, typename VECT3>
380  void asm_nonlinear_elasticity_rhs
381  (const VECT1 &R_, const mesh_im &mim, const getfem::mesh_fem &mf,
382  const VECT2 &U, const getfem::mesh_fem *mf_data, const VECT3 &PARAMS,
383  const abstract_hyperelastic_law &AHL,
384  const mesh_region &rg = mesh_region::all_convexes()) {
385  VECT1 &R = const_cast<VECT1 &>(R_);
386  GMM_ASSERT1(mf.get_qdim() >= mf.linked_mesh().dim(),
387  "wrong qdim for the mesh_fem");
388 
389  elasticity_nonlinear_term<VECT2, VECT3>
390  nterm(mf, U, mf_data, PARAMS, AHL, 1);
391 
392  getfem::generic_assembly assem;
393  if (mf_data)
394  assem.set("t=comp(NonLin(#1,#2).vGrad(#1)); V(#1) += t(i,j,:,i,j)");
395  else
396  assem.set("t=comp(NonLin(#1).vGrad(#1)); V(#1) += t(i,j,:,i,j)");
397  // comp() to be optimized ?
398  assem.push_mi(mim);
399  assem.push_mf(mf);
400  if (mf_data) assem.push_mf(*mf_data);
401  assem.push_nonlinear_term(&nterm);
402  assem.push_vec(R);
403  assem.assembly(rg);
404  }
405 
406  // added by Jean-Yves Heddebaut <jyhed@alpino.be>
407  int levi_civita(int i,int j,int k);
408 
409 
410  /**@ingroup asm
411  */
412  template<typename VECT2, typename VECT3>
413  scalar_type asm_elastic_strain_energy
414  (const mesh_im &mim, const getfem::mesh_fem &mf,
415  const VECT2 &U, const getfem::mesh_fem *mf_data, const VECT3 &PARAMS,
416  const abstract_hyperelastic_law &AHL,
417  const mesh_region &rg = mesh_region::all_convexes()) {
418 
419  GMM_ASSERT1(mf.get_qdim() >= mf.linked_mesh().dim(),
420  "wrong qdim for the mesh_fem");
421 
422  elasticity_nonlinear_term<VECT2, VECT3>
423  nterm(mf, U, mf_data, PARAMS, AHL, 2);
424  std::vector<scalar_type> V(1);
425 
426  getfem::generic_assembly assem;
427  if (mf_data)
428  assem.set("V() += comp(NonLin(#1,#2))(i)");
429  else
430  assem.set("V() += comp(NonLin(#1))(i)");
431 
432  assem.push_mi(mim);
433  assem.push_mf(mf);
434  if (mf_data) assem.push_mf(*mf_data);
435  assem.push_nonlinear_term(&nterm);
436  assem.push_vec(V);
437  assem.assembly(rg);
438 
439  return V[0];
440  }
441 
442 
443 
444 
445 
446 
447 
448 
449 
450  /* ******************************************************************** */
451  /* Mixed nonlinear incompressibility assembly procedure */
452  /* ******************************************************************** */
453 
454  template<typename VECT1> class incomp_nonlinear_term
455  : public getfem::nonlinear_elem_term {
456 
457  const mesh_fem &mf;
458  std::vector<scalar_type> U;
459  size_type N;
460  base_vector coeff;
461  base_matrix gradPhi;
462  bgeot::multi_index sizes_;
463  int version;
464 
465  public:
466  incomp_nonlinear_term(const mesh_fem &mf_, const VECT1 &U_,
467  int version_)
468  : mf(mf_), U(mf_.nb_basic_dof()),
469  N(mf_.get_qdim()),
470  gradPhi(N, N), sizes_(N, N),
471  version(version_) {
472  if (version == 1) { sizes_.resize(1); sizes_[0] = 1; }
473  mf.extend_vector(U_, U);
474  }
475 
476  const bgeot::multi_index &sizes(size_type) const { return sizes_; }
477 
478  virtual void compute(getfem::fem_interpolation_context& ctx,
479  bgeot::base_tensor &t) {
480  size_type cv = ctx.convex_num();
481  slice_vector_on_basic_dof_of_element(mf, U, cv, coeff);
482  ctx.pf()->interpolation_grad(ctx, coeff, gradPhi, mf.get_qdim());
483  gmm::add(gmm::identity_matrix(), gradPhi);
484  scalar_type det = gmm::lu_inverse(gradPhi);
485 
486  if (version != 1) {
487  if (version == 2) det = sqrt(gmm::abs(det));
488  for (size_type i = 0; i < N; ++i)
489  for (size_type j = 0; j < N; ++j) {
490  t(i,j) = - det * gradPhi(j,i);
491  }
492  }
493  else t[0] = scalar_type(1) - det;
494 
495  }
496  };
497 
498  /**@ingroup asm*/
499  template<typename MAT1, typename MAT2, typename VECT1, typename VECT2>
500  void asm_nonlinear_incomp_tangent_matrix(const MAT1 &K_, const MAT2 &B_,
501  const mesh_im &mim,
502  const mesh_fem &mf_u,
503  const mesh_fem &mf_p,
504  const VECT1 &U, const VECT2 &P,
505  const mesh_region &rg = mesh_region::all_convexes()) {
506  MAT1 &K = const_cast<MAT1 &>(K_);
507  MAT2 &B = const_cast<MAT2 &>(B_);
508  GMM_ASSERT1(mf_u.get_qdim() == mf_u.linked_mesh().dim(),
509  "wrong qdim for the mesh_fem");
510 
511  incomp_nonlinear_term<VECT1> ntermk(mf_u, U, 0);
512  incomp_nonlinear_term<VECT1> ntermb(mf_u, U, 2);
513  getfem::generic_assembly
514  assem("P=data(#2);"
515  "t=comp(NonLin$1(#1).vGrad(#1).Base(#2));"
516  "M$2(#1,#2)+= t(i,j,:,i,j,:);"
517  /*"w=comp(NonLin$2(#1).vGrad(#1).NonLin$2(#1).vGrad(#1).Base(#2));"
518  "M$1(#1,#1)+= w(j,i,:,j,k, m,k,:,m,i,p).P(p)"
519  "-w(i,j,:,i,j, k,l,:,k,l,p).P(p)"*/
520  /*
521  "w=comp(vGrad(#1).NonLin$2(#1).vGrad(#1).NonLin$2(#1).Base(#2));"
522  "M$1(#1,#1)+= w(:,j,k, j,i, :,m,i, m,k, p).P(p)"
523  "-w(:,j,i, j,i, :,m,l, m,l, p).P(p)"
524  */
525  "w1=comp(vGrad(#1)(:,j,k).NonLin$2(#1)(j,i).vGrad(#1)(:,m,i).NonLin$2(#1)(m,k).Base(#2)(p).P(p));"
526  "w2=comp(vGrad(#1)(:,j,i).NonLin$2(#1)(j,i).vGrad(#1)(:,m,l).NonLin$2(#1)(m,l).Base(#2)(p).P(p));"
527  "M$1(#1,#1)+= w1-w2"
528  );
529 
530  assem.push_mi(mim);
531  assem.push_mf(mf_u);
532  assem.push_mf(mf_p);
533  assem.push_nonlinear_term(&ntermk);
534  assem.push_nonlinear_term(&ntermb);
535  assem.push_mat(K);
536  assem.push_mat(B);
537  assem.push_data(P);
538  assem.assembly(rg);
539  }
540 
541 
542  /**@ingroup asm
543  */
544  template<typename VECT1, typename VECT2, typename VECT3>
545  void asm_nonlinear_incomp_rhs
546  (const VECT1 &R_U_, const VECT1 &R_P_, const mesh_im &mim,
547  const getfem::mesh_fem &mf_u, const getfem::mesh_fem &mf_p,
548  const VECT2 &U, const VECT3 &P,
549  const mesh_region &rg = mesh_region::all_convexes()) {
550  VECT1 &R_U = const_cast<VECT1 &>(R_U_);
551  VECT1 &R_P = const_cast<VECT1 &>(R_P_);
552  GMM_ASSERT1(mf_u.get_qdim() == mf_u.linked_mesh().dim(),
553  "wrong qdim for the mesh_fem");
554 
555  incomp_nonlinear_term<VECT2> nterm_tg(mf_u, U, 0);
556  incomp_nonlinear_term<VECT2> nterm(mf_u, U, 1);
557 
558  getfem::generic_assembly
559  assem("P=data(#2); "
560  "t=comp(NonLin$1(#1).vGrad(#1).Base(#2));"
561  "V$1(#1) += t(i,j,:,i,j,k).P(k);"
562  "w=comp(NonLin$2(#1).Base(#2)); V$2(#2) += w(1,:)");
563  // assem() to be optimized ?
564 
565  assem.push_mi(mim);
566  assem.push_mf(mf_u);
567  assem.push_mf(mf_p);
568  assem.push_nonlinear_term(&nterm_tg);
569  assem.push_nonlinear_term(&nterm);
570  assem.push_vec(R_U);
571  assem.push_vec(R_P);
572  assem.push_data(P);
573  assem.assembly(rg);
574  }
575 
576 
577 
578  //===========================================================================
579  //
580  // Bricks
581  //
582  //===========================================================================
583 
584 
585  /** Add a nonlinear (large strain) elasticity term to the model with
586  respect to the variable
587  `varname` (deprecated brick, use add_finite_strain_elaticity instead).
588  Note that the constitutive law is described by `AHL` which
589  should not be freed while the model is used. `dataname` described the
590  parameters of the constitutive laws. It could be a vector of value
591  of length the number of parameter of the constitutive law or a vector
592  field described on a finite element method.
593  */
594  size_type add_nonlinear_elasticity_brick
595  (model &md, const mesh_im &mim, const std::string &varname,
596  const phyperelastic_law &AHL, const std::string &dataname,
597  size_type region = size_type(-1));
598 
599 
600 
601  void compute_Von_Mises_or_Tresca
602  (model &md, const std::string &varname, const phyperelastic_law &AHL,
603  const std::string &dataname, const mesh_fem &mf_vm,
604  model_real_plain_vector &VM, bool tresca);
605 
606 
607  void compute_sigmahathat(model &md,
608  const std::string &varname,
609  const phyperelastic_law &AHL,
610  const std::string &dataname,
611  const mesh_fem &mf_sigma,
612  model_real_plain_vector &SIGMA);
613 
614 
615  /**
616  Compute the Von-Mises stress or the Tresca stress of a field
617  with respect to the constitutive elasticity law AHL (only valid in 3D).
618  */
619  template <class VECTVM> void compute_Von_Mises_or_Tresca
620  (model &md, const std::string &varname, const phyperelastic_law &AHL,
621  const std::string &dataname, const mesh_fem &mf_vm,
622  VECTVM &VM, bool tresca) {
623  model_real_plain_vector VMM(mf_vm.nb_dof());
624  compute_Von_Mises_or_Tresca
625  (md, varname, AHL, dataname, mf_vm, VMM, tresca);
626  gmm::copy(VMM, VM);
627  }
628 
629  /** Add a nonlinear incompressibility term (for large strain elasticity)
630  to the model with respect to the variable
631  `varname` (the displacement) and `multname` (the pressure).
632  */
633  size_type add_nonlinear_incompressibility_brick
634  (model &md, const mesh_im &mim, const std::string &varname,
635  const std::string &multname, size_type region = size_type(-1));
636 
637  //===========================================================================
638  // High-level generic assembly bricks
639  //===========================================================================
640 
641  /** Add a finite strain elasticity brick
642  to the model with respect to the variable
643  `varname` (the displacement).
644  For 2D meshes, switch automatically to plane strain elasticity.
645  High-level generic assembly version.
646  */
647  size_type add_finite_strain_elasticity_brick
648  (model &md, const mesh_im &mim, const std::string &lawname,
649  const std::string &varname, const std::string &params,
650  size_type region = size_type(-1));
651 
652 
653  /** Add a finite strain incompressibility term (for large strain elasticity)
654  to the model with respect to the variable
655  `varname` (the displacement) and `multname` (the pressure).
656  High-level generic assembly version.
657  */
658  size_type add_finite_strain_incompressibility_brick
659  (model &md, const mesh_im &mim, const std::string &varname,
660  const std::string &multname, size_type region = size_type(-1));
661 
662  /**
663  Interpolate the Von-Mises stress of a field `varname`
664  with respect to the nonlinear elasticity constitutive law `lawname`
665  with parameters `params` (only valid in 3D).
666  */
667  void compute_finite_strain_elasticity_Von_Mises
668  (model &md, const std::string &lawname, const std::string &varname,
669  const std::string &params, const mesh_fem &mf_vm,
670  model_real_plain_vector &VM,
671  const mesh_region &rg=mesh_region::all_convexes());
672 
673  IS_DEPRECATED inline void finite_strain_elasticity_Von_Mises
674  (model &md, const std::string &varname, const std::string &lawname,
675  const std::string &params, const mesh_fem &mf_vm,
676  model_real_plain_vector &VM,
677  const mesh_region &rg=mesh_region::all_convexes()) {
678  compute_finite_strain_elasticity_Von_Mises(md, varname, lawname, params,
679  mf_vm, VM, rg);
680  }
681 
682 } /* end of namespace getfem. */
683 
684 
685 #endif /* GETFEM_NONLINEAR_ELASTICITY_H__ */
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::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::membrane_elastic_law
Linear law for a membrane (plane stress), orthotropic material caracterized by Ex,...
Definition: getfem_nonlinear_elasticity.h:131
getfem::Neo_Hookean_hyperelastic_law
Neo-Hookean hyperelastic law variants.
Definition: getfem_nonlinear_elasticity.h:170
getfem_assembling_tensors.h
Generic assembly implementation.
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
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
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::mesh_im
Describe an integration method linked to a mesh.
Definition: getfem_mesh_im.h:47
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::abstract_hyperelastic_law
Base class for material law.
Definition: getfem_nonlinear_elasticity.h:59
getfem::generic_assembly::push_mf
void push_mf(const mesh_fem &mf_)
Add a new mesh_fem.
Definition: getfem_assembling_tensors.h:650
getfem::SaintVenant_Kirchhoff_hyperelastic_law
Saint-Venant Kirchhoff hyperelastic law.
Definition: getfem_nonlinear_elasticity.h:107
getfem::generic_assembly::push_mat
void push_mat(const MAT &m)
Add a new output matrix (fake const version..)
Definition: getfem_assembling_tensors.h:671
getfem::mesh_fem
Describe a finite element method linked to a mesh.
Definition: getfem_mesh_fem.h:148
getfem::fem_interpolation_context::convex_num
size_type convex_num() const
get the current convex number
Definition: getfem_fem.cc:52
getfem::mesh_region::all_convexes
static mesh_region all_convexes()
provide a default value for the mesh_region parameters of assembly procedures etc.
Definition: getfem_mesh_region.h:142
getfem::model
`‘Model’' variables store the variables, the data and the description of a model.
Definition: getfem_models.h: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_assembly::push_data
void push_data(const VEC &d)
Add a new data (dense array)
Definition: getfem_assembling_tensors.h:659
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
getfem::fem_interpolation_context
structure passed as the argument of fem interpolation functions.
Definition: getfem_fem.h:749
getfem::generic_assembly::assembly
void assembly(const mesh_region &region=mesh_region::all_convexes())
do the assembly on the specified region (boundary or set of convexes)
Definition: getfem_assembling_tensors.cc:1814
getfem::plane_strain_hyperelastic_law
Plane strain hyperelastic law (takes another law as a parameter)
Definition: getfem_nonlinear_elasticity.h:212
getfem_models.h
Model representation in Getfem.
getfem::mesh_fem::nb_basic_dof_of_element
virtual size_type nb_basic_dof_of_element(size_type cv) const
Return the number of degrees of freedom attached to a given convex.
Definition: getfem_mesh_fem.h:494
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.
getfem::generic_assembly
Generic assembly of vectors, matrices.
Definition: getfem_assembling_tensors.h:607
getfem::Ciarlet_Geymonat_hyperelastic_law
Ciarlet-Geymonat hyperelastic law.
Definition: getfem_nonlinear_elasticity.h:198
gmm_inoutput.h
Input/output on sparse matrices.
getfem::Mooney_Rivlin_hyperelastic_law
Mooney-Rivlin hyperelastic law.
Definition: getfem_nonlinear_elasticity.h:153
getfem_derivatives.h
Compute the gradient of a field on a getfem::mesh_fem.
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
getfem::generic_assembly::push_vec
void push_vec(VEC &v)
Add a new output vector.
Definition: getfem_assembling_tensors.h:663
getfem::generic_assembly::push_mi
void push_mi(const mesh_im &im_)
Add a new mesh_im.
Definition: getfem_assembling_tensors.h:652
bgeot::alpha
size_type alpha(short_type n, short_type d)
Return the value of which is the number of monomials of a polynomial of variables and degree .
Definition: bgeot_poly.cc:47
getfem::generalized_Blatz_Ko_hyperelastic_law
Blatz_Ko hyperelastic law.
Definition: getfem_nonlinear_elasticity.h:185
getfem::slice_vector_on_basic_dof_of_element
void slice_vector_on_basic_dof_of_element(const mesh_fem &mf, const VEC1 &vec, size_type cv, VEC2 &coeff, size_type qmult1=size_type(-1), size_type qmult2=size_type(-1))
Given a mesh_fem.
Definition: getfem_mesh_fem.h:659
getfem_interpolation.h
Interpolation of fields from a mesh_fem onto another.
getfem::nonlinear_elem_term
abstract class for integration of non-linear terms into the mat_elem computations the nonlinear term ...
Definition: getfem_mat_elem_type.h:67
getfem::mesh_fem::ind_basic_dof_of_element
virtual ind_dof_ct ind_basic_dof_of_element(size_type cv) const
Give an array of the dof numbers a of convex.
Definition: getfem_mesh_fem.h:450
getfem::mesh_region
structure used to hold a set of convexes and/or convex faces.
Definition: getfem_mesh_region.h:55
getfem::fem_interpolation_context::pf
const pfem pf() const
get the current FEM descriptor
Definition: getfem_fem.h:782
getfem::generic_assembly::push_nonlinear_term
void push_nonlinear_term(pnonlinear_elem_term net)
Add a new non-linear term.
Definition: getfem_assembling_tensors.h:655
getfem::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
getfem::mesh_fem::nb_dof
virtual size_type nb_dof() const
Return the total number of degrees of freedom.
Definition: getfem_mesh_fem.h:562