GetFEM++  5.3
plasticity.cc
Go to the documentation of this file.
1 /*===========================================================================
2 
3  Copyright (C) 2000-2017 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 /**
24  @file plasticity.cc
25  @brief Small deformation plasticity problem.
26 
27  This program is used to check that getfem++ is working.
28  This is also a good example of use of GetFEM++.
29 */
30 
35 #include "getfem/getfem_export.h"
36 
37 using std::endl; using std::cout; using std::cerr;
38 using std::ends; using std::cin;
39 
40 /* some GetFEM++ types that we will be using */
41 
42 /* special class for small (dim<16) vectors */
44 /* geometrical nodes(derived from base_small_vector)*/
45 using bgeot::base_node;
46 using bgeot::scalar_type; /* = double */
47 using bgeot::size_type; /* = unsigned long */
48 using bgeot::base_matrix; /* small dense matrix. */
49 
50 /* definition of some matrix/vector types. These ones are built
51  using the predefined types in Gmm++ */
53 typedef getfem::modeling_standard_sparse_matrix sparse_matrix;
54 typedef getfem::modeling_standard_plain_vector plain_vector;
55 
56 template<typename VEC>
57 static void vecsave(std::string fname, const VEC& V);
58 
59 //function to save a vector
60 template<typename VEC>
61 static void vecsave( std::string fname, const VEC& V) {
62 
63  std::ofstream f(fname.c_str()); f.precision(16);
64  for (size_type i=0; i < V.size(); ++i)
65  f << V[i] << "\n";
66 }
67 
68 
69 
70 //=================================================================
71 // structure for the elastoplastic problem
72 //=================================================================
73 
74 struct elastoplasticity_problem {
75 
76  enum { DIRICHLET_BOUNDARY_NUM = 0,
77  NEUMANN_BOUNDARY_NUM = 1};
78 
79  getfem::mesh mesh; /* the mesh */
80  getfem::mesh_im mim; /* integration methods. */
81  getfem::im_data mim_data; /* Mim data for the pastic strain. */
82  getfem::mesh_fem mf_u; /* main mesh_fem, for the
83  elastoplastic displacement */
84  getfem::mesh_fem mf_xi; /* mesh_fem, for the plastic multiplier. */
85  getfem::mesh_fem mf_rhs; /* mesh_fem for the right hand side
86  (f(x),..) */
87  scalar_type lambda, mu; /* Lamé coefficients. */
88 
89  scalar_type residual; /* max residual for the iterative solvers */
90  scalar_type sigma_y;
91  size_type flag_hyp;
92  std::string datafilename;
93  bgeot::md_param PARAM;
94  bool do_export;
95 
96  bool solve(plain_vector &U);
97  void init(void);
98 
99  elastoplasticity_problem(void) : mim(mesh), mim_data(mim), mf_u(mesh),
100  mf_xi(mesh), mf_rhs(mesh) {}
101 
102 };
103 
104 
105 
106 
107 /* Read parameters from the .param file, build the mesh,
108  set finite element and integration methods
109  and selects the boundaries.
110 */
111 void elastoplasticity_problem::init(void) {
112 
113  std::string MESH_TYPE =
114  PARAM.string_value("MESH_TYPE","Mesh type ");
115  std::string FEM_TYPE =
116  PARAM.string_value("FEM_TYPE","FEM name");
117  std::string FEM_TYPE_XI =
118  PARAM.string_value("FEM_TYPE_XI","FEM name");
119  std::string INTEGRATION =
120  PARAM.string_value("INTEGRATION",
121  "Name of integration method");
122  do_export = (PARAM.int_value("EXPORT", "Perform or not the vtk export") != 0);
123 
124  cout << "MESH_TYPE=" << MESH_TYPE << "\n";
125  cout << "FEM_TYPE=" << FEM_TYPE << "\n";
126  cout << "INTEGRATION=" << INTEGRATION << "\n";
127 
128  residual = PARAM.real_value("RESIDUAL", "residual");
129 
130  // file to save the mesh
131  datafilename = PARAM.string_value("ROOTFILENAME",
132  "Filename for saving");
133 
134  /* First step : build the mesh */
135  size_type N;
136  bgeot::pgeometric_trans pgt = 0;
137 
138  if (MESH_TYPE != "load") {
139  std::cout << "created getfem mesh" << "\n";
140  pgt = bgeot::geometric_trans_descriptor(MESH_TYPE);
141  N = pgt->dim();
142  std::vector<size_type> nsubdiv(N);
143  nsubdiv[0]=PARAM.int_value
144  ("NX", "Number of space steps in x direction ");
145  nsubdiv[1]=PARAM.int_value
146  ("NY", "Number of space steps in y direction ");
147 
148  if(N==3)
149  nsubdiv[2]=PARAM.int_value
150  ("NZ", "Number of space steps in z direction ");
151  getfem::regular_unit_mesh(mesh, nsubdiv, pgt,
152  PARAM.int_value("MESH_NOISED")!= 0);
153 
154  bgeot::base_matrix M(N,N);
155 
156  for (size_type i=0; i < N; ++i) {
157  static const char *t[] = {"LX","LY","LZ"};
158  M(i,i) = (i<3) ? PARAM.real_value(t[i],t[i]) : 1.0;
159  }
160 
161  if (N>1) {
162  M(0,1) = PARAM.real_value("INCLINE") *
163  PARAM.real_value("LY");
164  }
165 
166  /* scale the unit mesh to [LX,LY,..] and incline it */
167  mesh.transformation(M);
168 
169  } else {
170  std ::cout << "mesh from pdetool" << "\n";
171  std::string MESH_FILE
172  = PARAM.string_value("MESH_FILE", "Mesh file name");
173 
174  mesh.read_from_file(MESH_FILE);
175 
176  N = mesh.dim();
177  pgt = mesh.trans_of_convex
178  (mesh.convex_index().first_true());
179  }
180 
181  mu = PARAM.real_value("MU", "Lamé coefficient mu");
182  /* muT = PARAM.real_value("MUT", "Lamé coefficient muT");
183  lambdaB = PARAM.real_value("LAMBDAB",
184  "Lamé coefficient lambdaB");
185  */
186  lambda = PARAM.real_value("LAMBDA",
187  "Lamé coefficient lambda");
188  mf_u.set_qdim(bgeot::dim_type(N));
189 
190 
191  /* set the finite element on the mf_u */
192  getfem::pfem pf_u =
193  getfem::fem_descriptor(FEM_TYPE);
194  getfem::pintegration_method ppi =
195  getfem::int_method_descriptor(INTEGRATION);
196 
197  mim.set_integration_method(mesh.convex_index(), ppi);
198  mim_data.set_tensor_size(bgeot::multi_index(N,N));
199  mf_u.set_finite_element(mesh.convex_index(), pf_u);
200 
201 
202  /* set the finite element on the mf_sigma */
203  getfem::pfem pf_xi = getfem::fem_descriptor(FEM_TYPE_XI);
204  mf_xi.set_finite_element(mesh.convex_index(), pf_xi);
205 
206 
207  /* set the finite element on mf_rhs
208  (same as mf_u is DATA_FEM_TYPE is not used in the .param file)*/
209  std::string data_fem_name
210  = PARAM.string_value("DATA_FEM_TYPE");
211 
212  if (data_fem_name.size() == 0) {
213  GMM_ASSERT1(pf_u->is_lagrange(),
214  "You are using a non-lagrange FEM. "
215  << "In that case you need to set "
216  << "DATA_FEM_TYPE in the .param file");
217  mf_rhs.set_finite_element(mesh.convex_index(), pf_u);
218 
219  } else {
220  mf_rhs.set_finite_element(mesh.convex_index(),
221  getfem::fem_descriptor(data_fem_name));
222  }
223 
224 
225  /* set boundary conditions
226  * (Neuman on the upper face, Dirichlet elsewhere) */
227  cout << "Selecting Neumann and Dirichlet boundaries\n";
228  getfem::mesh_region border_faces;
229  getfem::outer_faces_of_mesh(mesh, border_faces);
230 
231  for (getfem::mr_visitor it(border_faces);
232  !it.finished(); ++it) {
233  assert(it.is_face());
234  base_node un
235  = mesh.normal_of_face_of_convex(it.cv(), it.f());
236  un /= gmm::vect_norm2(un);
237 
238  if (gmm::abs(un[0] - 1.0) < 1.0E-7)
239  mesh.region(NEUMANN_BOUNDARY_NUM).add(it.cv(), it.f());
240  else if (gmm::abs(un[0] + 1.0) < 1.0E-7)
241  mesh.region(DIRICHLET_BOUNDARY_NUM).add(it.cv(),it.f());
242 
243  }
244 
245  // Plasticity part
246  sigma_y = PARAM.real_value("SIGMA_Y", "plasticity yield stress");
247  flag_hyp=PARAM.int_value("FLAG_HYP");
248 }
249 
250 
251 
252 //==================================================================
253 // Model.
254 //==================================================================
255 
256 bool elastoplasticity_problem::solve(plain_vector &U) {
257 
258  size_type nb_dof_rhs = mf_rhs.nb_dof();
259  size_type N = mesh.dim();
260  getfem::model model;
261 
262  gmm::set_traces_level(1);
263 
264  // Main unknown of the problem.
265  model.add_fem_variable("u", mf_u);
266  model.add_fem_data("Previous_u", mf_u);
267 
268  model.add_initialized_scalar_data("lambda", lambda);
269  model.add_initialized_scalar_data("mu", mu);
270  model.add_initialized_scalar_data("sigma_y", sigma_y);
271 
272  model.add_fem_data("xi", mf_xi);
273  model.add_fem_data("Previous_xi", mf_xi);
274  model.add_im_data("Previous_Ep", mim_data);
275 
276  /* choose the projection type */
277  getfem::pconstraints_projection
278  proj = std::make_shared<getfem::VM_projection>(0);
279 
280  std::vector<std::string> plastic_variables = {"u", "xi", "Previous_Ep"};
281  std::vector<std::string> plastic_data = {"lambda", "mu", "sigma_y"};
282 
283 
285  (model, mim, "Prandtl Reuss", getfem::DISPLACEMENT_ONLY,
286  plastic_variables, plastic_data);
287 
288  plain_vector F(nb_dof_rhs * N);
289  model.add_initialized_fem_data("NeumannData", mf_rhs, F);
291  (model, mim, "u", "NeumannData", NEUMANN_BOUNDARY_NUM);
292 
293  model.add_initialized_fem_data("DirichletData", mf_rhs, F);
295  (model, mim, "u", mf_u, DIRICHLET_BOUNDARY_NUM,
296  "DirichletData");
297 
298  const size_type Nb_t = 19;
299  std::vector<scalar_type> T(19);
300  T[0] = 0; T[1] = 0.9032; T[2] = 1; T[3] = 1.1; T[4] = 1.3;
301  T[5] = 1.5; T[6] = 1.7; T[7] = 1.74; T[8] = 1.7; T[9] = 1.5;
302  T[10] = 1.3; T[11] = 1.1; T[12] = 1; T[13] = 0.9032; T[14] = 0.7;
303  T[15] = 0.5; T[16] = 0.3; T[17] = 0.1; T[18] = 0;
304 
305 
306  getfem::mesh_fem mf_vm(mesh);
307  mf_vm.set_classical_discontinuous_finite_element(1);
308  getfem::base_vector VM(mf_vm.nb_dof());
309  getfem::base_vector plast(mf_vm.nb_dof());
310 
311  for (size_type nb = 0; nb < Nb_t; ++nb) {
312  cout << "=============iteration number : " << nb << "==========" << endl;
313 
314  scalar_type t = T[nb];
315 
316  // Defining the Neumann condition right hand side.
317  base_small_vector v(N);
318  v[N-1] = -PARAM.real_value("FORCE");
319  gmm::scale(v,t);
320 
321  for (size_type i = 0; i < nb_dof_rhs; ++i)
322  gmm::copy(v, gmm::sub_vector
323  (F, gmm::sub_interval(i*N, N)));
324 
325  gmm::copy(F, model.set_real_variable("NeumannData"));
326 
327  // Generic solve.
328  cout << "Number of variables : "
329  << model.nb_dof() << endl;
330 
331  getfem::newton_search_with_step_control ls;
332  // getfem::simplest_newton_line_search ls;
333  gmm::iteration iter(residual, 2, 40000);
334  getfem::standard_solve(model, iter,
335  getfem::rselect_linear_solver(model, "superlu"), ls);
336 
338  (model, mim, "Prandtl Reuss", getfem::DISPLACEMENT_ONLY,
339  plastic_variables, plastic_data);
340 
341  // Get the solution and save it
342  gmm::copy(model.real_variable("u"), U);
343 
344 
346  (model, mim, "Prandtl Reuss", getfem::DISPLACEMENT_ONLY,
347  plastic_variables, plastic_data, mf_vm, VM);
348 
349  std::stringstream fname; fname << datafilename << "_" << nb << ".vtk";
350 
351  if (do_export) {
352  getfem::vtk_export exp(fname.str());
353  exp.exporting(mf_vm);
354  exp.write_point_data(mf_vm,VM, "Von Mises stress");
355  exp.write_point_data(mf_u, U, "displacement");
356  }
357 
358  }
359 
360  if (do_export) {
361  cout << "export done, you can view the data file with "
362  "(for example)\n"
363  "mayavi2 -d " << datafilename << "_1.vtk -f "
364  "WarpVector -m Surface -m Outline\n";
365  }
366 
367  return true;
368 }
369 
370 
371 //==================================================================
372 // main program.
373 //==================================================================
374 
375 int main(int argc, char *argv[]) {
376 
377  GMM_SET_EXCEPTION_DEBUG;
378  // Exceptions make a memory fault, to debug.
379 
380  FE_ENABLE_EXCEPT;
381  // Enable floating point exception for Nan.
382 
383  elastoplasticity_problem p;
384  p.PARAM.read_command_line(argc, argv);
385  p.init();
386  p.mesh.write_to_file(p.datafilename + ".mesh");
387  plain_vector U(p.mf_u.nb_dof());
388  if (!p.solve(U)) GMM_ASSERT1(false, "Solve has failed");
389  vecsave(p.datafilename + ".U", U);
390  cout << "Resultats dans fichier : "
391  <<p.datafilename<<".U \n";
392  p.mf_u.write_to_file(p.datafilename + ".meshfem",true);
393  scalar_type t[2]={p.mu,p.lambda};
394  vecsave(p.datafilename+".coef",
395  std::vector<scalar_type>(t, t+2));
396 
397  return 0;
398 }
structure used to hold a set of convexes and/or convex faces.
size_type nb_dof() const
Total number of degrees of freedom in the model.
Export solutions to various formats.
size_type add_small_strain_elastoplasticity_brick(model &md, const mesh_im &mim, std::string lawname, plasticity_unknowns_type unknowns_type, const std::vector< std::string > &varnames, const std::vector< std::string > &params, size_type region=size_type(-1))
Adds a small strain plasticity term to the model md.
pfem fem_descriptor(const std::string &name)
get a fem descriptor from its string name.
Definition: getfem_fem.cc:3950
The Iteration object calculates whether the solution has reached the desired accuracy, or whether the maximum number of iterations has been reached.
Definition: gmm_iter.h:53
size_type APIDECL add_source_term_brick(model &md, const mesh_im &mim, const std::string &varname, const std::string &dataexpr, size_type region=size_type(-1), const std::string &directdataname=std::string())
Add a source term on the variable varname.
Describe a mesh (collection of convexes (elements) and points).
Definition: getfem_mesh.h:95
void small_strain_elastoplasticity_next_iter(model &md, const mesh_im &mim, std::string lawname, plasticity_unknowns_type unknowns_type, const std::vector< std::string > &varnames, const std::vector< std::string > &params, size_type region=size_type(-1))
Function that allows to pass from a time step to another for the small strain plastic brick...
void standard_solve(model &md, gmm::iteration &iter, rmodel_plsolver_type lsolver, abstract_newton_line_search &ls)
A default solver for the model brick system.
pintegration_method int_method_descriptor(std::string name, bool throw_if_not_found=true)
Get an integration method from its name .
Describe an integration method linked to a mesh.
void exporting(const mesh &m)
should be called before write_*_data
Build regular meshes.
pgeometric_trans geometric_trans_descriptor(std::string name)
Get the geometric transformation from its string name.
Plasticty bricks.
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
``Model&#39;&#39; variables store the variables, the data and the description of a model. ...
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
void add_initialized_scalar_data(const std::string &name, T e)
Add a scalar data (i.e.
Standard solvers for model bricks.
"iterator" class for regions.
sparse vector built upon std::vector.
Definition: gmm_def.h:488
std::shared_ptr< const getfem::virtual_fem > pfem
type of pointer on a fem description
Definition: getfem_fem.h:239
Describe a finite element method linked to a mesh.
model_real_plain_vector & set_real_variable(const std::string &name, size_type niter) const
Gives the write access to the vector value of a variable.
void compute_small_strain_elastoplasticity_Von_Mises(model &md, const mesh_im &mim, std::string lawname, plasticity_unknowns_type unknowns_type, const std::vector< std::string > &varnames, const std::vector< std::string > &params, const mesh_fem &mf_vm, model_real_plain_vector &VM, size_type region=size_type(-1))
This function computes the Von Mises stress field with respect to a small strain elastoplasticity ter...
void add_initialized_fem_data(const std::string &name, const mesh_fem &mf, const VECT &v)
Add an initialized fixed size data to the model, assumed to be a vector field if the size of the vect...
void add_fem_data(const std::string &name, const mesh_fem &mf, dim_type qdim=1, size_type niter=1)
Add a data being the dofs of a finite element method to the model.
void add_im_data(const std::string &name, const im_data &im_data, size_type niter=1)
Add data, defined at integration points.
void add_fem_variable(const std::string &name, const mesh_fem &mf, size_type niter=1)
Add a variable being the dofs of a finite element method to the model.
const model_real_plain_vector & real_variable(const std::string &name, size_type niter) const
Gives the access to the vector value of a variable.
VTK export.
Definition: getfem_export.h:67
im_data provides indexing to the integration points of a mesh im object.
Miscelleanous assembly routines for common terms. Use the low-level generic assembly. Prefer the high-level one.
void regular_unit_mesh(mesh &m, std::vector< size_type > nsubdiv, bgeot::pgeometric_trans pgt, bool noised=false)
Build a regular mesh of the unit square/cube/, etc.
std::shared_ptr< const bgeot::geometric_trans > pgeometric_trans
pointer type for a geometric transformation
size_type APIDECL add_Dirichlet_condition_with_multipliers(model &md, const mesh_im &mim, const std::string &varname, const std::string &multname, size_type region, const std::string &dataname=std::string())
Add a Dirichlet condition on the variable varname and the mesh region region.
void APIDECL outer_faces_of_mesh(const mesh &m, const dal::bit_vector &cvlst, convex_face_ct &flist)
returns a list of "exterior" faces of a mesh (i.e.
Definition: getfem_mesh.cc:780