FEMPACK
Finite Element Routines

FEMPACK is a FORTRAN90 library, using double precision arithmetic, which can help in implementing the finite element method.

The emphasis is on simplicity and clarity. Primarily, the 2D planar case is handled, with a choice of low order triangular ("T") and quadrilateral ("Q") elements.

A few routines are included for computing a "sphere grid", that is, a finite element mesh on the surface of a sphere.

Related Data and Programs:

DUNAVANT is a MATLAB library which defines Dunavant rules for quadrature on a triangle.

FEKETE is a FORTRAN90 library of routines for defining a Fekete rule for quadrature or interpolation over a triangle.

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 FORTRAN90 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 an executable FORTRAN90 program which applies the finite element method, with piecewise linear basis functions, to a linear two point boundary value problem;

FEM1D_ADAPTIVE is a FORTRAN90 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 FORTRAN90 program that applies the finite element method to a nonlinear two point boundary value problem in a 1D region.

FEM1D_PMETHOD is a FORTRAN90 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 FORTRAN90 program that solves the time dependent heat equation in the unit square.

FEM2D_POISSON is a FORTRAN90 program for solving Poisson's equation on a square, using the finite element method.

FEMPACK is also available in a C++ version and a MATLAB version.

FREE_FEM_HEAT is a FORTRAN90 program which uses the finite element method and the backward Euler method to solve the time-dependent heat equation on an arbitrary triangulated region in 2D. (The user supplies some information in some problem-dependent subroutines.)

FREE_FEM_POISSON is a FORTRAN90 program for solving Poisson's equation on a triangulated region, using the finite element method.

FREE_FEM_STOKES is a FORTRAN90 program for solving the Stokes equations on a triangulated region, using the finite element method.

HOT_PIPE is a sample problem that can be run with the FEM_50_HEAT MATLAB script.

HOT_POINT is a sample problem that can be run with the FEM_50_HEAT MATLAB script.

MESH_BANDWIDTH is an interactive executable FORTRAN90 program which returns the geometric bandwidth associated with a mesh of elements of any order and in a space of arbitrary dimension.

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.

POLYGONAL_SURFACE is a data directory which contains a description of the data files that can be used to describe a polygonal surface. FEMPACK includes routines for creating a polygonal surface grid of a sphere.

TRIANGULATION_ORDER3 is a directory which contains a description and examples of how an order 3 (linear) triangulation is defined by a node and element file.

TRIANGULATION_ORDER6 is a directory which contains a description and examples of how an order 6 (quadratic) triangulation is defined by a node and element file.

TRIANGULATION_ORDER6_CONTOUR is a MATLAB script which can make contours of a scalar quantity defined on the nodes of an order 6 triangulation.

WANDZURA is a FORTRAN90 library of routines for defining Dunavant rules for quadrature on a triangle.

Reference:

1. Milton Abramowitz, Irene Stegun,
   Handbook of Mathematical Functions,
   National Bureau of Standards, 1964,
   ISBN: 0-486-61272-4,
   LC: QA47.A34.
2. Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
   LINPACK User's Guide,
   SIAM, 1979,
   ISBN13: 978-0-898711-72-1.
3. Vladimir Krylov,
   Approximate Calculation of Integrals,
   Dover, 2006,
   ISBN: 0486445798.
4. Hans Rudolf Schwarz,
   Methode der Finiten Elemente,
   Teubner Studienbuecher, 1980,
   ISBN: 3-519-02349-0.
5. Gilbert Strang, George Fix,
   An Analysis of the Finite Element Method,
   Cambridge, 1973,
   ISBN: 096140888X,
   LC: TA335.S77.
6. Arthur Stroud,
   Approximate Calculation of Multiple Integrals,
   Prentice Hall, 1971,
   ISBN: 0130438936,
   LC: QA311.S85.
7. Arthur Stroud, Don Secrest,
   Gaussian Quadrature Formulas,
   Prentice Hall, 1966,
   LC: QA299.4G3S7.
