FREE_FEM_STOKES is a MATLAB 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.
FFS_SPARSE is a version of this program which uses MATLAB's sparse matrix storage, factorization and solution techniques.
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 MATLAB M file of the form:
function [ node_u_condition, node_v_condition, node_p_condition ] = boundary_type ( node_num, node_xy, node_boundary, node_type, node_u_condition, node_v_condition, node_p_condition )
The value of the kinematic viscosity is determined by a user-supplied MATLAB M file of the form
function nu = constants ( 'DUMMY' ).
The right hand side of any Dirichlet boundary conditions are determined by a user-supplied MATLAB M file of the form
function [ u_bc, v_bc, p_bc ] = dirichlet_condition ( n, xy ).
The right hand side of any Neumann boundary conditions are determined by a user-supplied MATLAB M file of the form
function [ u_bc, v_bc, p_bc ] = neumann_condition ( n, xy ).
The right hand side source term functions are determined by a user-supplied MATLAB M file of the form
function [ u_rhs, v_rhs, p_rhs ] = rhs ( m, xy ).
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 write the four user files described above, make all the files associated with free_fem_stokes available in the same directory, or in the user's MATLAB path, and supplying the names of the node and triangle files to the main program:
BUMP is an executable FORTRAN90 program which solves a fluid flow problem in a channel including a bump which obstructs and redirects the flow.
CHANNEL is an executable FORTRAN90 program which solves a fluid flow problem in a channel.
DIRECTION_ARROWS is a MATLAB program that can be used to plot the vector direction field.
DIRECTION_ARROWS_GRID is a MATLAB program which reads files of node and velocity data, and, using interpolation, creates a velocity direction plot with arrows place on a uniform grid of the user's specification.
FEM is a data directory which contains a description of the data files that can be used to describe a finite element model.
FEM_50 is a MATLAB finite element program in just 50 lines of code.
FEM_50_HEAT is a modified version of FEM_50 suitable for solving the heat equation.
FEM_BASIS_T3_DISPLAY is a MATLAB program which displays a basis function associated with a linear triangle ("T3") mesh.
FEM_BASIS_T6_DISPLAY is a MATLAB program which reads a quadratic triangle mesh and displays any associated basis function.
FEM_IO is a library of MATLAB routines for reading or writing the node, element and data files that define a finite element model.
FEM_SAMPLE is a MATLAB library of routines for evaluating a finite element function defined on an order 3 or order 6 triangulation.
FEM_TO_TEC is a MATLAB program that can convert an FEM model into a TEC graphics file.
FEM1D is a MATLAB program which applies the finite element method, with piecewise linear basis functions, to a linear two point boundary value problem;
FEM1D_ADAPTIVE is a MATLAB program that applies the finite element method to a linear two point boundary value problem in a 1D region, using adaptive refinement to improve the solution.
FEM1D_NONLINEAR is a MATLAB program that applies the finite element method to a nonlinear two point boundary value problem in a 1D region.
FEM1D_PMETHOD is a MATLAB program that applies the p-method version of the finite element method to a linear two point boundary value problem in a 1D region.
FEM2D_HEAT is a MATLAB program that solves the time dependent heat equation in the unit square.
FEM2D_POISSON is a MATLAB program for solving Poisson's equation on a square, using the finite element method.
FEMPACK is a MATLAB library of routines for finite element calculations.
FFNS_SPARSE is a MATLAB program for solving the steady incompressible Navier Stokes equations on an arbitrary triangulated region, using the finite element method and MATLAB's sparse facility.
FFS_SPARSE is a MATLAB version of FREE_FEM_STOKES which uses MATLAB's sparse matrix storage, factorization and solution techniques.
FLOW3 is an executable FORTRAN90 code for solving steady incompressible Navier Stokes equations in 2D using the finite element method.
FLOW5 is an executable FORTRAN90 code for solving steady incompressible Navier Stokes equations in 2D using the finite element method.
FLOW7 is an executable FORTRAN90 code for solving steady incompressible Navier Stokes equations in 2D using the finite element method.
FREE_FEM_HEAT is a MATLAB program that solves the time dependent heat equation in an arbitrary triangulated 2D region.
FREE_FEM_NAVIER_STOKES is a MATLAB program for solving the steady incompressible Navier Stokes equations on an arbitrary triangulated region, using the finite element method.
FREE_FEM_POISSON is a MATLAB program for solving Poisson's equation on a triangulated region, using the finite element method.
FREE_FEM_STOKES is also available in a C++ version and a FORTRAN90 version.
HCELL is a FORTRAN77 program which computes the pressure and velocity for a Navier Stokes flow in an "H"-shaped region.
HOT_PIPE is a sample problem that can be run with FEM_50_HEAT.
HOT_POINT is a sample problem that can be run with FEM_50_HEAT.
INOUT is a FORTRAN77 program which computes the pressure and velocity for a Navier Stokes flow in a square region with an inlet and an outlet.
MHD_CONTROL is an executable FORTRAN90 program which tries to control the evolution of an MHD system so that a particular state is achieved.
MHD_FLOW is an executable FORTRAN90 program for the evolution of an MHD system.
NAST2D is an executable C++ program which uses the finite volume method to set up and solve the 2D incompressible Navier Stokes equations with heat.
NAST2D_F90 is an executable FORTRAN90 program which uses the finite volume method to set up and solve the 2D incompressible Navier Stokes equations with heat.
NSASM is a C library of routines, intended to be used with a MATLAB calling program, and which set up the sparse matrix needed for a Newton iteration to solve a finite element formulation of the steady incompressible 2D Navier Stokes equations.
PLTMG_SINGLE is a FORTRAN77 library of routines for solving elliptic partial differential equations using the finite element method with piecewise linear triangles and the multigrid approach.
PLOT_POINTS is an executable FORTRAN90 program that can make a plot of the nodes that define the region.
TABLE is a format used to store the input and output files used by the program.
TABLE_IO is a MATLAB library of routines needed to read the data files.
TCELL is an executable FORTRAN77 program which computes the pressure and velocity for a Navier Stokes flow in a "T"-shaped region.
TRIANGULATION_ORDER3 is a data directory that describes the format for an order 3 triangulation and includes examples.
TRIANGULATION_ORDER3_CONTOUR is a MATLAB program that can make contour plots of the pressure field computed by this program.
TRIANGULATION_ORDER6 is a data directory that describes the format for an order 6 triangulation and includes examples.
TRIANGULATION_PLOT is a MATLAB program which can be used to plot the elements that triangulate the region.
VECTOR_PLOT is an executable FORTRAN90 program that can plot the velocity field and the velocity direction field.
VECTOR_STREAM_GRID is a MATLAB program which reads node and vector data from a file, computes an interpolatory function, evaluates on a uniform grid of points specified by the user, and displays a streamline plot of the vector field.
VELOCITY_ARROWS is a MATLAB program that can plot the velocity field.
VELOCITY_ARROWS_GRID is a MATLAB program which reads files of node and velocity data, and, using interpolation, creates a vector plot with arrows place on a uniform grid of the user's specification.
A GZIP'ed TAR file of the contents of this directory is available. This is only done as a convenience for users who want ALL the files, and don't want to download them individually. This is not a convenience for me, so don't be surprised if the tar file is somewhat out of date.
You can go up one level to the MATLAB source codes.