FEM2D_STOKES is a C++ program which applies the finite element method to solve a form of the steady incompressible Stokes's equations over an arbitrary triangulated region in 2D.
The geometry is entirely external to the program. The user specifies one file of nodal coordinates, and one file that describes the triangles in terms of six node coordinates.
The program makes a default assumption that all boundary velocities correspond to Dirichlet boundary conditions, and that one pressure is specified (for uniqueness of the pressure system). The user can adjust these boundary conditions (and even specify Dirichlet constraints on any variable at any node) by setting the appropriate data in certain user routines.
At the moment, Neumann conditions, if specified, must have a zero right hand side. The machinery to integrate a nonzero Neumann condition has not been set up yet.
The linear system is created, stored, and solved using a form of the LINPACK/LAPACK "general band" storage and versions of LINPACK's DGBFA and DGBSL factorization and solution routines.
The computational region is initially unknown by the program. The user specifies it by preparing a file containing the coordinates of the nodes, and a file containing the indices of nodes that make up triangles that form a triangulation of the region. For the following ridiculously small example:
10-11-12
|\ |\
| \ | \
6 7 8 9
| \ | \
| \| \
1--2--3--4--5
the node file could be:
0.0 0.0
1.0 0.0
2.0 0.0
3.0 0.0
0.0 1.0
1.0 1.0
2.0 1.0
3.0 1.0
0.0 2.0
1.0 2.0
2.0 2.0
and the triangle file would be
1 3 10 2 7 6
12 10 3 11 7 8
3 5 12 4 9 8
The state variables are a velocity vector (U,V)(X,Y) and a scalar pressure P(X,Y). The state variables are constrained by the momentum and continuity equations, which apply inside the region:
- nu * ( Uxx + Uyy ) + dP/dx = U_RHS(x,y)
- nu * ( Vxx + Vyy ) + dP/dy = V_RHS(x,y)
dU/dx + dV/dy = P_RHS(x,y)
where, typically, the right hand side functions are zero. However,
the user is free to assign nonzero values to these functions
through a user routine.
At every point on the boundary of the region, the program assumes that both components of the velocity are specified.
U(node) = U_BC(node)
V(node) = V_BC(node)
This is known as a "Dirichlet boundary condition". The right hand side
of the boundary condition is left unspecified until the user
routine is called. If a wall is intended, then the appropriate
condition has U_BC and V_BC zero. An inlet might have a
particular flow profile function used to assign nonzero values.
At one point in the region, the program assumes that the value of the pressure is specified.
P(node) = P_BC(node)
Such a condition must be supplied; otherwise the pressure cannot
be uniquely determined, since it is essentially a potential function,
and so is unique only "up to a constant". Note that the program allows
the user to specify pressure conditions anyway, and these can be
of Dirichlet or Neumann type. In general, however, this is not
a physically or mathematically or computationally good thing to do!
The user routine boundary_type can be used to modify the type of the boundary conditions associated with a degree of freedom at a boundary node - or even to add constraints to variables associated with nodes in the interior.
One choice that the user can make is to reset certain boundary conditions to be of Neumann type:
dU/dn(node) = U_BC(node)
dV/dn(node) = V_BC(node)
The right hand side of the boundary condition is left unspecified
until the user routine is called. As mentioned earlier, the program
cannot currently handle Neumann conditions with nonzero right hand side.
(A nonzero value is simply ignored, but won't actually cause the
program to fail.)
We use linear finite elements for the pressure function, and to generate these, we only need the nodes that are the vertices of the triangles. (We will call these vertices "pressure nodes.") Because quadratic basis functions are to be used for the velocity, however, each triangle will also have three extra midside nodes available for that.
We now assume that the unknown velocity component functions U(x,y) and V(x,y) can be represented as linear combinations of the quadratic basis functions associated with each node, and that the scalar pressure P(x,y) can be represented as a linear combination of the linear basis functions associated with each pressure node.
For every node, we can determine a quadratic velocity basis function PSI(I)(x,y). For every pressure node I, we can determine a linear basis function PHI(I)(x,y). If we assume that our solutions are linear combinations of these basis functions, then we need to solve for the coefficients.
To do so, we apply the Galerkin-Petrov method. For each pressure node, and its corresponding basis function PHI(I), we generate a copy of the continuity equation, multiplied by that basis function, and integrated over the region:
Integral ( Ux(x,y) + Vy(x,y) ) * PHI(I)(x,y) dx dy = Integral ( P_RHS(x,y) * PHI(I)(x,y) dx dy )
Similarly, for each node and its corresponding velocity basis function PSI(I), we generate two copies of the momentum equation, one for each component. We multiply by PSI(I), integrate over the region, and use integration by parts to lower the order of differentiation. This gives us:
Integral nu * ( Ux(x,y) * PSIx(I)(x,y) + Uy(x,y) * PSIy(I)(x,y) ) + Px(x,y) * PSI(I)(x,y) dx dy = Integral ( U_RHS(x,y) * PSI(I)(x,y) dx dy )
Integral nu * ( Vx(x,y) * PSIx(I)(x,y) + Vy(x,y) * PSIy(I)(x,y) ) + Py(x,y) * PSI(I)(x,y) dx dy = Integral ( V_RHS(x,y) * PSI(I)(x,y) dx dy )
After adjusting for the boundary conditions, the set of all such equations yields a linear system for the coefficients of the finite element representation of the solution.
The program requires the user to supply the following routines:
The default boundary condition types are passed to the user, along with other information. The user modifies any data as necessary, and returns it. This is done by a user-supplied routine of the form
void boundary_type ( int node_num, double node_xy[], bool node_boundary[], int node_type[], double node_u_condition[], node_v_condition[], node_p_condition[] )
The right hand sides of any Dirichlet boundary conditions are determined by a user-supplied routine of the form
void dirichlet_condition ( int node_num, double node_xy[], double u_bc[], double v_bc[], double p_bc[] )
The right hand sides of any Neumann boundary conditions are determined by a user-supplied routine of the form
void neumann_condition ( int node_num, double node_xy[], double u_bc[], double v_bc[], double p_bc[] )
The right hand sides (or "source terms") of the state equations are determined by a user-supplied routine of the form
void rhs ( int node_num, double node_xy[], double u_rhs[], double v_rhs[], double p_rhs[] )
The program writes out various node, triangle, pressure and velocity data files that can be used to create plots of the geometry and the solution.
Graphics files created include:
Data files created include:
To run the program, the user must first write a file, called, perhaps, myprog.C, containing routines defining certain data, compile this file with fem2d_stokes.C, and create an executable. It would be sensible to rename the executable to fem2d_stokes.
g++ fem2d_stokes.C myprog.C
mv a.out fem2d_stokes
Once the executable is created, it can be run as follows:
fem2d_stokes nodes.txt triangles.txtwhere
The computer code and data files described and made available on this web page are distributed under the GNU LGPL license.
FEM2D_STOKES is available in a C++ version and a FORTRAN90 version and a MATLAB version.
FEM2D_STOKES_CAVITY, a C++ library which contains the user-supplied routines necessary to run fem2d_stokes on the "cavity" problem.
FEM2D_STOKES_CHANNEL. a C++ library which contains the user-supplied routines necessary to run fem2d_stokes on the "channel" problem.
FEM2D_STOKES_INOUT. a C++ library which contains the user-supplied routines necessary to run fem2d_stokes on the "inout" problem.
You can go up one level to the C++ source codes.