Go to the documentation of this file.
39 using std::endl;
using std::cout;
using std::cerr;
40 using std::ends;
using std::cin;
45 using bgeot::scalar_type;
47 using bgeot::base_matrix;
53 typedef getfem::modeling_standard_sparse_matrix sparse_matrix;
54 typedef getfem::modeling_standard_plain_vector plain_vector;
59 struct stokes_problem {
61 enum { DIRICHLET_BOUNDARY_NUM = 0, NEUMANN_BOUNDARY_NUM = 1};
71 std::string datafilename;
72 bgeot::md_param PARAM;
74 bool solve(plain_vector &U);
76 stokes_problem(
void) : mim(mesh), mf_u(mesh), mf_p(mesh),
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";
97 std::vector<size_type> nsubdiv(N);
98 std::fill(nsubdiv.begin(),nsubdiv.end(),
99 PARAM.int_value(
"NX",
"Nomber of space steps "));
101 PARAM.int_value(
"MESH_NOISED") != 0);
104 bgeot::base_matrix M(N,N);
106 static const char *t[] = {
"LX",
"LY",
"LZ"};
107 M(i,i) = (i<3) ? PARAM.real_value(t[i],t[i]) : 1.0;
109 mesh.transformation(M);
111 datafilename = PARAM.string_value(
"ROOTFILENAME",
"Base name of data files.");
112 residual = PARAM.real_value(
"RESIDUAL");
if (residual == 0.) residual = 1e-10;
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));
120 getfem::pintegration_method ppi =
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(),
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);
137 mf_rhs.set_finite_element(mesh.convex_index(),
143 cout <<
"Selecting Neumann and Dirichlet boundaries\n";
147 assert(it.is_face());
148 base_node un = mesh.normal_of_face_of_convex(it.cv(), it.f());
150 if (gmm::abs(un[N-1] - 1.0) < 1.0E-7) {
151 mesh.region(NEUMANN_BOUNDARY_NUM).add(it.cv(), it.f());
153 mesh.region(DIRICHLET_BOUNDARY_NUM).add(it.cv(), it.f());
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;
169 bool stokes_problem::solve(plain_vector &U) {
179 (
"lambda", plain_vector(1, 0.0));
182 (model, mim,
"u",
"lambda",
"nu");
189 std::vector<scalar_type> F(mf_rhs.nb_dof()*N);
198 (model, mim,
"u", mf_u, DIRICHLET_BOUNDARY_NUM,
"DirichletData");
206 return (iter.converged());
213 int main(
int argc,
char *argv[]) {
215 GETFEM_MPI_INIT(argc, argv);
216 GMM_SET_EXCEPTION_DEBUG;
221 p.PARAM.read_command_line(argc, argv);
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";
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";
237 GMM_STANDARD_CATCH_ERROR;
Export solutions to various formats.
Standard solvers for model bricks.
size_t size_type
used as the common size type in the library
pintegration_method int_method_descriptor(std::string name, bool throw_if_not_found=true)
Get an integration method from its name .
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Describe an integration method linked to a mesh.
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 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...
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.
Describe a finite element method linked to a mesh.
`‘Model’' variables store the variables, the data and the description of a model.
The Iteration object calculates whether the solution has reached the desired accuracy,...
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.
"iterator" class for regions.
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.
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norm2(const V &v)
Euclidean norm of a vector.
std::shared_ptr< const getfem::virtual_fem > pfem
type of pointer on a fem description
container for small vectors of POD (Plain Old Data) types.
sparse vector built upon std::vector.
Describe a mesh (collection of convexes (elements) and points).
std::shared_ptr< const bgeot::geometric_trans > pgeometric_trans
pointer type for a geometric transformation
Include common gmm files.
structure used to hold a set of convexes and/or convex faces.
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.
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 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.
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.
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.
pfem fem_descriptor(const std::string &name)
get a fem descriptor from its string name.
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.
Miscelleanous assembly routines for common terms. Use the low-level generic assembly....