GetFEM  5.4.2
plate.cc
Go to the documentation of this file.
1 /*===========================================================================
2 
3  Copyright (C) 2002-2020 Yves Renard, Michel Salaün.
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 plate.cc
24  @brief Linear Elastostatic plate problem.
25 
26  This program is used to check that getfem++ is working. This is
27  also a good example of use of GetFEM.
28 */
29 
30 #include "getfem/getfem_assembling.h" /* import assembly methods (and norms comp.) */
31 #include "getfem/getfem_linearized_plates.h"
32 #include "getfem/getfem_export.h" /* export functions (save solution in a file) */
33 #include "getfem/getfem_regular_meshes.h"
34 #include "getfem/getfem_model_solvers.h"
35 #include "gmm/gmm.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 using bgeot::base_small_vector; /* special class for small (dim<16) vectors */
42 using bgeot::base_node; /* geometrical nodes(derived from base_small_vector)*/
43 using bgeot::scalar_type; /* = double */
44 using bgeot::size_type; /* = unsigned long */
45 using bgeot::dim_type;
46 using bgeot::base_matrix; /* small dense matrix. */
47 
48 /* definition of some matrix/vector types. These ones are built
49  * using the predefined types in Gmm++
50  */
51 typedef getfem::modeling_standard_sparse_vector sparse_vector;
52 typedef getfem::modeling_standard_sparse_matrix sparse_matrix;
53 typedef getfem::modeling_standard_plain_vector plain_vector;
54 
55 /*
56  structure for the elastostatic problem
57 */
58 struct plate_problem {
59 
60  enum { SIMPLY_FIXED_BOUNDARY_NUM = 0 };
61  getfem::mesh mesh; /* the mesh */
62  getfem::mesh_im mim, mim_subint;
63  getfem::mesh_fem mf_ut;
64  getfem::mesh_fem mf_u3;
65  getfem::mesh_fem mf_theta;
66  getfem::mesh_fem mf_rhs; /* mesh_fem for the right hand side (f(x),..) */
67  getfem::mesh_fem mf_coef; /* mesh_fem used to represent pde coefficients */
68  scalar_type lambda, mu; /* Lamé coefficients. */
69  scalar_type E, nu; /* Lamé coefficients. */
70  scalar_type epsilon; /* thickness of the plate. */
71  scalar_type pressure;
72  scalar_type residual; /* max residual for the iterative solvers */
73  scalar_type LX , LY ; // default : LX = LY = 1
74  bool mitc;
75  int sol_ref; // sol_ref = 0 : simple support on the vertical edges
76  // sol_ref = 1 : homogeneous on the vertical edges
77  // sol_ref = 2 : homogeneous on the 4 vertical
78  // edges with solution u3 = sin²(x)*sin²(y)
79  scalar_type eta; // useful only if sol_ref == 2 :
80  // eta = 0 => Kirchhoff-Love
81  // eta = small => Mindlin
82  size_type N_Four ;
83  base_matrix theta1_Four, theta2_Four, u3_Four ;
84 
85  int study_flag; // if studyflag = 1, then the loadings applied are chosen
86  // in order to have a maximal vertical displacement equal to one.
87  // Nothing is done if study_flag has another value.
88 
89  std::string datafilename;
90  bgeot::md_param PARAM;
91 
92  base_small_vector theta_exact(base_node P);
93  scalar_type u3_exact(base_node P);
94 
95  bool solve(plain_vector &Ut, plain_vector &U3, plain_vector &THETA);
96  void init(void);
97  void compute_error(plain_vector &Ut, plain_vector &U3, plain_vector &THETA);
98  plate_problem(void) : mim(mesh), mim_subint(mesh), mf_ut(mesh), mf_u3(mesh),
99  mf_theta(mesh), mf_rhs(mesh), mf_coef(mesh) {}
100 };
101 
102 /* Read parameters from the .param file, build the mesh, set finite element
103  * and integration methods and selects the boundaries.
104  */
105 void plate_problem::init(void) {
106  std::string MESH_FILE = PARAM.string_value("MESH_FILE");
107  std::string MESH_TYPE = PARAM.string_value("MESH_TYPE","Mesh type ");
108  std::string FEM_TYPE_UT = PARAM.string_value("FEM_TYPE_UT","FEM name");
109  std::string FEM_TYPE_U3 = PARAM.string_value("FEM_TYPE_U3","FEM name");
110  std::string FEM_TYPE_THETA = PARAM.string_value("FEM_TYPE_THETA","FEM name");
111  std::string INTEGRATION = PARAM.string_value("INTEGRATION",
112  "Name of integration method");
113  std::string INTEGRATION_CT = PARAM.string_value("INTEGRATION_CT",
114  "Name of integration method");
115  cout << "MESH_TYPE=" << MESH_TYPE << "\n";
116  cout << "FEM_TYPE_UT=" << FEM_TYPE_UT << "\n";
117  cout << "INTEGRATION=" << INTEGRATION << "\n";
118  cout << "INTEGRATION_CT=" << INTEGRATION_CT << "\n";
119 
120  /* First step : build the mesh */
121  size_type N;
122  bgeot::pgeometric_trans pgt =
123  bgeot::geometric_trans_descriptor(MESH_TYPE);
124  if (!MESH_FILE.empty()) {
125  cout << "MESH_FILE=" << MESH_FILE << "\n";
126  mesh.read_from_file(MESH_FILE);
127  MESH_TYPE = bgeot::name_of_geometric_trans
128  (mesh.trans_of_convex(mesh.convex_index().first_true()));
129  cout << "MESH_TYPE=" << MESH_TYPE << "\n";
130  N = mesh.dim();
131  } else {
132  N = pgt->dim();
133  GMM_ASSERT1(N == 2, "For a plate problem, N should be 2");
134  std::vector<size_type> nsubdiv(N);
135  std::fill(nsubdiv.begin(),nsubdiv.end(),
136  PARAM.int_value("NX", "Number of space steps "));
137  getfem::regular_unit_mesh(mesh, nsubdiv, pgt,
138  PARAM.int_value("MESH_NOISED") != 0);
139  }
140 
141  datafilename = PARAM.string_value("ROOTFILENAME","Base name of data files.");
142  residual = PARAM.real_value("RESIDUAL"); if (residual == 0.) residual = 1e-10;
143  mitc = (PARAM.int_value("MITC", "Mitc version ?") != 0);
144 
145 
146  cout << "MITC = " ;
147  if (mitc) cout << "true \n" ; else cout << "false \n" ;
148  sol_ref = int(PARAM.int_value("SOL_REF") ) ;
149  study_flag = int(PARAM.int_value("STUDY_FLAG") ) ;
150  eta = (PARAM.real_value("ETA") );
151  N_Four = (PARAM.int_value("N_Four") ) ;
152 
153 
154  LX = PARAM.real_value("LX");
155  LY = PARAM.real_value("LY");
156  mu = PARAM.real_value("MU", "Lamé coefficient mu");
157  lambda = PARAM.real_value("LAMBDA", "Lamé coefficient lambda");
158  epsilon = PARAM.real_value("EPSILON", "thickness of the plate");
159  pressure = PARAM.real_value("PRESSURE",
160  "pressure on the top surface of the plate.");
161 
162 
163  cout << "SOL_REF = " ;
164  if (sol_ref==0) cout << "appui simple aux 2 bords verticaux\n" ;
165  if (sol_ref==1) cout << "encastrement aux 2 bords verticaux\n" ;
166  if (sol_ref==2) {
167  cout << "encastrement aux 4 bords verticaux, solution en sin(x)^2*sin(y)^2\n" ;
168  cout << "eta = " << eta <<"\n";
169  }
170  if (sol_ref==4) {
171  cout << "bord en appuis simple\n" ;
172  cout << "nombre de terme pour calcul sol exacte : " << N_Four << " \n" ;
173  // Calcul des coeeficients de Fourier de la solution exacte :
174  // Cas où le chargement est seulement vertical (pas de moment appliqué)
175  gmm::resize( theta1_Four, N_Four, N_Four) ;
176  gmm::resize( theta2_Four, N_Four, N_Four) ;
177  gmm::resize( u3_Four, N_Four, N_Four) ;
178  base_matrix Jmn(3, 3) ;
179  base_small_vector Bmn(3), Xmn(3) ;
180  scalar_type /*det_Jmn, */ A, B, e2, Pmn ;
181  E = 4.*mu*(mu+lambda) / (2. * mu + lambda);
182  nu = lambda / (2. * mu + lambda);
183  e2 = epsilon * epsilon / 4. ;
184  for(size_type i = 0 ; i < N_Four ; i++) {
185  for(size_type j = 0 ; j < N_Four ; j++) {
186  A = scalar_type(j + 1) * M_PI / LX ;
187  B = scalar_type(i + 1) * M_PI / LY ;
188  Jmn(0, 0) = 2. * A * A / (1. - nu) + B * B + 3. / e2 ;
189  Jmn(0, 1) = A * B * (1. +nu) / (1. - nu) ;
190  Jmn(0, 2) = A * 3. / e2 ;
191  Jmn(1, 0) = A * B * (1. +nu) / (1. - nu) ;
192  Jmn(1, 1) = 2. * B * B / (1. - nu) + A * A + 3. / e2 ;
193  Jmn(1, 2) = B * 3. / e2 ;
194  Jmn(2, 0) = - A ;
195  Jmn(2, 1) = - B ;
196  Jmn(2, 2) = A * A + B * B ;
197  gmm::scale(Jmn, - E*(epsilon/2.) / (1. + nu) ) ;
198 
199  // calcul du développement de Fourrier du chargement :
200  if ( ( (i + 1) % 2 == 1 ) && ( (j + 1) % 2 == 1) ) {
201  Pmn = 16. * pressure / ( scalar_type(i + 1) * scalar_type(j + 1) * M_PI * M_PI) ; }
202  else {
203  Pmn = 0. ; }
204  Bmn[0] = 0. ;
205  Bmn[1] = 0. ;
206  Bmn[2] = Pmn ;
207  gmm::lu_solve(Jmn, Xmn, Bmn) ;
208  theta1_Four(i, j) = Xmn[0] ;
209  theta1_Four(i, j) = Xmn[1] ;
210  u3_Four(i, j) = Xmn[2] ;
211  }
212  }
213  }
214 
215  mf_ut.set_qdim(dim_type(N));
216  mf_theta.set_qdim(dim_type(N));
217 
218  /* set the finite element on the mf_u */
219  getfem::pfem pf_ut = getfem::fem_descriptor(FEM_TYPE_UT);
220  getfem::pfem pf_u3 = getfem::fem_descriptor(FEM_TYPE_U3);
221  getfem::pfem pf_theta = getfem::fem_descriptor(FEM_TYPE_THETA);
222 
223  getfem::pintegration_method ppi =
224  getfem::int_method_descriptor(INTEGRATION);
225  getfem::pintegration_method ppi_ct =
226  getfem::int_method_descriptor(INTEGRATION_CT);
227 
228  mim.set_integration_method(mesh.convex_index(), ppi);
229  mim_subint.set_integration_method(mesh.convex_index(), ppi_ct);
230  mf_ut.set_finite_element(mesh.convex_index(), pf_ut);
231  mf_u3.set_finite_element(mesh.convex_index(), pf_u3);
232  mf_theta.set_finite_element(mesh.convex_index(), pf_theta);
233 
234  /* set the finite element on mf_rhs (same as mf_u is DATA_FEM_TYPE is
235  not used in the .param file */
236  std::string data_fem_name = PARAM.string_value("DATA_FEM_TYPE");
237  if (data_fem_name.size() == 0) {
238  GMM_ASSERT1(pf_ut->is_lagrange(), "You are using a non-lagrange FEM. "
239  << "In that case you need to set "
240  << "DATA_FEM_TYPE in the .param file");
241  mf_rhs.set_finite_element(mesh.convex_index(), pf_ut);
242  } else {
243  mf_rhs.set_finite_element(mesh.convex_index(),
244  getfem::fem_descriptor(data_fem_name));
245  }
246 
247  /* set the finite element on mf_coef. Here we use a very simple element
248  * since the only function that need to be interpolated on the mesh_fem
249  * is f(x)=1 ... */
250  mf_coef.set_finite_element(mesh.convex_index(),
251  getfem::classical_fem(pgt,0));
252 
253  /* set boundary conditions
254  * (Neuman on the upper face, Dirichlet elsewhere) */
255  cout << "Selecting Neumann and Dirichlet boundaries\n";
256  getfem::mesh_region border_faces;
257  getfem::outer_faces_of_mesh(mesh, border_faces);
258  for (getfem::mr_visitor i(border_faces); !i.finished(); ++i) {
259  assert(i.is_face());
260  base_node un = mesh.normal_of_face_of_convex(i.cv(), i.f());
261  un /= gmm::vect_norm2(un);
262  switch(sol_ref){
263  case 0 :
264  if (gmm::abs(un[1]) <= 1.0E-7) // new Neumann face
265  mesh.region(SIMPLY_FIXED_BOUNDARY_NUM).add(i.cv(), i.f());
266  break ;
267  case 1 :
268  if (gmm::abs(un[1]) <= 1.0E-7) // new Neumann face
269  mesh.region(SIMPLY_FIXED_BOUNDARY_NUM).add(i.cv(), i.f());
270  break ;
271  case 2 :
272  if ( (gmm::abs(un[0]) <= 1.0E-7) || (gmm::abs(un[1]) <= 1.0E-7) )
273  mesh.region(SIMPLY_FIXED_BOUNDARY_NUM).add(i.cv(), i.f());
274  break ;
275  case 3 :
276  if (un[0] <= (- 1. + 1.0E-7)) // new Neumann face
277  mesh.region(SIMPLY_FIXED_BOUNDARY_NUM).add(i.cv(), i.f());
278  break ;
279  case 4 :
280  if ( (gmm::abs(un[0]) <= 1.0E-7) || (gmm::abs(un[1]) <= 1.0E-7) )
281  mesh.region(SIMPLY_FIXED_BOUNDARY_NUM).add(i.cv(), i.f());
282  break ;
283  default :
284  GMM_ASSERT1(false, "SOL_REF parameter is undefined");
285  break ;
286  }
287  }
288 }
289 
290 base_small_vector plate_problem::theta_exact(base_node P) {
291  base_small_vector theta(2);
292  if (sol_ref == 0) { // appui simple aux 2 bords
293  theta[0] = - (-pressure / (32. * mu * epsilon * epsilon * epsilon / 8.))
294  * (4. * pow(P[0] - .5, 3.) - 3 * (P[0] - .5));
295  theta[1] = 0.;
296  }
297  if (sol_ref == 1) { // encastrement aux 2 bords
298  theta[0] = - (-pressure / (16. * mu * epsilon * epsilon * epsilon / 8.))
299  * P[0] * ( 2.*P[0]*P[0] - 3.* P[0] + 1. ) ;
300  theta[1] = 0.;
301  }
302  if (sol_ref == 2) { // encastrement aux 4 bords et sols vraiment 2D, non polynomiale
303  theta[0] = (1. + eta) * M_PI * sin(2.*M_PI*P[0]) * sin(M_PI*P[1])* sin(M_PI*P[1]) ;
304  theta[1] = (1. + eta) * M_PI * sin(2.*M_PI*P[1]) * sin(M_PI*P[0])* sin(M_PI*P[0]) ;
305  }
306  if (sol_ref == 3) { // plaque cantilever
307  theta[0] = - (- 3. * pressure / (8. * mu * epsilon * epsilon * epsilon / 8.))
308  * P[0] * ( 0.25 * P[0] * P[0] - P[0] + 1. ) ;
309  theta[1] = 0.;
310  }
311  if (sol_ref == 4) { // bord entier en appui simple
312  theta[0] = 0. ;
313  theta[1] = 0. ;
314  for(size_type i = 0 ; i < N_Four ; i ++) {
315  for(size_type j = 0 ; j < N_Four ; j ++) {
316  theta[0] -= theta1_Four(i, j) * cos( scalar_type(j + 1) * M_PI * P[0] / LX ) * sin( scalar_type(i + 1) * M_PI * P[1] / LY ) ;
317  theta[0] -= theta2_Four(i, j) * sin( scalar_type(j + 1) * M_PI * P[0] / LX ) * cos( scalar_type(i + 1) * M_PI * P[1] / LY ) ;
318  }
319  }
320  }
321  return theta;
322 }
323 
324 scalar_type plate_problem::u3_exact(base_node P) {
325  switch(sol_ref) {
326  case 0 : return (pressure / (32. * mu * epsilon * epsilon * epsilon / 8.))
327  * P[0] * (P[0] - 1.)
328  * (gmm::sqr(P[0] - .5) -1.25-(8.* epsilon*epsilon / 4.));
329  break ;
330  case 1 : return (pressure /(32.* mu * epsilon * epsilon * epsilon / 8.))
331  * P[0] * (P[0] - 1.)
332  * ( P[0] * P[0] - P[0] - 8. * epsilon *epsilon / 4.) ;
333  break ;
334  case 2 : return gmm::sqr(sin(M_PI*P[0])) * gmm::sqr(sin(M_PI*P[1]));
335  break ;
336  case 3 : return (3. * pressure / (4. * mu * epsilon * epsilon * epsilon / 8. ))
337  * P[0] * ( P[0] * P[0] * P[0] / 24. - P[0] * P[0] / 6. + P[0] / 4.
338  - (epsilon * epsilon / 4.) * P[0] / 3.
339  + 2. * (epsilon * epsilon / 4.) / 3.) ;
340  break ;
341  case 4 :
342  scalar_type u3_local ;
343  u3_local = 0. ;
344  for(size_type i = 0 ; i < N_Four ; i ++) {
345  for(size_type j = 0 ; j < N_Four ; j ++)
346  u3_local += u3_Four(i, j) * sin( scalar_type(j + 1) * M_PI * P[0] / LX ) * sin( scalar_type(i + 1) * M_PI * P[1] / LY ) ;
347  }
348  return (u3_local) ;
349  break ;
350  default : GMM_ASSERT1(false, "indice de solution de référence incorrect");
351  }
352 }
353 
354 
355 /* compute the error with respect to the exact solution */
356 void plate_problem::compute_error(plain_vector &Ut, plain_vector &U3, plain_vector &THETA) {
357  cout.precision(16);
358  if (PARAM.int_value("SOL_EXACTE") == 1) {
359  gmm::clear(Ut); gmm::clear(U3); gmm::clear(THETA);
360  }
361 
362  std::vector<scalar_type> V(mf_rhs.nb_dof()*2);
363 
364  getfem::interpolation(mf_ut, mf_rhs, Ut, V);
365  mf_rhs.set_qdim(2);
366  scalar_type l2 = gmm::sqr(getfem::asm_L2_norm(mim, mf_rhs, V));
367  scalar_type h1 = gmm::sqr(getfem::asm_H1_norm(mim, mf_rhs, V));
368  scalar_type linf = gmm::vect_norminf(V);
369  mf_rhs.set_qdim(1);
370  cout << "L2 error = " << sqrt(l2) << endl
371  << "H1 error = " << sqrt(h1) << endl
372  << "Linfty error = " << linf << endl;
373 
374  getfem::interpolation(mf_theta, mf_rhs, THETA, V);
375  GMM_ASSERT1(!mf_rhs.is_reduced(),
376  "To be adapted, use interpolation_function");
377  for (size_type i = 0; i < mf_rhs.nb_dof(); ++i) {
378  gmm::add(gmm::scaled(theta_exact(mf_rhs.point_of_basic_dof(i)), -1.0),
379  gmm::sub_vector(V, gmm::sub_interval(i*2, 2)));
380  }
381  mf_rhs.set_qdim(2);
382  l2 += gmm::sqr(getfem::asm_L2_norm(mim, mf_rhs, V));
383  h1 += gmm::sqr(getfem::asm_H1_semi_norm(mim, mf_rhs, V));
384  linf = std::max(linf, gmm::vect_norminf(V));
385  mf_rhs.set_qdim(1);
386  cout << "L2 error theta:" << sqrt(l2) << endl
387  << "H1 error theta:" << sqrt(h1) << endl
388  << "Linfty error = " << linf << endl;
389 
390  gmm::resize(V, mf_rhs.nb_dof());
391  getfem::interpolation(mf_u3, mf_rhs, U3, V);
392 
393  for (size_type i = 0; i < mf_rhs.nb_dof(); ++i)
394  V[i] -= u3_exact(mf_rhs.point_of_basic_dof(i));
395 
396  l2 = gmm::sqr(getfem::asm_L2_norm(mim, mf_rhs, V));
397  h1 = gmm::sqr(getfem::asm_H1_semi_norm(mim, mf_rhs, V));
398  linf = std::max(linf, gmm::vect_norminf(V));
399 
400  cout.precision(16);
401  cout << "L2 error u3:" << sqrt(l2) << endl
402  << "H1 error u3:" << sqrt(h1) << endl
403  << "Linfty error = " << linf << endl;
404 
405  // stockage de l'erreur H1 :
406  if (PARAM.int_value("SAUV")){
407  std::ofstream f_out("errH1.don");
408  if (!f_out) throw std :: runtime_error("Impossible to open file") ;
409  f_out << sqrt(h1) << "\n" ;
410  }
411 
412 
413 }
414 
415 /**************************************************************************/
416 /* Model. */
417 /**************************************************************************/
418 
419 bool plate_problem::solve(plain_vector &Ut, plain_vector &U3, plain_vector &THETA) {
420  size_type nb_dof_rhs = mf_rhs.nb_dof();
421 
422  cout << "Number of dof for ut: " << mf_ut.nb_dof() << endl;
423  cout << "Number of dof for u3: " << mf_u3.nb_dof() << endl;
424  cout << "Number of dof for theta: " << mf_theta.nb_dof() << endl;
425 
426  E = 4.*mu*(mu+lambda) / (2. * mu + lambda);
427  nu = lambda / (2. * mu + lambda);
428  scalar_type kappa = 5./6.;
429 
430  getfem::model md;
431  md.add_fem_variable("ut", mf_ut);
432  md.add_fem_variable("u3", mf_u3);
433  md.add_fem_variable("theta", mf_theta);
434 
435  // Linearized plate brick.
436  md.add_initialized_scalar_data("E", E);
437  md.add_initialized_scalar_data("nu", nu);
438  md.add_initialized_scalar_data("lambda", lambda);
439  md.add_initialized_scalar_data("mu", mu);
440  md.add_initialized_scalar_data("epsilon", epsilon);
441  md.add_initialized_scalar_data("kappa", kappa);
442  getfem::add_Mindlin_Reissner_plate_brick(md, mim, mim_subint, "u3", "theta",
443  "E", "nu", "epsilon", "kappa",
444  (mitc) ? 2 : 1);
445  getfem::add_isotropic_linearized_elasticity_brick(md, mim, "ut", "lambda", "mu");
446 
447  // Defining the surface source term.
448  if (study_flag == 1 ){
449  cout << "Attention : l'intensité de la pression verticale " ;
450  cout << "a été choisie pour que le déplacement maximal soit unitaire." ;
451  cout << "Pour annuler cette option, faire STUDY_FLAG = 0\n" ;
452  switch(sol_ref) {
453  case 0 :
454  pressure = 128. * mu * epsilon * epsilon * epsilon / 8. ;
455  pressure /= 1.25 + 8. * epsilon * epsilon / 4. ;
456  break ;
457  case 1 :
458  pressure = 128. * mu * epsilon * epsilon * epsilon / 8. ;
459  pressure /= 0.25 + 8. * epsilon * epsilon / 4. ;
460  break ;
461  case 3 :
462  pressure = 32. * mu * epsilon * epsilon * epsilon / 8.;
463  pressure /= 3. + 8. * epsilon * epsilon / 4.;
464  default :
465  break ;
466  }
467  }
468  plain_vector F(nb_dof_rhs);
469  plain_vector M(nb_dof_rhs * 2);
470  if (sol_ref == 2) {
471  base_small_vector P(2) ;
472  scalar_type sx, sy, cx, cy, s2x, s2y, c2x, c2y ;
473  E = 4.*mu*(mu+lambda) / (2. * mu + lambda);
474  nu = lambda / (2. * mu + lambda);
475  for (size_type i = 0; i < nb_dof_rhs; ++i) {
476  P = mf_rhs.point_of_basic_dof(i);
477  sx = sin(M_PI*P[0]) ;
478  cx = cos(M_PI*P[0]) ;
479  sy = sin(M_PI*P[1]) ;
480  cy = cos(M_PI*P[1]) ;
481  c2x = cos(2.*M_PI*P[0]) ;
482  c2y = cos(2.*M_PI*P[1]) ;
483  s2x = sin(2.*M_PI*P[0]) ;
484  s2y = sin(2.*M_PI*P[1]) ;
485  F[i] = 2. * (epsilon / 2.) * E * M_PI * M_PI * eta *
486  ( sy * sy * c2x + sx * sx * c2y ) / ( 1. + nu ) ;
487  M[2*i] = -((epsilon * epsilon * epsilon / 8.) * E * M_PI * s2x / 3. / (1. + nu))
488  * ( (4. * M_PI * M_PI * (1. + eta) * (2. * c2y - 1.) / (1.- nu))
489  - 3. * eta * sy * sy / (epsilon/2.) / (epsilon/2.) ) ;
490  M[2*i+1] = -((epsilon * epsilon * epsilon/8.) * E * M_PI * s2y / 3. / (1. + nu))
491  * ( (4. * M_PI * M_PI * (1. + eta) * (2. * c2x - 1.) / (1.- nu))
492  - 3. * eta * sx * sx / (epsilon/2.) / (epsilon/2.) ) ;
493  }
494  }
495  else // sol_ref = 0 ou 1 ou 3 ou 4: pression verticale uniforme
496  for (size_type i = 0; i < nb_dof_rhs; ++i) {
497  F[i] = pressure;
498  }
499 
500  md.add_initialized_fem_data("VF", mf_rhs, F);
501  getfem::add_source_term_brick(md, mim, "u3", "VF");
502  md.add_initialized_fem_data("VM", mf_rhs, M);
503  getfem::add_source_term_brick(md, mim, "theta", "VM");
504 
505  getfem::add_Dirichlet_condition_with_multipliers
506  (md, mim, "u3", mf_u3, SIMPLY_FIXED_BOUNDARY_NUM);
507  getfem::add_Dirichlet_condition_with_multipliers
508  (md, mim, "ut", mf_ut, SIMPLY_FIXED_BOUNDARY_NUM);
509 
510  if (sol_ref == 1 || sol_ref == 2 || sol_ref == 3)
511  getfem::add_Dirichlet_condition_with_multipliers
512  (md, mim, "theta", mf_ut, SIMPLY_FIXED_BOUNDARY_NUM);
513 
514 
515  // Generic solve.
516  gmm::iteration iter(residual, 1, 40000);
517  getfem::standard_solve(md, iter);
518 
519 
520  gmm::resize(Ut, mf_ut.nb_dof());
521  gmm::copy(md.real_variable("ut"), Ut);
522  gmm::resize(U3, mf_u3.nb_dof());
523  gmm::copy(md.real_variable("u3"), U3);
524  gmm::resize(THETA, mf_theta.nb_dof());
525  gmm::copy(md.real_variable("theta"), THETA);
526 
527  if (PARAM.int_value("VTK_EXPORT")) {
528  cout << "export to " << datafilename + ".vtk" << "..\n";
529  getfem::vtk_export exp(datafilename + ".vtk",
530  PARAM.int_value("VTK_EXPORT")==1);
531  exp.exporting(mf_u3);
532  exp.write_point_data(mf_u3, U3, "plate_normal_displacement");
533  cout << "export done, you can view the data file with (for example)\n"
534  "mayavi2 -d " << datafilename << ".vtk -f "
535  "WarpScalar -m Surface -m Outline\n";
536 // cout << "export done, you can view the data file with (for example)\n"
537 // "mayavi -d " << datafilename << ".vtk -f ExtractVectorNorm -f "
538 // "WarpVector -m BandedSurfaceMap -m Outline\n";
539  }
540  if (PARAM.int_value("DX_EXPORT")) {
541  cout << "export to " << datafilename + ".dx" << ".\n";
542  getfem::dx_export exp(datafilename + ".dx",
543  PARAM.int_value("DX_EXPORT")==1);
544  exp.exporting(mf_u3);
545  exp.write_point_data(mf_u3, U3, "plate_normal_displacement");
546  }
547 
548 
549  return (iter.converged());
550 }
551 
552 /**************************************************************************/
553 /* main program. */
554 /**************************************************************************/
555 
556 int main(int argc, char *argv[]) {
557 
558  GETFEM_MPI_INIT(argc, argv);
559  GMM_SET_EXCEPTION_DEBUG; // Exceptions make a memory fault, to debug.
560  FE_ENABLE_EXCEPT; // Enable floating point exception for Nan.
561 
562  try {
563  plate_problem p;
564  p.PARAM.read_command_line(argc, argv);
565  p.init();
566  if ((p.study_flag != 1)&&((p.sol_ref == 0) || (p.sol_ref ==1)))
567  p.pressure *= p.epsilon * p.epsilon * p.epsilon / 8.;
568  p.mesh.write_to_file(p.datafilename + ".mesh");
569  plain_vector Ut, U3, THETA;
570  bool ok = p.solve(Ut, U3, THETA);
571  p.compute_error(Ut, U3, THETA);
572  GMM_ASSERT1(ok, "Solve has failed");
573 
574  }
575  GMM_STANDARD_CATCH_ERROR;
576 
577  GETFEM_MPI_FINALIZE;
578 
579  return 0;
580 }
getfem_export.h
Export solutions to various formats.
getfem_model_solvers.h
Standard solvers for model bricks.
getfem::add_Mindlin_Reissner_plate_brick
size_type add_Mindlin_Reissner_plate_brick(model &md, const mesh_im &mim, const mesh_im &mim_reduced, const std::string &u3, const std::string &Theta, const std::string &param_E, const std::string &param_nu, const std::string &param_epsilon, const std::string &param_kappa, size_type variant=size_type(2), size_type region=size_type(-1))
Add a term corresponding to the classical Reissner-Mindlin plate model for which u3 is the transverse...
Definition: getfem_linearized_plates.cc:28
bgeot::name_of_geometric_trans
std::string name_of_geometric_trans(pgeometric_trans p)
Get the string name of a geometric transformation.
Definition: bgeot_geometric_trans.cc:1173
getfem::asm_H1_norm
scalar_type asm_H1_norm(const mesh_im &mim, const mesh_fem &mf, const VEC &U, const mesh_region &rg=mesh_region::all_convexes())
compute the H1 norm of U.
Definition: getfem_assembling.h:302
gmm::resize
void resize(M &v, size_type m, size_type n)
*â€/
Definition: gmm_blas.h:231
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
getfem::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.
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:59
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::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::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
gmm::iteration
The Iteration object calculates whether the solution has reached the desired accuracy,...
Definition: gmm_iter.h:53
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_linearized_plates.h
Reissner-Mindlin plate model brick.
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::asm_L2_norm
scalar_type asm_L2_norm(const mesh_im &mim, const mesh_fem &mf, const VEC &U, const mesh_region &rg=mesh_region::all_convexes())
compute , U might be real or complex
Definition: getfem_assembling.h:56
gmm::rsvector
sparse vector built upon std::vector.
Definition: gmm_def.h:488
getfem::model::add_initialized_scalar_data
void add_initialized_scalar_data(const std::string &name, T e)
Add a scalar data (i.e.
Definition: getfem_models.h:782
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::add_isotropic_linearized_elasticity_brick
size_type APIDECL add_isotropic_linearized_elasticity_brick(model &md, const mesh_im &mim, const std::string &varname, const std::string &dataname_lambda, const std::string &dataname_mu, size_type region=size_type(-1), const std::string &dataname_preconstraint=std::string())
Linear elasticity brick ( ).
Definition: getfem_models.cc:6106
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
gmm::vect_norminf
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norminf(const V &v)
Infinity norm of a vector.
Definition: gmm_blas.h:693
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::asm_H1_semi_norm
scalar_type asm_H1_semi_norm(const mesh_im &mim, const mesh_fem &mf, const VEC &U, const mesh_region &rg=mesh_region::all_convexes())
compute , U might be real or complex
Definition: getfem_assembling.h:179
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::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_assembling.h
Miscelleanous assembly routines for common terms. Use the low-level generic assembly....
getfem::vtk_export
VTK/VTU export.
Definition: getfem_export.h:68