8. Olgierd Zienkiewicz,
   The Finite Element Method,
   Sixth Edition,
   Butterworth-Heinemann, 2005,
   ISBN: 0750663200,
   TA640.2.Z54
9. Daniel Zwillinger, editor,
   CRC Standard Mathematical Tables and Formulae,
   30th Edition,
   CRC Press, 1996,
   ISBN: 0-8493-2479-3.

Source Code:

• fempack.f90, the source code;
• fempack.csh, commands to compile the source code.

Examples and Tests:

• fempack_prb.f90, the calling program;
• fempack_prb.csh, commands to compile, link and run the calling program;
• fempack_prb_output.txt, the output file.
• fempack_prb_q4.png, a PNG image of a Q4 mesh.
• fempack_prb_t3.png, a PNG image of a T3 mesh.
• fempack_prb_t4.png, a PNG image of a T4 mesh.
• fempack_prb_t6.png, a PNG image of a T6 mesh.
• sphere_q4_nodes.txt, a node file for a Q4 mesh, which is also a polygonal surface, that approximates a sphere.
• sphere_q4_elements.txt, an element file for a Q4 mesh, which is also a polygonal surface, that approximates a sphere.
• sphere_q9_nodes.txt, a node file for a Q9 mesh that approximates a sphere.
• sphere_q9_elements.txt, an element file for a Q9 mesh that approximates a sphere.
• sphere_q16_nodes.txt, a node file for a Q16 mesh that approximates a sphere.
• sphere_q16_elements.txt, an element file for a Q16 mesh that approximates a sphere.
• sphere_t3_nodes.txt, a node file for a T3 mesh, which is also a polygonal surface, that approximates a sphere.
• sphere_t3_elements.txt, an element file for a T3 mesh, which is also a polygonal surface, that approximates a sphere.
• sphere_t6_nodes.txt, a node file for a T6 mesh that approximates a sphere.
• sphere_t6_elements.txt, an element file for a T6 mesh that approximates a sphere.

List of Routines:

