GetFEM  5.4.2
nonlinear_elastostatic.cc
Go to the documentation of this file.
1 /*===========================================================================
2 
3  Copyright (C) 2002-2020 Yves Renard, Julien Pommier.
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  @file nonlinear_elastostatic.cc
24  @brief Nonlinear Elastostatic problem (large strain).
25 
26  A rubber bar is submitted to a large torsion.
27 
28  This program is used to check that getfem++ is working. This is also
29  a good example of use of GetFEM.
30 */
31 
32 #include "getfem/getfem_export.h"
33 #include "getfem/getfem_regular_meshes.h"
34 #include "getfem/getfem_model_solvers.h"
35 #include "getfem/getfem_nonlinear_elasticity.h"
36 #include "getfem/getfem_superlu.h"
37 #include "gmm/gmm.h"
38 
39 using std::endl; using std::cout; using std::cerr;
40 using std::ends; using std::cin;
41 
42 /* some GetFEM types that we will be using */
43 using bgeot::base_small_vector; /* special class for small (dim<16) vectors */
44 using bgeot::base_node; /* geometrical nodes(derived from base_small_vector)*/
45 using bgeot::base_vector;
46 using bgeot::scalar_type; /* = double */
47 using bgeot::size_type; /* = unsigned long */
48 using bgeot::dim_type;
49 using bgeot::base_matrix; /* small dense matrix. */
50 
51 /* definition of some matrix/vector types. These ones are built
52  * using the predefined types in Gmm++
53  */
54 typedef getfem::modeling_standard_sparse_vector sparse_vector;
55 typedef getfem::modeling_standard_sparse_matrix sparse_matrix;
56 typedef getfem::modeling_standard_plain_vector plain_vector;
57 
58 /*
59  structure for the elastostatic problem
60 */
61 struct elastostatic_problem {
62 
63  enum { DIRICHLET_BOUNDARY_NUM = 0, NEUMANN_BOUNDARY_NUM = 1};
64  getfem::mesh mesh; /* the mesh */
65  getfem::mesh_im mim; /* the integration methods */
66  getfem::mesh_fem mf_u; /* main mesh_fem, for the elastostatic solution */
67  getfem::mesh_fem mf_p; /* mesh_fem for the pressure. */
68  getfem::mesh_fem mf_rhs; /* mesh_fem for the right hand side (f(x),..) */
69  getfem::mesh_fem mf_coef; /* mesh_fem used to represent pde coefficients */
70  scalar_type p1, p2, p3; /* elastic coefficients. */
71 
72  scalar_type residual; /* max residual for the iterative solvers */
73 
74  std::string datafilename;
75  bgeot::md_param PARAM;
76 
77  bool solve(plain_vector &U);
78  void init(void);
79  elastostatic_problem(void) : mim(mesh), mf_u(mesh), mf_p(mesh), mf_rhs(mesh),
80  mf_coef(mesh) {}
81 };
82 
83 
84 /* Read parameters from the .param file, build the mesh, set finite element
85  and integration methods and selects the boundaries.
86 
87  (this is boilerplate code, not very interesting)
88  */
89 void elastostatic_problem::init(void) {
90  std::string MESH_TYPE = PARAM.string_value("MESH_TYPE","Mesh type ");
91  std::string FEM_TYPE = PARAM.string_value("FEM_TYPE","FEM name");
92  std::string FEM_TYPE_P= PARAM.string_value("FEM_TYPE_P",
93  "FEM name for the pressure");
94  std::string INTEGRATION = PARAM.string_value("INTEGRATION",
95  "Name of integration method");
96  cout << "MESH_TYPE=" << MESH_TYPE << "\n";
97  cout << "FEM_TYPE=" << FEM_TYPE << "\n";
98  cout << "INTEGRATION=" << INTEGRATION << "\n";
99 
100  /* First step : build the mesh */
101  bgeot::pgeometric_trans pgt =
102  bgeot::geometric_trans_descriptor(MESH_TYPE);
103  size_type N = pgt->dim();
104  std::vector<size_type> nsubdiv(N);
105  std::fill(nsubdiv.begin(),nsubdiv.end(),
106  PARAM.int_value("NX", "Nomber of space steps "));
107  nsubdiv[1] = PARAM.int_value("NY") ? PARAM.int_value("NY") : nsubdiv[0];
108  if (N>2)
109  nsubdiv[2] = PARAM.int_value("NZ") ? PARAM.int_value("NZ") : nsubdiv[0];
110  getfem::regular_unit_mesh(mesh, nsubdiv, pgt,
111  PARAM.int_value("MESH_NOISED") != 0);
112 
113  bgeot::base_matrix M(N,N);
114  for (size_type i=0; i < N; ++i) {
115  static const char *t[] = {"LX","LY","LZ"};
116  M(i,i) = (i<3) ? PARAM.real_value(t[i],t[i]) : 1.0;
117  }
118  if (N>1) { M(0,1) = PARAM.real_value("INCLINE") * PARAM.real_value("LY"); }
119 
120  /* scale the unit mesh to [LX,LY,..] and incline it */
121  mesh.transformation(M);
122 
123  datafilename = PARAM.string_value("ROOTFILENAME","Base name of data files.");
124  residual = PARAM.real_value("RESIDUAL"); if (residual == 0.) residual = 1e-10;
125 
126  p1 = PARAM.real_value("P1", "First Elastic coefficient");
127  p2 = PARAM.real_value("P2", "Second Elastic coefficient");
128  p3 = PARAM.real_value("P3", "Third Elastic coefficient");
129 
130  mf_u.set_qdim(dim_type(N));
131 
132  /* set the finite element on the mf_u */
133  getfem::pfem pf_u =
134  getfem::fem_descriptor(FEM_TYPE);
135  getfem::pintegration_method ppi =
136  getfem::int_method_descriptor(INTEGRATION);
137 
138  mim.set_integration_method(ppi);
139  mf_u.set_finite_element(pf_u);
140 
141  mf_p.set_finite_element(getfem::fem_descriptor(FEM_TYPE_P));
142 
143  /* set the finite element on mf_rhs (same as mf_u is DATA_FEM_TYPE is
144  not used in the .param file */
145  std::string data_fem_name = PARAM.string_value("DATA_FEM_TYPE");
146  if (data_fem_name.size() == 0) {
147  GMM_ASSERT1(pf_u->is_lagrange(), "You are using a non-lagrange FEM"
148  ". In that case you need to set "
149  << "DATA_FEM_TYPE in the .param file");
150  mf_rhs.set_finite_element(mesh.convex_index(), pf_u);
151  } else {
152  mf_rhs.set_finite_element(mesh.convex_index(),
153  getfem::fem_descriptor(data_fem_name));
154  }
155 
156  /* set the finite element on mf_coef. Here we use a very simple element
157  * since the only function that need to be interpolated on the mesh_fem
158  * is f(x)=1 ... */
159  mf_coef.set_finite_element(mesh.convex_index(),
160  getfem::classical_fem(pgt,0));
161 
162  /* set boundary conditions
163  * (Neuman on the upper face, Dirichlet elsewhere) */
164  cout << "Selecting Neumann and Dirichlet boundaries\n";
165  getfem::mesh_region border_faces;
166  getfem::outer_faces_of_mesh(mesh, border_faces);
167  for (getfem::mr_visitor it(border_faces); !it.finished(); ++it) {
168  assert(it.is_face());
169  base_node un = mesh.normal_of_face_of_convex(it.cv(), it.f());
170  un /= gmm::vect_norm2(un);
171  if (gmm::abs(un[N-1] - 1.0) < 1.0E-7) {
172  mesh.region(DIRICHLET_BOUNDARY_NUM).add(it.cv(),it.f());
173  } else if (gmm::abs(un[N-1] + 1.0) < 1.0E-7) {
174  mesh.region(DIRICHLET_BOUNDARY_NUM).add(it.cv(),it.f());
175  }
176  }
177 }
178 
179 /**************************************************************************/
180 /* Model. */
181 /**************************************************************************/
182 
183 bool elastostatic_problem::solve(plain_vector &U) {
184  size_type nb_dof_rhs = mf_rhs.nb_dof();
185  size_type N = mesh.dim();
186  size_type law_num = PARAM.int_value("LAW");
187  // Linearized elasticity brick.
188  base_vector p(3); p[0] = p1; p[1] = p2; p[2] = p3;
189 
190  /* choose the material law */
191  std::string lawname;
192  switch (law_num) {
193  case 0:
194  case 1: lawname = "Saint_Venant_Kirchhoff"; p.resize(2); break;
195  case 2: lawname = "Ciarlet_Geymonat"; p.resize(3); break;
196  case 3: lawname = "Incompressible_Mooney_Rivlin"; p.resize(2); break;
197  default: GMM_ASSERT1(false, "no such law");
198  }
199 
200  if (N == 2) {
201  cout << "2D plane strain hyper-elasticity\n";
202  lawname = "Plane_Strain_"+lawname;
203  }
204 
205  getfem::model model;
206 
207  // Main unknown of the problem (displacement).
208  model.add_fem_variable("u", mf_u);
209 
210  // Nonlinear elasticity brick
211  model.add_initialized_fixed_size_data("params", p);
212  getfem::add_finite_strain_elasticity_brick(model, mim, lawname, "u","params");
213 
214  // Incompressibility
215  if (law_num == 1 || law_num == 3) {
216  model.add_fem_variable("p", mf_p);
217  getfem::add_finite_strain_incompressibility_brick(model, mim, "u", "p");
218  }
219 
220  // Defining the volumic source term.
221  base_vector f(N);
222  f[0] = PARAM.real_value("FORCEX","Amplitude of the gravity");
223  f[1] = PARAM.real_value("FORCEY","Amplitude of the gravity");
224  if (N>2)
225  f[2] = PARAM.real_value("FORCEZ","Amplitude of the gravity");
226  plain_vector F(nb_dof_rhs * N);
227  for (size_type i = 0; i < nb_dof_rhs; ++i) {
228  gmm::copy(f, gmm::sub_vector(F, gmm::sub_interval(i*N, N)));
229  }
230  // Volumic source term brick.
231  model.add_initialized_fem_data("VolumicData", mf_rhs, F);
232  getfem::add_source_term_brick(model, mim, "u", "VolumicData");
233 
234  // Dirichlet condition
235  plain_vector F2(nb_dof_rhs * N);
236  model.add_initialized_fem_data("DirichletData", mf_rhs, F2);
237  if (PARAM.int_value("DIRICHLET_VERSION") == 0)
238  getfem::add_Dirichlet_condition_with_multipliers
239  (model, mim, "u", mf_u, DIRICHLET_BOUNDARY_NUM, "DirichletData");
240  else
241  getfem::add_Dirichlet_condition_with_penalization
242  (model, mim, "u", 1E15, DIRICHLET_BOUNDARY_NUM, "DirichletData");
243 
244 
245  gmm::iteration iter;
246  gmm::set_traces_level(1); // To avoid some trace messages
247 
248 
249  /* prepare the export routine for OpenDX (all time steps will be exported)
250  (can be viewed with "dx -edit nonlinear_elastostatic.net")
251  */
252  getfem::dx_export exp(datafilename + ".dx",
253  PARAM.int_value("VTK_EXPORT")==1);
254  getfem::stored_mesh_slice sl;
255  sl.build(mesh, getfem::slicer_boundary(mesh),8);
256  exp.exporting(sl,true); exp.exporting_mesh_edges();
257  //exp.begin_series("deformationsteps");
258  exp.write_point_data(mf_u, U, "stepinit");
259  exp.serie_add_object("deformationsteps");
260 
261  GMM_ASSERT1(!mf_rhs.is_reduced(), "To be adapted for reduced mesh_fems");
262 
263  int nb_step = int(PARAM.int_value("NBSTEP"));
264  size_type maxit = PARAM.int_value("MAXITER");
265 
266  for (int step = 0; step < nb_step; ++step) {
267  plain_vector DF(F);
268 
269  gmm::copy(gmm::scaled(F, (step+1.)/(scalar_type)nb_step), DF);
270  gmm::copy(DF, model.set_real_variable("VolumicData"));
271 
272  if (N>2) {
273  /* Apply the gradual torsion/extension */
274  scalar_type torsion = PARAM.real_value("TORSION",
275  "Amplitude of the torsion");
276  torsion *= (step+1)/scalar_type(nb_step);
277  scalar_type extension = PARAM.real_value("EXTENSION",
278  "Amplitude of the extension");
279  extension *= (step+1)/scalar_type(nb_step);
280  base_node G(N); G[0] = G[1] = 0.5;
281  for (size_type i = 0; i < nb_dof_rhs; ++i) {
282  const base_node P = mf_rhs.point_of_basic_dof(i) - G;
283  scalar_type r = sqrt(P[0]*P[0]+P[1]*P[1]),
284  theta = atan2(P[1],P[0]);
285  F2[i*N+0] = r*cos(theta + (torsion*P[2])) - P[0];
286  F2[i*N+1] = r*sin(theta + (torsion*P[2])) - P[1];
287  F2[i*N+2] = extension * P[2];
288  }
289  }
290  /* update the imposed displacement */
291  gmm::copy(F2, model.set_real_variable("DirichletData"));
292 
293  cout << "step " << step << ", number of variables : "
294  << model.nb_dof() << endl;
295  iter = gmm::iteration(residual, int(PARAM.int_value("NOISY")),
296  maxit ? maxit : 40000);
297 
298  /* let the non-linear solve (Newton) do its job */
299  getfem::standard_solve(model, iter);
300 
301  gmm::copy(model.real_variable("u"), U);
302  //char s[100]; sprintf(s, "step%d", step+1);
303 
304  /* append the new displacement to the exported opendx file */
305  exp.write_point_data(mf_u, U); //, s);
306  exp.serie_add_object("deformationsteps");
307  }
308 
309  // Solution extraction
310  gmm::copy(model.real_variable("u"), U);
311 
312  return (iter.converged());
313 }
314 
315 /**************************************************************************/
316 /* main program. */
317 /**************************************************************************/
318 
319 int main(int argc, char *argv[]) {
320 
321  GETFEM_MPI_INIT(argc, argv);
322  GMM_SET_EXCEPTION_DEBUG; // Exceptions make a memory fault, to debug.
323  FE_ENABLE_EXCEPT; // Enable floating point exception for Nan.
324 
325  try {
326  elastostatic_problem p;
327  p.PARAM.read_command_line(argc, argv);
328  p.init();
329  p.mesh.write_to_file(p.datafilename + ".mesh");
330  p.mf_u.write_to_file(p.datafilename + ".mf", true);
331  p.mf_rhs.write_to_file(p.datafilename + ".mfd", true);
332  plain_vector U(p.mf_u.nb_dof());
333  GMM_ASSERT1(p.solve(U), "Solve has failed");
334  if (p.PARAM.int_value("VTK_EXPORT")) {
335  gmm::vecsave(p.datafilename + ".U", U);
336  cout << "export to " << p.datafilename + ".vtk" << "..\n";
337  getfem::vtk_export exp(p.datafilename + ".vtk",
338  p.PARAM.int_value("VTK_EXPORT")==1);
339  exp.exporting(p.mf_u);
340  exp.write_point_data(p.mf_u, U, "elastostatic_displacement");
341  cout << "export done, you can view the data file with (for example)\n"
342  "mayavi2 -d " << p.datafilename << ".vtk -f ExtractVectorNorm -f "
343  "WarpVector -m Surface -m Outline\n";
344  }
345  }
346  GMM_STANDARD_CATCH_ERROR;
347 
348  GETFEM_MPI_FINALIZE;
349 
350  return 0;
351 }
getfem_export.h
Export solutions to various formats.
getfem_model_solvers.h
Standard solvers for model bricks.
getfem_nonlinear_elasticity.h
Non-linear elasticty and incompressibility bricks.
getfem::stored_mesh_slice
The output of a getfem::mesh_slicer which has been recorded.
Definition: getfem_mesh_slice.h:47
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
getfem::int_method_descriptor
pintegration_method int_method_descriptor(std::string name, bool throw_if_not_found=true)
Get an integration method from its name .
Definition: getfem_integration.cc:1130
getfem_regular_meshes.h
Build regular meshes.
getfem::mesh_im
Describe an integration method linked to a mesh.
Definition: getfem_mesh_im.h:47
getfem::add_Dirichlet_condition_with_multipliers
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.
Definition: getfem_models.cc:4671
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::model::add_initialized_fem_data
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...
Definition: getfem_models.h:835
bgeot::geometric_trans_descriptor
pgeometric_trans geometric_trans_descriptor(std::string name)
Get the geometric transformation from its string name.
Definition: bgeot_geometric_trans.cc:1163
getfem::classical_fem
pfem classical_fem(bgeot::pgeometric_trans pgt, short_type k, bool complete=false)
Give a pointer on the structures describing the classical polynomial fem of degree k on a given conve...
Definition: getfem_fem.cc:4141
getfem::mesh_fem
Describe a finite element method linked to a mesh.
Definition: getfem_mesh_fem.h:148
getfem::model
`‘Model’' variables store the variables, the data and the description of a model.
Definition: getfem_models.h:114
getfem_superlu.h
SuperLU interface for getfem.
gmm::iteration
The Iteration object calculates whether the solution has reached the desired accuracy,...
Definition: gmm_iter.h:53
getfem::model::add_initialized_fixed_size_data
void add_initialized_fixed_size_data(const std::string &name, const VECT &v)
Add a fixed size data (assumed to be a vector) to the model and initialized with v.
Definition: getfem_models.h:743
getfem::mesh_region::visitor
"iterator" class for regions.
Definition: getfem_mesh_region.h:237
getfem::model::add_fem_variable
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.
Definition: getfem_models.cc:834
gmm::vect_norm2
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norm2(const V &v)
Euclidean norm of a vector.
Definition: gmm_blas.h:557
getfem::pfem
std::shared_ptr< const getfem::virtual_fem > pfem
type of pointer on a fem description
Definition: getfem_fem.h:244
bgeot::small_vector
container for small vectors of POD (Plain Old Data) types.
Definition: bgeot_small_vector.h:205
getfem::add_finite_strain_elasticity_brick
size_type add_finite_strain_elasticity_brick(model &md, const mesh_im &mim, const std::string &lawname, const std::string &varname, const std::string &params, size_type region=size_type(-1))
Add a finite strain elasticity brick to the model with respect to the variable varname (the displacem...
Definition: getfem_nonlinear_elasticity.cc:2300
gmm::rsvector
sparse vector built upon std::vector.
Definition: gmm_def.h:488
getfem::stored_mesh_slice::build
void build(const getfem::mesh &m, const slicer_action &a, size_type nrefine=1)
Build the slice, by applying a slicer_action operation.
Definition: getfem_mesh_slice.h:197
getfem::add_Dirichlet_condition_with_penalization
size_type APIDECL add_Dirichlet_condition_with_penalization(model &md, const mesh_im &mim, const std::string &varname, scalar_type penalization_coeff, size_type region, const std::string &dataname=std::string(), const mesh_fem *mf_mult=0)
Add a Dirichlet condition on the variable varname and the mesh region region.
Definition: getfem_models.cc:4710
getfem::mesh
Describe a mesh (collection of convexes (elements) and points).
Definition: getfem_mesh.h:95
bgeot::pgeometric_trans
std::shared_ptr< const bgeot::geometric_trans > pgeometric_trans
pointer type for a geometric transformation
Definition: bgeot_geometric_trans.h:186
gmm.h
Include common gmm files.
getfem::mesh_region
structure used to hold a set of convexes and/or convex faces.
Definition: getfem_mesh_region.h:55
getfem::regular_unit_mesh
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.
Definition: getfem_regular_meshes.cc:238
getfem::add_source_term_brick
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.
Definition: getfem_models.cc:4127
getfem::slicer_boundary
Extraction of the boundary of a slice.
Definition: getfem_mesh_slicers.h:253
getfem::outer_faces_of_mesh
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:822
getfem::standard_solve
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.
Definition: getfem_model_solvers.cc:410
getfem::dx_export
A (quite large) class for exportation of data to IBM OpenDX.
Definition: getfem_export.h:318
getfem::fem_descriptor
pfem fem_descriptor(const std::string &name)
get a fem descriptor from its string name.
Definition: getfem_fem.cc:4232
getfem::model::set_real_variable
model_real_plain_vector & set_real_variable(const std::string &name, size_type niter) const
Gives the write access to the vector value of a variable.
Definition: getfem_models.cc:3013
getfem::model::nb_dof
size_type nb_dof(bool with_internal=false) const
Total number of degrees of freedom in the model.
Definition: getfem_models.cc:295
getfem::model::real_variable
const model_real_plain_vector & real_variable(const std::string &name, size_type niter) const
Gives the access to the vector value of a variable.
Definition: getfem_models.cc:2959
getfem::vtk_export
VTK/VTU export.
Definition: getfem_export.h:68