GetFEM  5.4.2
stokes.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 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[]) { GETFEM_MPI_INIT(argc, 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; GETFEM_MPI_FINALIZE; 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  GETFEM_MPI_INIT(argc, argv);
216  GMM_SET_EXCEPTION_DEBUG; // Exceptions make a memory fault, to debug.
217  FE_ENABLE_EXCEPT; // Enable floating point exception for Nan.
218 
219  try {
220  stokes_problem p;
221  p.PARAM.read_command_line(argc, argv);
222  p.init();
223  // p.mesh.write_to_file(p.datafilename + ".mesh");
224  plain_vector U(p.mf_u.nb_dof());
225  if (!p.solve(U)) GMM_ASSERT1(false, "Solve has failed");
226  if (p.PARAM.int_value("VTK_EXPORT")) {
227  cout << "export to " << p.datafilename + ".vtk" << "..\n";
228  getfem::vtk_export exp(p.datafilename + ".vtk",
229  p.PARAM.int_value("VTK_EXPORT")==1);
230  exp.exporting(p.mf_u);
231  exp.write_point_data(p.mf_u, U, "stokes_displacement");
232  cout << "export done, you can view the data file with (for example)\n"
233  "mayavi2 -d " << p.datafilename << ".vtk -f ExtractVectorNorm -f "
234  "WarpVector -m Surface -m Outline\n";
235  }
236  }
237  GMM_STANDARD_CATCH_ERROR;
238  GETFEM_MPI_FINALIZE;
239 
240  return 0;
241 }
getfem_export.h
Export solutions to various formats.
getfem_model_solvers.h
Standard solvers for model bricks.
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::add_linear_incompressibility
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.
Definition: getfem_models.cc:6377
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
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
gmm::rsvector
sparse vector built upon std::vector.
Definition: gmm_def.h:488
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
getfem::interpolation_function
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.
Definition: getfem_interpolation.h:186
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::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