GetFEM++  5.3
stokes.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 stokes.cc
24  @brief Solve the stokes problem (incompressible viscuous fluid).
25 
26  This program is used to check that getfem++ is working. This is also
27  a good example of use of GetFEM++.
28 
29  This file is almost identical to @link elastostatic.cc
30  tests/elastostatic.cc@endlink, except than a
31  linear incompressibility brick is inserted.
32 */
33 
34 #include "getfem/getfem_assembling.h" /* import assembly methods (and norms comp.) */
35 #include "getfem/getfem_export.h" /* export functions (save solution in a file) */
38 #include "gmm/gmm.h"
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::scalar_type; /* = double */
46 using bgeot::size_type; /* = unsigned long */
47 using bgeot::base_matrix; /* small dense matrix. */
48 
49 /* definition of some matrix/vector types. These ones are built
50  * using the predefined types in Gmm++
51  */
53 typedef getfem::modeling_standard_sparse_matrix sparse_matrix;
54 typedef getfem::modeling_standard_plain_vector plain_vector;
55 
56 /*
57  * structure for the stokes problem
58  */
59 struct stokes_problem {
60 
61  enum { DIRICHLET_BOUNDARY_NUM = 0, NEUMANN_BOUNDARY_NUM = 1};
62  getfem::mesh mesh; /* the mesh */
63  getfem::mesh_im mim; /* integration methods. */
64  getfem::mesh_fem mf_u; /* main mesh_fem, for the velocity */
65  getfem::mesh_fem mf_p; /* mesh_fem for the pressure */
66  getfem::mesh_fem mf_rhs; /* mesh_fem for the right hand side (f(x),..) */
67  scalar_type nu; /* Lamé coefficients. */
68 
69  scalar_type residual; /* max residual for the iterative solvers */
70 
71  std::string datafilename;
72  bgeot::md_param PARAM;
73 
74  bool solve(plain_vector &U);
75  void init(void);
76  stokes_problem(void) : mim(mesh), mf_u(mesh), mf_p(mesh),
77  mf_rhs(mesh) {}
78 };
79 
80 /* Read parameters from the .param file, build the mesh, set finite element
81  * and integration methods and selects the boundaries.
82  */
83 void stokes_problem::init(void) {
84  std::string MESH_TYPE = PARAM.string_value("MESH_TYPE","Mesh type ");
85  std::string FEM_TYPE = PARAM.string_value("FEM_TYPE","FEM name");
86  std::string FEM_TYPE_P = PARAM.string_value("FEM_TYPE_P","FEM name P");
87  std::string INTEGRATION = PARAM.string_value("INTEGRATION",
88  "Name of integration method");
89  cout << "MESH_TYPE=" << MESH_TYPE << "\n";
90  cout << "FEM_TYPE=" << FEM_TYPE << "\n";
91  cout << "INTEGRATION=" << INTEGRATION << "\n";
92 
93  /* First step : build the mesh */
96  size_type N = pgt->dim();
97  std::vector<size_type> nsubdiv(N);
98  std::fill(nsubdiv.begin(),nsubdiv.end(),
99  PARAM.int_value("NX", "Nomber of space steps "));
100  getfem::regular_unit_mesh(mesh, nsubdiv, pgt,
101  PARAM.int_value("MESH_NOISED") != 0);
102 
103  /* scale the unit mesh to [LX,LY,..] and incline it */
104  bgeot::base_matrix M(N,N);
105  for (size_type i=0; i < N; ++i) {
106  static const char *t[] = {"LX","LY","LZ"};
107  M(i,i) = (i<3) ? PARAM.real_value(t[i],t[i]) : 1.0;
108  }
109  mesh.transformation(M);
110 
111  datafilename = PARAM.string_value("ROOTFILENAME","Base name of data files.");
112  residual = PARAM.real_value("RESIDUAL"); if (residual == 0.) residual = 1e-10;
113 
114  nu = PARAM.real_value("NU", "Viscosité); mf_u.set_qdim(bgeot::dim_type(N)); /* set the finite element on the mf_u */ getfem::pfem pf_u = getfem::fem_descriptor(FEM_TYPE); getfem::pintegration_method ppi = getfem::int_method_descriptor(INTEGRATION); mim.set_integration_method(mesh.convex_index(), ppi); mf_u.set_finite_element(mesh.convex_index(), pf_u); mf_p.set_finite_element(mesh.convex_index(), getfem::fem_descriptor(FEM_TYPE_P)); /* set the finite element on mf_rhs (same as mf_u is DATA_FEM_TYPE is not used in the .param file */ std::string data_fem_name = PARAM.string_value("DATA_FEM_TYPE"); if (data_fem_name.size() == 0) { GMM_ASSERT1(pf_u->is_lagrange(), "You are using a non-lagrange FEM " << data_fem_name << ". In that case you need to set " << "DATA_FEM_TYPE in the .param file"); mf_rhs.set_finite_element(mesh.convex_index(), pf_u); } else { mf_rhs.set_finite_element(mesh.convex_index(), getfem::fem_descriptor(data_fem_name)); } /* set boundary conditions * (Neuman on the upper face, Dirichlet elsewhere) */ cout << "Selecting Neumann and Dirichlet boundaries\n"; getfem::mesh_region border_faces; getfem::outer_faces_of_mesh(mesh, border_faces); for (getfem::mr_visitor it(border_faces); !it.finished(); ++it) { assert(it.is_face()); base_node un = mesh.normal_of_face_of_convex(it.cv(), it.f()); un /= gmm::vect_norm2(un); if (gmm::abs(un[N-1] - 1.0) < 1.0E-7) { // new Neumann face mesh.region(NEUMANN_BOUNDARY_NUM).add(it.cv(), it.f()); } else { mesh.region(DIRICHLET_BOUNDARY_NUM).add(it.cv(), it.f()); } } } /**************************************************************************/ /* Model. */ /**************************************************************************/ base_small_vector sol_f(const base_small_vector &P) { base_small_vector res(P.size()); res[P.size()-1] = -1.0; return res; } bool stokes_problem::solve(plain_vector &U) { size_type N = mesh.dim(); getfem::model model; // Main unknown of the problem. model.add_fem_variable("u", mf_u); // Linearized elasticity brick. model.add_initialized_fixed_size_data ("lambda", plain_vector(1, 0.0)); model.add_initialized_fixed_size_data("nu", plain_vector(1, nu)); getfem::add_isotropic_linearized_elasticity_brick (model, mim, "u", "lambda", "nu"); // Linearized incompressibility condition brick. model.add_fem_variable("p", mf_p); // Adding the pressure as a variable add_linear_incompressibility(model, mim, "u", "p"); // Volumic source term. std::vector<scalar_type> F(mf_rhs.nb_dof()*N); getfem::interpolation_function(mf_rhs, F, sol_f); model.add_initialized_fem_data("VolumicData", mf_rhs, F); getfem::add_source_term_brick(model, mim, "u", "VolumicData"); // Dirichlet condition. gmm::clear(F); model.add_initialized_fem_data("DirichletData", mf_rhs, F); getfem::add_Dirichlet_condition_with_multipliers (model, mim, "u", mf_u, DIRICHLET_BOUNDARY_NUM, "DirichletData"); gmm::iteration iter(residual, 1, 40000); getfem::standard_solve(model, iter); // Solution extraction gmm::copy(model.real_variable("u"), U); return (iter.converged()); } /**************************************************************************/ /* main program. */ /**************************************************************************/ int main(int argc, char *argv[]) { GMM_SET_EXCEPTION_DEBUG; // Exceptions make a memory fault, to debug. FE_ENABLE_EXCEPT; // Enable floating point exception for Nan. try { stokes_problem p; p.PARAM.read_command_line(argc, argv); p.init(); // p.mesh.write_to_file(p.datafilename + ".mesh"); plain_vector U(p.mf_u.nb_dof()); if (!p.solve(U)) GMM_ASSERT1(false, "Solve has failed"); if (p.PARAM.int_value("VTK_EXPORT")) { cout << "export to " << p.datafilename + ".vtk" << "..\n"; getfem::vtk_export exp(p.datafilename + ".vtk", p.PARAM.int_value("VTK_EXPORT")==1); exp.exporting(p.mf_u); exp.write_point_data(p.mf_u, U, "stokes_displacement"); cout << "export done, you can view the data file with (for example)\n" "mayavi2 -d " << p.datafilename << ".vtk -f ExtractVectorNorm -f " "WarpVector -m Surface -m Outline\n"; } } GMM_STANDARD_CATCH_ERROR; return 0; } ");
115  mf_u.set_qdim(bgeot::dim_type(N));
116 
117  /* set the finite element on the mf_u */
118  getfem::pfem pf_u =
119  getfem::fem_descriptor(FEM_TYPE);
120  getfem::pintegration_method ppi =
121  getfem::int_method_descriptor(INTEGRATION);
122 
123  mim.set_integration_method(mesh.convex_index(), ppi);
124  mf_u.set_finite_element(mesh.convex_index(), pf_u);
125  mf_p.set_finite_element(mesh.convex_index(),
126  getfem::fem_descriptor(FEM_TYPE_P));
127 
128  /* set the finite element on mf_rhs (same as mf_u is DATA_FEM_TYPE is
129  not used in the .param file */
130  std::string data_fem_name = PARAM.string_value("DATA_FEM_TYPE");
131  if (data_fem_name.size() == 0) {
132  GMM_ASSERT1(pf_u->is_lagrange(), "You are using a non-lagrange FEM "
133  << data_fem_name << ". In that case you need to set "
134  << "DATA_FEM_TYPE in the .param file");
135  mf_rhs.set_finite_element(mesh.convex_index(), pf_u);
136  } else {
137  mf_rhs.set_finite_element(mesh.convex_index(),
138  getfem::fem_descriptor(data_fem_name));
139  }
140 
141  /* set boundary conditions
142  * (Neuman on the upper face, Dirichlet elsewhere) */
143  cout << "Selecting Neumann and Dirichlet boundaries\n";
144  getfem::mesh_region border_faces;
145  getfem::outer_faces_of_mesh(mesh, border_faces);
146  for (getfem::mr_visitor it(border_faces); !it.finished(); ++it) {
147  assert(it.is_face());
148  base_node un = mesh.normal_of_face_of_convex(it.cv(), it.f());
149  un /= gmm::vect_norm2(un);
150  if (gmm::abs(un[N-1] - 1.0) < 1.0E-7) { // new Neumann face
151  mesh.region(NEUMANN_BOUNDARY_NUM).add(it.cv(), it.f());
152  } else {
153  mesh.region(DIRICHLET_BOUNDARY_NUM).add(it.cv(), it.f());
154  }
155  }
156 }
157 
158 /**************************************************************************/
159 /* Model. */
160 /**************************************************************************/
161 
162 base_small_vector sol_f(const base_small_vector &P) {
163  base_small_vector res(P.size());
164  res[P.size()-1] = -1.0;
165  return res;
166 }
167 
168 
169 bool stokes_problem::solve(plain_vector &U) {
170  size_type N = mesh.dim();
171 
172  getfem::model model;
173 
174  // Main unknown of the problem.
175  model.add_fem_variable("u", mf_u);
176 
177  // Linearized elasticity brick.
179  ("lambda", plain_vector(1, 0.0));
180  model.add_initialized_fixed_size_data("nu", plain_vector(1, nu));
182  (model, mim, "u", "lambda", "nu");
183 
184  // Linearized incompressibility condition brick.
185  model.add_fem_variable("p", mf_p); // Adding the pressure as a variable
186  add_linear_incompressibility(model, mim, "u", "p");
187 
188  // Volumic source term.
189  std::vector<scalar_type> F(mf_rhs.nb_dof()*N);
190  getfem::interpolation_function(mf_rhs, F, sol_f);
191  model.add_initialized_fem_data("VolumicData", mf_rhs, F);
192  getfem::add_source_term_brick(model, mim, "u", "VolumicData");
193 
194  // Dirichlet condition.
195  gmm::clear(F);
196  model.add_initialized_fem_data("DirichletData", mf_rhs, F);
198  (model, mim, "u", mf_u, DIRICHLET_BOUNDARY_NUM, "DirichletData");
199 
200  gmm::iteration iter(residual, 1, 40000);
201  getfem::standard_solve(model, iter);
202 
203  // Solution extraction
204  gmm::copy(model.real_variable("u"), U);
205 
206  return (iter.converged());
207 }
208 
209 /**************************************************************************/
210 /* main program. */
211 /**************************************************************************/
212 
213 int main(int argc, char *argv[]) {
214 
215  GMM_SET_EXCEPTION_DEBUG; // Exceptions make a memory fault, to debug.
216  FE_ENABLE_EXCEPT; // Enable floating point exception for Nan.
217 
218  try {
219  stokes_problem p;
220  p.PARAM.read_command_line(argc, argv);
221  p.init();
222  // p.mesh.write_to_file(p.datafilename + ".mesh");
223  plain_vector U(p.mf_u.nb_dof());
224  if (!p.solve(U)) GMM_ASSERT1(false, "Solve has failed");
225  if (p.PARAM.int_value("VTK_EXPORT")) {
226  cout << "export to " << p.datafilename + ".vtk" << "..\n";
227  getfem::vtk_export exp(p.datafilename + ".vtk",
228  p.PARAM.int_value("VTK_EXPORT")==1);
229  exp.exporting(p.mf_u);
230  exp.write_point_data(p.mf_u, U, "stokes_displacement");
231  cout << "export done, you can view the data file with (for example)\n"
232  "mayavi2 -d " << p.datafilename << ".vtk -f ExtractVectorNorm -f "
233  "WarpVector -m Surface -m Outline\n";
234  }
235  }
236  GMM_STANDARD_CATCH_ERROR;
237 
238  return 0;
239 }
structure used to hold a set of convexes and/or convex faces.
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 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.
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.
pgeometric_trans geometric_trans_descriptor(std::string name)
Get the geometric transformation from its string name.
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. ...
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
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.
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:59
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_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
Include common gmm files.
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