• BANDWIDTH_MESH: bandwidth of finite element mesh.
• BANDWIDTH_VAR determines the bandwidth for finite element variables.
• BASIS_11_T3: one basis at one point for the T3 element.
• BASIS_11_T3_TEST verifies BASIS_11_T3.
• BASIS_11_T6: one basis at one point for the T6 element.
• BASIS_11_T6_TEST verifies BASIS_11_T6.
• BASIS_MN_Q4: all bases at N points for a Q4 element.
• BASIS_MN_Q4_TEST verifies BASIS_MN_Q4.
• BASIS_MN_T3: all bases at N points for a T3 element.
• BASIS_MN_T3_TEST verifies BASIS_MN_T3.
• BASIS_MN_T4: all bases at N points for a T4 element.
• BASIS_MN_T4_TEST verifies BASIS_MN_T4.
• BASIS_MN_T6: all bases at N points for a T6 element.
• BASIS_MN_T6_TEST verifies BASIS_MN_T6.
• CH_CAP capitalizes a single character.
• DEGREES_TO_RADIANS converts an angle from degrees to radians.
• DERIVATIVE_AVERAGE_T3 averages derivatives at the nodes of a T3 mesh.
• DERIVATIVE_AVERAGE_T6 averages derivatives at the nodes of a T6 mesh.
• DIV_Q4 estimates the divergence and vorticity of a discrete field.
• DIV_T3 estimates the divergence and vorticity of a discrete field.
• ELEMENT_CODE returns the code for each element.
• ELEMENTS_EPS creates an EPS file image of the elements of a grid.
• GET_UNIT returns a free FORTRAN unit number.
• GRID_ELEMENT returns the element grid associated with any available element.
• GRID_ELEMENT_NUM returns the number of elements in a grid.
• GRID_NODE_NUM returns the number of nodes in a grid.
• GRID_PRINT prints the elements that form a grid.
• GRID_Q4_ELEMENT produces an element grid of 4 node quadrilaterals.
• GRID_Q4_ELEMENT_NUM counts the elements in a grid of 4 node quadrilaterals.
• GRID_Q4_NODE_NUM counts the nodes in a grid of 4 node quadrilaterals.
• GRID_Q8_ELEMENT produces an element grid of 8 node quadrilaterals.
• GRID_Q8_ELEMENT_NUM counts the elements in a grid of 8 node quadrilaterals.
• GRID_Q8_NODE_NUM counts the nodes in a grid of 8 node quadrilaterals.
• GRID_Q9_ELEMENT produces an element grid of 9 node quadrilaterals.
• GRID_Q9_ELEMENT_NUM counts the elements in a grid of 9 node quadrilaterals.
• GRID_Q9_NODE_NUM counts the nodes in a grid of 9 node quadrilaterals.
• GRID_Q12_ELEMENT produces an element grid of 12 node quadrilaterals.
• GRID_Q12_ELEMENT_NUM counts the elements in a grid of 12 node quadrilaterals.
• GRID_Q12_NODE_NUM counts the nodes in a grid of 12 node quadrilaterals.
• GRID_Q16_ELEMENT produces an element grid of 16 node quadrilaterals.
• GRID_Q16_ELEMENT_NUM counts the elements in a grid of 16 node quadrilaterals.
• GRID_Q16_NODE_NUM counts the nodes in a grid of 16 node quadrilaterals.
• GRID_QL_ELEMENT produces an element grid of 6 node quadratics/linears.
• GRID_QL_ELEMENT_NUM counts the elements in a grid of QL quadrilaterals.
• GRID_QL_NODE_NUM counts the nodes in a grid of QL quadrilaterals.
• GRID_SHAPE_2D guesses the shape N1 by N2 of a vector of data.
• GRID_T3_ELEMENT produces an element grid of pairs of 3 node triangles.
• GRID_T3_ELEMENT_NUM counts the elements in a grid of 3 node triangles.
• GRID_T3_NODE_NUM counts the nodes in a grid of 3 node triangles.
• GRID_T4_ELEMENT produces an element grid of pairs of 4 node triangles.
• GRID_T4_ELEMENT_NUM counts the elements in a grid of 4 node triangles.
• GRID_T4_NODE_NUM counts the nodes in a grid of 4 node triangles.
• GRID_T6_ELEMENT produces an element grid of pairs of 6 node triangles.
• GRID_T6_ELEMENT_NUM counts the elements in a grid of 6 node triangles.
• GRID_T6_NODE_NUM counts the nodes in a grid of 6 node triangles.
• GRID_T10_ELEMENT produces an element grid of pairs of 10 node triangles.
• GRID_T10_ELEMENT_NUM counts the elements in a grid of 10 node triangles.
• GRID_T10_NODE_NUM counts the nodes in a grid of 10 node triangles.
• GRID_TEST tests the grid element routines.
• GRID_WIDTH computes the width of a given grid.
• I4_MODP returns the nonnegative remainder of integer division.
• I4_WRAP forces an integer to lie between given limits by wrapping.
• INTERP interpolates a quantity in an element from basis node values.
• INTERP_TEST tests the interpolation property of an element.
• ITABLE_DATA_WRITE writes data to an integer table file.
• ITABLE_HEADER_WRITE writes the header to an integer table file.
• ITABLE_WRITE writes an integer table file.
• I4VEC_PRINT prints an I4VEC.
• LEGENDRE_COM computes abscissas and weights for Gauss-Legendre quadrature.
• LEGENDRE_SET sets abscissas and weights for Gauss-Legendre quadrature.
• MAP returns the interpolation matrix for any available element.
• MAP_TEST tests the map routines.
• MASS_MATRIX_T6 computes the mass matrix, using 6-node triangles.
• NEXT_BOUNDARY_NODE returns the next boundary node in any element.
• NEXT_BOUNDARY_NODE_Q4 returns the next boundary node in a Q4 element.
• NEXT_BOUNDARY_NODE_Q8 returns the next boundary node in a Q8 element.
• NEXT_BOUNDARY_NODE_Q9 returns the next boundary node in a Q9 element.
• NEXT_BOUNDARY_NODE_Q12 returns the next boundary node in a Q12 element.
• NEXT_BOUNDARY_NODE_Q16 returns the next boundary node in a Q16 element.
• NEXT_BOUNDARY_NODE_QL returns the next boundary node in a QL element.
• NEXT_BOUNDARY_NODE_T3 returns the next boundary node in a T3 element.
• NEXT_BOUNDARY_NODE_T4 returns the next boundary node in a T4 element.
• NEXT_BOUNDARY_NODE_T6 returns the next boundary node in a T6 element.
• NEXT_BOUNDARY_NODE_T10 returns the next boundary node in a T10 element.
• NODE_REFERENCE returns the basis nodes for any available element.
• NODE_REFERENCE_Q4 returns the basis nodes for a 4 node quadrilateral.
• NODE_REFERENCE_Q8 returns the basis nodes for an 8 node quadrilateral.
• NODE_REFERENCE_Q9 returns the basis nodes for a 9 node quadrilateral.
• NODE_REFERENCE_Q12 returns the basis nodes for a 12 node quadrilateral.
• NODE_REFERENCE_Q16 returns the basis nodes for a 16 node quadrilateral.
• NODE_REFERENCE_QL returns the basis nodes for a quadratic/linear.
• NODE_REFERENCE_T3 returns the basis nodes for the 3 node triangle.
• NODE_REFERENCE_T4 returns the basis nodes for the 4 node triangle.
• NODE_REFERENCE_T6 returns the basis nodes for a 6 node triangle.
• NODE_REFERENCE_T10 returns the basis nodes for a 10 node triangle.
• NS_T6_VAR_COUNT counts the Navier Stokes variables on a T6 grid.
• NS_T6_VAR_SET sets the Navier Stokes variables on a T6 grid.
• ORDER_CODE returns the order for each element.
• PHYSICAL_TO_REFERENCE_T3 maps physical points to reference points.
• POINTS_PLOT plots a pointset.
• POLY returns the polynomial terms associated with any available element.
• POLY_Q4 returns the monomials associated with a 4 node quadrilateral.
• POLY_Q8 returns the monomials associated with an 8 node quadrilateral.
• POLY_Q9 returns the monomials associated with a 9 node quadrilateral.
• POLY_Q12 returns the monomials associated with a 12 node quadrilateral.
• POLY_Q16 returns the monomials associated with a 16 node quadrilateral.
• POLY_QL returns the monomials for a quadratic/linear quadrilateral.
• POLY_T3 returns the monomials associated with a 3 node triangle.
• POLY_T6 returns the monomials associated with a 6 node triangle.
• POLY_T10 returns the monomials associated with a 10 node triangle.
• R8_HUGE returns the largest legal R8.
• R8_POWER computes the P-th power of R.
• R8_SWAP switches two R8's.
• R8_UNIFORM_01 returns a unit pseudorandom R8.
• R8GE_FA performs a LINPACK style PLU factorization of an R8GE matrix.
• R8GE_INVERSE computes the inverse of a matrix factored by R8GE_FA.
• R8MAT_PRINT prints an R8MAT.
• R8MAT_PRINT_SOME prints some of an R8MAT.
• R8TABLE_DATA_WRITE writes data to an R8TABLE file.
• R8TABLE_HEADER_WRITE writes the header to an R8TABLE file.
• R8TABLE_WRITE writes an R8TABLE file.
• REFERENCE_SAMPLE samples a reference element.
• REFERENCE_TO_PHYSICAL_T3 maps T3 reference points to physical points.
• REFERENCE_TO_PHYSICAL_T6 maps T6 reference points to physical points.
• S_BLANK_DELETE removes blanks from a string, left justifying the remainder.
• S_EQI is a case insensitive comparison of two strings for equality.
• S_L2NORM computes the "big" L2 norm of a scalar function over a region.
• SERENE interpolates data using a Q8 element.
• SHAPE evaluates shape functions for any available reference element.
• SHAPE_Q4 evaluates shape functions for a 4 node reference quadrilateral.
• SHAPE_Q8 evaluates shape functions for an 8 node reference quadrilateral.
• SHAPE_Q9 evaluates shape functions for a 9 node reference quadrilateral.
• SHAPE_Q12 evaluates shape functions for a 12 node reference quadrilateral.
• SHAPE_Q16 evaluates shape functions for a 16 node reference quadrilateral.
• SHAPE_QL evaluates shape functions for a 6 node quadratic/linear.
• SHAPE_T3 evaluates shape functions for a 3 node reference triangle.
• SHAPE_T4 evaluates shape functions for a 4 node reference triangle.
• SHAPE_T6 evaluates shape functions for a 6 node reference triangle.
• SHAPE_T10 evaluates shape functions for a 10 node reference triangle.
• SHAPE_TEST verifies the shape function values at the basis nodes.
• SPHERE_GRID_ELEMENT_NUM returns the number of elements in a sphere grid.
• SPHERE_GRID_NODE_NUM returns the number of nodes in a sphere grid.
• SPHERE_GRID_Q4_ELEMENT produces a Q4 sphere grid.
• SPHERE_GRID_Q4_ELEMENT_NUM counts the elements in a Q4 sphere grid.
• SPHERE_GRID_Q4_NODE_NUM counts nodes in a Q4 sphere grid.
• SPHERE_GRID_Q4_NODE_XYZ produces node coordinates for a Q4 sphere grid.
• SPHERE_GRID_Q9_ELEMENT produces a Q9 sphere grid.
• SPHERE_GRID_Q9_ELEMENT_NUM counts the elements in a Q9 sphere grid.
• SPHERE_GRID_Q9_NODE_NUM counts nodes in a Q9 sphere grid.
• SPHERE_GRID_Q9_NODE_XYZ produces node coordinates for a Q9 sphere grid.
• SPHERE_GRID_Q16_ELEMENT produces a Q16 sphere grid.
• SPHERE_GRID_Q16_ELEMENT_NUM counts the elements in a Q16 sphere grid.
• SPHERE_GRID_Q16_NODE_NUM counts nodes in a Q16 sphere grid.
• SPHERE_GRID_Q16_NODE_XYZ produces node coordinates for a Q16 sphere grid.
• SPHERE_GRID_T3_ELEMENT produces a T3 sphere grid.
• SPHERE_GRID_T3_ELEMENT_NUM counts the elements in a T3 sphere grid.
• SPHERE_GRID_T3_NODE_NUM counts nodes in a T3 sphere grid.
• SPHERE_GRID_T3_NODE_XYZ produces node coordinates for a T3 sphere grid.
• SPHERE_GRID_T6_ELEMENT produces a T6 sphere grid.
• SPHERE_GRID_T6_ELEMENT_NUM counts the elements in a T6 sphere grid.
• SPHERE_GRID_T6_NODE_NUM counts nodes in a T6 sphere grid.
• SPHERE_GRID_T6_NODE_XYZ produces node coordinates for a T6 sphere grid.
• TIMESTAMP prints the current YMDHMS date as a time stamp.
• TIMESTRING writes the current YMDHMS date into a string.
• TRIANGLE_UNIT_SET sets a quadrature rule in the unit triangle.
• TRIANGLE_UNIT_SIZE returns the "size" of a unit triangle quadrature rule.

You can go up one level to the FORTRAN90 source codes.