GetFEM++  5.3
elastostatic.cc
Go to the documentation of this file.
1 /*===========================================================================
2 
3  Copyright (C) 2002-2017 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 elastostatic.cc
24  @brief Linear Elastostatic problem. A dummy linear
25  elastotatic problem is solved on a regular mesh, and is compared to
26  the analytical solution.
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  @see laplacian.cc
32  @see nonlinear_elastostatic.cc
33 */
34 
35 #include "getfem/getfem_config.h"
36 #include "getfem/getfem_assembling.h" /* import assembly methods (and norms comp.) */
37 #include "getfem/getfem_export.h" /* export functions (save solution in a file) */
40 #include "gmm/gmm.h"
43 #include "getfem/getfem_import.h"
44 
45 using std::endl; using std::cout; using std::cerr;
46 using std::ends; using std::cin;
47 
48 
49 /* some GetFEM++ types that we will be using */
50 using bgeot::base_small_vector; /* special class for small (dim<16) vectors */
51 using bgeot::base_node; /* geometrical nodes(derived from base_small_vector)*/
52 using bgeot::scalar_type; /* = double */
53 using bgeot::size_type; /* = unsigned long */
54 using bgeot::dim_type;
55 using bgeot::base_matrix; /* small dense matrix. */
56 
57 /* definition of some matrix/vector types.
58  * default types of getfem_model_solvers.h
59  */
61 typedef getfem::modeling_standard_sparse_matrix sparse_matrix;
62 typedef getfem::modeling_standard_plain_vector plain_vector;
63 
64 /**************************************************************************/
65 /* Exact solution. */
66 /**************************************************************************/
67 
68 gmm::row_matrix<base_small_vector> sol_K;
69 static scalar_type sol_lambda, sol_mu, alph = 0.3;
70 int sol_sing;
71 
72 base_small_vector sol_u(const base_node &x) {
73  int N = x.size(); base_small_vector res(N);
74  switch(sol_sing) {
75  case 0 :
76  for (int i = 0; i < N; ++i)
77  res[i] = alph * sin(gmm::vect_sp(sol_K.row(i), x));
78  break;
79  case 1 :
80  {
81  base_small_vector trans(x.size());
82  gmm::fill(trans, M_PI / 10.0);
83  base_node y = x - trans;
84  scalar_type a = gmm::vect_norm2(y);
85  for (int i = 0; i < N; ++i) res[i] = a;
86  break;
87  }
88  case 2 :
89  {
90  base_small_vector trans(x.size());
91  gmm::fill(trans, M_PI / 10.0);
92  base_node y = x - trans;
93  scalar_type a = gmm::sqrt(gmm::vect_norm2(y));
94  for (int i = 0; i < N; ++i) res[i] = a;
95  break;
96  }
97  }
98  return res;
99 }
100 
101 base_small_vector sol_f(const base_node &x) {
102  int N = x.size();
103  base_small_vector res(N);
104  switch (sol_sing) {
105  case 0 :
106  for (int i = 0; i < N; i++) {
107  res[i] = alph * ( sol_mu * gmm::vect_sp(sol_K.row(i), sol_K.row(i)) )
108  * sin(gmm::vect_sp(sol_K.row(i), x));
109  for (int j = 0; j < N; j++)
110  res[i] += alph * ( (sol_lambda + sol_mu) * sol_K(j,j) * sol_K(j,i))
111  * sin(gmm::vect_sp(sol_K.row(j), x));
112  }
113  break;
114  case 1 :
115  {
116  base_small_vector trans(x.size());
117  gmm::fill(trans, M_PI / 10.0);
118  base_node y = x - trans;
119  scalar_type r = gmm::vect_norm2(y) + 1e-100;
120  scalar_type r2 = r*r;
121  scalar_type tr(0); tr = std::accumulate(y.begin(), y.end(), tr);
122 
123  for (int i = 0; i < N; i++)
124  res[i] = sol_lambda * (y[i]*tr / r2 - 1.0) / r
125  + sol_mu * (y[i]*tr/r2 - N) / r;
126  }
127  break;
128  case 2 :
129  {
130  base_small_vector trans(x.size());
131  gmm::fill(trans, M_PI / 10.0);
132  base_node y = x - trans;
133 
134  scalar_type r = gmm::vect_norm2(y) + 1e-100;
135  scalar_type a = gmm::sqrt(1.0/r);
136  scalar_type b = a*a*a, c = b*b*a;
137  scalar_type tr(0); tr = std::accumulate(y.begin(), y.end(), tr);
138  for (int i = 0; i < N; ++i)
139  res[i] -= b * (sol_lambda + sol_mu * (N+1-3.0/2.0)) * 0.5
140  - c * 3.0 * (sol_lambda + sol_mu) * (y[i]*tr) / 4.0;
141  }
142  break;
143  }
144  return res;
145 }
146 
147 base_matrix sol_sigma(const base_node &x) {
148  int N = x.size();
149  base_matrix res(N,N);
150  switch (sol_sing) {
151  case 0 :
152  for (int i = 0; i < N; i++)
153  for (int j = 0; j <= i; j++) {
154  res(j,i) = res(i,j) = alph * sol_mu *
155  ( sol_K(i,j) * cos(gmm::vect_sp(sol_K.row(i), x))
156  + sol_K(j,i) * cos(gmm::vect_sp(sol_K.row(j), x))
157  );
158  if (i == j)
159  for (int k = 0; k < N; k++)
160  res(i,j) += alph * sol_lambda * sol_K(k,k)
161  * cos(gmm::vect_sp(sol_K.row(k), x));
162  }
163  break;
164  case 1 :
165  {
166  base_small_vector trans(x.size());
167  gmm::fill(trans, M_PI / 10.0);
168  base_node y = x - trans;
169  scalar_type r = gmm::vect_norm2(y) + 1e-100;
170  scalar_type tr(0); tr = std::accumulate(y.begin(), y.end(), tr);
171  for (int i = 0; i < N; i++) {
172  res(i, i) += sol_lambda * tr / r;
173  for (int j = 0; j < N; j++)
174  res(i, j) += sol_mu * (y[i] + y[j]) / r;
175  }
176  }
177 
178  break;
179  case 2 :
180  {
181  base_small_vector trans(x.size());
182  gmm::fill(trans, M_PI / 10.0);
183  base_node y = x - trans;
184  scalar_type r = gmm::vect_norm2(y) + 1e-100;
185  scalar_type a = gmm::sqrt(1.0/r);
186  scalar_type b = a*a*a;
187  scalar_type tr(0); tr = std::accumulate(y.begin(), y.end(), tr);
188  for (int i = 0; i < N; i++) {
189  res(i, i) += sol_lambda * tr * b * 0.5;
190  for (int j = 0; j < N; j++)
191  res(i, j) += sol_mu * b * (y[i] + y[j]) * 0.5;
192  }
193  }
194  }
195 
196  return res;
197 }
198 
199 /*
200  structure for the elastostatic problem
201 */
202 struct elastostatic_problem {
203 
204  enum { DIRICHLET_BOUNDARY_NUM = 0, NEUMANN_BOUNDARY_NUM = 1};
205  getfem::mesh mesh; /* the mesh */
206  getfem::mesh_im mim; /* the integration methods. */
207  getfem::mesh_fem mf_u; /* main mesh_fem, for the elastostatic solution */
208  getfem::mesh_fem mf_mult; /* mesh_fem for the Dirichlet condition. */
209  getfem::mesh_fem mf_rhs; /* mesh_fem for the right hand side (f(x),..) */
210  getfem::mesh_fem mf_p; /* mesh_fem for the pressure for mixed form */
211  scalar_type lambda, mu; /* Lamé coefficients. */
212 
213  scalar_type residual; /* max residual for iterative solvers */
214  bool mixed_pressure, refine;
215  size_type dirichlet_version;
216 
217  std::string datafilename;
218  bgeot::md_param PARAM;
219 
220  void select_boundaries(void);
221  bool solve(plain_vector &U);
222  void init(void);
223  void compute_error(plain_vector &U);
224  elastostatic_problem(void) : mim(mesh),mf_u(mesh), mf_mult(mesh),
225  mf_rhs(mesh),mf_p(mesh) {}
226 };
227 
228 /* Selects the boundaries */
229 
230 void elastostatic_problem::select_boundaries(void) {
231  size_type N = mesh.dim();
232  getfem::mesh_region border_faces;
233  getfem::outer_faces_of_mesh(mesh, border_faces);
234  for (getfem::mr_visitor i(border_faces); !i.finished(); ++i) {
235  base_node un = mesh.normal_of_face_of_convex(i.cv(), i.f());
236  un /= gmm::vect_norm2(un);
237  if (gmm::abs(un[N-1] - 1.0) < 0.5) { // new Neumann face
238  mesh.region(NEUMANN_BOUNDARY_NUM).add(i.cv(), i.f());
239  } else {
240  mesh.region(DIRICHLET_BOUNDARY_NUM).add(i.cv(), i.f());
241  }
242  }
243 }
244 
245 
246 /* Read parameters from the .param file, build the mesh, set finite element
247  * and integration methods and selects the boundaries.
248  */
249 void elastostatic_problem::init(void) {
250  std::string MESH_FILE = PARAM.string_value("MESH_FILE", "Mesh file");
251  std::string FEM_TYPE = PARAM.string_value("FEM_TYPE","FEM name");
252  std::string INTEGRATION = PARAM.string_value("INTEGRATION",
253  "Name of integration method");
254 
255  cout << "MESH_FILE=" << MESH_FILE << "\n";
256  cout << "FEM_TYPE=" << FEM_TYPE << "\n";
257  cout << "INTEGRATION=" << INTEGRATION << "\n";
258 
259 #if GETFEM_PARA_LEVEL > 1
260  double t_init=MPI_Wtime();
261 #endif
262 
263  size_type NX = PARAM.int_value("NX");
264  size_type N = PARAM.int_value("N");
265  std::stringstream filename; filename << MESH_FILE;
266  if ((MESH_FILE.compare(0,11,"structured:") == 0) && NX > 0) {
267  filename << ";NSUBDIV=[" << NX;
268  for (size_type i = 1; i < N; ++i) filename << "," << NX;
269  filename << "];";
270  }
271  getfem::import_mesh(filename.str(), mesh);
272 
273  GMM_ASSERT1(N == mesh.dim(), "The mesh has not the right dimension");
274 
275 #if GETFEM_PARA_LEVEL > 1
276  cout<<"temps creation maillage "<< MPI_Wtime()-t_init<<endl;
277 #endif
278 
279  dirichlet_version
280  = size_type(PARAM.int_value("DIRICHLET_VERSION",
281  "Dirichlet version"));
282  datafilename = PARAM.string_value("ROOTFILENAME","Base name of data files.");
283  scalar_type FT = PARAM.real_value("FT", "parameter for exact solution");
284  residual = PARAM.real_value("RESIDUAL");
285  if (residual == 0.) residual = 1e-10;
286  gmm::resize(sol_K, N, N);
287  for (size_type i = 0; i < N; i++)
288  for (size_type j = 0; j < N; j++)
289  sol_K(i,j) = (i == j) ? FT : -FT;
290 
291  mu = PARAM.real_value("MU", "Lamé coefficient mu");
292  lambda = PARAM.real_value("LAMBDA", "Lamé coefficient lambda");
293  sol_sing = int(PARAM.int_value("SOL_SING", "Optional singular solution"));
294  refine = (PARAM.int_value("REFINE", "Optional refinement") != 0);
295  sol_lambda = lambda; sol_mu = mu;
296  mf_u.set_qdim(dim_type(N));
297  mf_mult.set_qdim(dim_type(N));
298 
299  /* set the finite element on the mf_u */
300  getfem::pfem pf_u =
301  getfem::fem_descriptor(FEM_TYPE);
302  getfem::pintegration_method ppi =
303  getfem::int_method_descriptor(INTEGRATION);
304 
305  mim.set_integration_method(ppi);
306  mf_u.set_finite_element(pf_u);
307 
308  std::string dirichlet_fem_name = PARAM.string_value("DIRICHLET_FEM_TYPE");
309  if (dirichlet_fem_name.size() == 0)
310  mf_mult.set_finite_element(pf_u);
311  else {
312  cout << "DIRICHLET_FEM_TYPE=" << dirichlet_fem_name << "\n";
313  mf_mult.set_finite_element(getfem::fem_descriptor(dirichlet_fem_name));
314  }
315 
316  mixed_pressure =
317  (PARAM.int_value("MIXED_PRESSURE","Mixed version or not.") != 0);
318  if (mixed_pressure) {
319  std::string FEM_TYPE_P = PARAM.string_value("FEM_TYPE_P","FEM name P");
320  mf_p.set_finite_element(getfem::fem_descriptor(FEM_TYPE_P));
321  }
322 
323  /* set the finite element on mf_rhs (same as mf_u is DATA_FEM_TYPE is
324  not used in the .param file */
325  std::string data_fem_name = PARAM.string_value("DATA_FEM_TYPE");
326  if (data_fem_name.size() == 0) {
327  GMM_ASSERT1(pf_u->is_lagrange(), "You are using a non-lagrange FEM. "
328  << "In that case you need to set "
329  << "DATA_FEM_TYPE in the .param file");
330  mf_rhs.set_finite_element(pf_u);
331  } else {
332  mf_rhs.set_finite_element(getfem::fem_descriptor(data_fem_name));
333  }
334 
335  /* set boundary conditions
336  * (Neuman on the upper face, Dirichlet elsewhere) */
337  cout << "Selecting Neumann and Dirichlet boundaries\n";
338  select_boundaries();
339 
340 #if GETFEM_PARA_LEVEL > 1
341 
342 
343  t_init = MPI_Wtime();
344 
345  mf_u.nb_dof(); mf_rhs.nb_dof(); mf_mult.nb_dof();
346 
347  cout<<"enumerate dof time "<< MPI_Wtime()-t_init<<endl;
348 #else
349  double t_init = gmm::uclock_sec();
350  mf_u.nb_dof(); mf_rhs.nb_dof(); mf_mult.nb_dof();
351  cout << "enumerate dof time " << gmm::uclock_sec() - t_init << endl;
352 #endif
353 }
354 
355 /* compute the error with respect to the exact solution */
356 void elastostatic_problem::compute_error(plain_vector &U) {
357  size_type N = mesh.dim();
358  std::vector<scalar_type> V(mf_rhs.nb_basic_dof()*N);
359  getfem::interpolation(mf_u, mf_rhs, U, V);
360  for (size_type i = 0; i < mf_rhs.nb_basic_dof(); ++i) {
361  gmm::add(gmm::scaled(sol_u(mf_rhs.point_of_basic_dof(i)), -1.0),
362  gmm::sub_vector(V, gmm::sub_interval(i*N, N)));
363  }
364 
365 
366 
367  cout.precision(16);
368  mf_rhs.set_qdim(dim_type(N));
369  scalar_type l2 = getfem::asm_L2_norm(mim, mf_rhs, V);
370  scalar_type h1 = getfem::asm_H1_norm(mim, mf_rhs, V);
371 
372  if (getfem::MPI_IS_MASTER())
373  cout << "L2 error = " << l2 << endl
374  << "H1 error = " << h1 << endl
375  << "Linfty error = " << gmm::vect_norminf(V) << endl;
376 
377 /* getfem::vtk_export exp(datafilename + "_err.vtk", */
378 /* PARAM.int_value("VTK_EXPORT")==1); */
379 /* exp.exporting(mf_rhs); */
380 /* exp.write_point_data(mf_rhs, V, "elastostatic_displacement"); */
381 
382  mf_rhs.set_qdim(1);
383 }
384 
385 /**************************************************************************/
386 /* Model. */
387 /**************************************************************************/
388 
389 
390 bool elastostatic_problem::solve(plain_vector &U) {
391 
392  size_type N = mesh.dim();
393 
394  if (mixed_pressure) cout << "Number of dof for P: " << mf_p.nb_dof() << endl;
395  cout << "Number of dof for u: " << mf_u.nb_dof() << endl;
396 
397  getfem::model model;
398 
399  // Main unknown of the problem.
400  model.add_fem_variable("u", mf_u);
401 
402  // Linearized elasticity brick.
403  model.add_initialized_scalar_data("lambda", mixed_pressure ? 0.0 : lambda);
404  model.add_initialized_scalar_data("mu", mu);
406  (model, mim, "u", "lambda", "mu");
407 
408  // Linearized incompressibility condition brick.
409  if (mixed_pressure) {
410  model.add_initialized_scalar_data("incomp_coeff", 1.0/lambda);
411  model.add_fem_variable("p", mf_p); // Adding the pressure as a variable
413  (model, mim, "u", "p", size_type(-1), "incomp_coeff");
414  }
415 
416  // Volumic source term.
417  std::vector<scalar_type> F(mf_rhs.nb_dof()*N);
418  getfem::interpolation_function(mf_rhs, F, sol_f);
419  model.add_initialized_fem_data("VolumicData", mf_rhs, F);
420  getfem::add_source_term_brick(model, mim, "u", "VolumicData");
421 
422  // Neumann condition.
423  gmm::resize(F, mf_rhs.nb_dof()*N*N);
424  getfem::interpolation_function(mf_rhs, F, sol_sigma, NEUMANN_BOUNDARY_NUM);
425  model.add_initialized_fem_data("NeumannData", mf_rhs, F);
427  (model, mim, "u", "NeumannData", NEUMANN_BOUNDARY_NUM);
428 
429  // Dirichlet condition.
430  gmm::resize(F, mf_rhs.nb_dof()*N);
431  getfem::interpolation_function(mf_rhs, F, sol_u);
432  model.add_initialized_fem_data("DirichletData", mf_rhs, F);
434  (model, mim, "u", mf_u, DIRICHLET_BOUNDARY_NUM, "DirichletData");
435 
436  gmm::iteration iter(residual, 1, 40000);
437 #if GETFEM_PARA_LEVEL > 1
438  double t_init=MPI_Wtime();
439 #endif
440  dal::bit_vector cvref;
441 
442  do { // solve with optional refinement
443 
444  cout << "Total number of variables : " << model.nb_dof() << endl;
445 
446  // Defining the volumic source term.
447  size_type nb_dof_rhs = mf_rhs.nb_dof();
448  gmm::resize(F, nb_dof_rhs * N);
449  getfem::interpolation_function(mf_rhs, F, sol_f);
450  gmm::copy(F, model.set_real_variable("VolumicData"));
451 
452  // Defining the Neumann source term.
453  gmm::resize(F, nb_dof_rhs * N * N);
454  getfem::interpolation_function(mf_rhs, F, sol_sigma, NEUMANN_BOUNDARY_NUM);
455  gmm::copy(F, model.set_real_variable("NeumannData"));
456 
457  // Defining the Dirichlet condition value.
458  gmm::resize(F, nb_dof_rhs * N);
459  getfem::interpolation_function(mf_rhs, F, sol_u, DIRICHLET_BOUNDARY_NUM);
460  gmm::copy(F, model.set_real_variable("DirichletData"));
461 
462  iter.init();
463  getfem::standard_solve(model, iter);
464  gmm::resize(U, mf_u.nb_dof());
465  gmm::copy(model.real_variable("u"), U);
466 
467  if (refine) {
468  plain_vector ERR(mesh.convex_index().last_true()+1);
469  getfem::error_estimate(mim, mf_u, U, ERR);
470 
471  cout << "max = " << gmm::vect_norminf(ERR) << endl;
472  // scalar_type threshold = gmm::vect_norminf(ERR) * 0.7;
473  scalar_type threshold = 0.0001, min_ = 1e18;
474  cvref.clear();
475  for (dal::bv_visitor i(mesh.convex_index()); !i.finished(); ++i) {
476  if (ERR[i] > threshold) cvref.add(i);
477  min_ = std::min(min_, ERR[i]);
478  }
479  cout << "min = " << min_ << endl;
480  cout << "Nb elt to be refined : " << cvref.card() << endl;
481 
482  mesh.Bank_refine(cvref);
483  }
484 
485  } while (refine && cvref.card() > 0);
486 
487 #if GETFEM_PARA_LEVEL > 1
488  cout<<"temps standard solve "<< MPI_Wtime()-t_init<<endl;
489 #endif
490 
491 
492  return (iter.converged());
493 }
494 
495 /**************************************************************************/
496 /* main program. */
497 /**************************************************************************/
498 
499 int main(int argc, char *argv[]) {
500 
501  GETFEM_MPI_INIT(argc, argv); // For parallelized version
502 
503  GMM_SET_EXCEPTION_DEBUG; // Exceptions make a memory fault, to debug.
504  FE_ENABLE_EXCEPT; // Enable floating point exception for Nan.
505 
506  // try {
507 
508 
509 
510  elastostatic_problem p;
511  p.PARAM.read_command_line(argc, argv);
512 #if GETFEM_PARA_LEVEL > 1
513  double t_ref=MPI_Wtime();
514 #endif
515  p.init();
516 #if GETFEM_PARA_LEVEL > 1
517  cout << "temps init "<< MPI_Wtime()-t_ref << endl;
518 #endif
519  if (getfem::MPI_IS_MASTER())
520  p.mesh.write_to_file(p.datafilename + ".mesh");
521 
522  plain_vector U;
523 
524 #if GETFEM_PARA_LEVEL > 1
525  t_ref=MPI_Wtime();
526  cout<<"begining resol"<<endl;
527 #endif
528  if (!p.solve(U)) GMM_ASSERT1(false, "Solve has failed");
529 
530 #if GETFEM_PARA_LEVEL > 1
531  cout << "temps Resol "<< MPI_Wtime()-t_ref << endl;
532  t_ref = MPI_Wtime();
533 #endif
534  p.compute_error(U);
535 #if GETFEM_PARA_LEVEL > 1
536  cout << "temps error "<< MPI_Wtime()-t_ref << endl;
537  t_ref = MPI_Wtime();
538 #endif
539 
540  // if (getfem::MPI_IS_MASTER()) { p.mesh.write_to_file("toto.mesh"); }
541 
542  if (p.PARAM.int_value("VTK_EXPORT") && getfem::MPI_IS_MASTER()) {
543  cout << "export to " << p.datafilename + ".vtk" << "..\n";
544  getfem::vtk_export exp(p.datafilename + ".vtk",
545  p.PARAM.int_value("VTK_EXPORT")==1);
546  exp.exporting(p.mf_u);
547  exp.write_point_data(p.mf_u, U, "elastostatic_displacement");
548  cout << "export done, you can view the data file with (for example)\n"
549  "mayavi2 -d " << p.datafilename << ".vtk -f ExtractVectorNorm -f "
550  "WarpVector -m Surface -m Outline\n";
551  }
552 
553  // } GMM_STANDARD_CATCH_ERROR;
554 
555  GETFEM_MPI_FINALIZE;
556 
557 
558 
559  return 0;
560 }
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.
void interpolation_function(mesh_fem &mf_target, const VECT &VV, F &f, mesh_region rg=mesh_region::all_convexes())
interpolation of a function f on mf_target.
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.
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 ( ).
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.
Describe a mesh (collection of convexes (elements) and points).
Definition: getfem_mesh.h:95
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.
Import mesh files from various formats.
Definition of a posteriori error estimates.
size_type APIDECL add_linear_incompressibility(model &md, const mesh_im &mim, const std::string &varname, const std::string &multname_pressure, size_type region=size_type(-1), const std::string &dataexpr_penal_coeff=std::string())
Mixed linear incompressibility condition brick.
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. ...
void import_mesh(const std::string &filename, const std::string &format, mesh &m)
imports a mesh file.
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.
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
Interpolation of fields from a mesh_fem onto another.
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 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...
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
size_type APIDECL add_normal_source_term_brick(model &md, const mesh_im &mim, const std::string &varname, const std::string &dataexpr, size_type region)
Add a source term on the variable varname on a boundary region.
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
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.
Include common gmm files.
Miscelleanous assembly routines for common terms. Use the low-level generic assembly. Prefer the high-level one.
defines and typedefs for namespace getfem
void resize(M &v, size_type m, size_type n)
*/
Definition: gmm_blas.h:231
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