FEM2D_HEAT\
Finite Element Solution of the Heat Equation\
on a Triangulated Region {#fem2d_heat-finite-element-solution-of-the-heat-equation-on-a-triangulated-region align=”center”}
=============================================
FEM2D_HEAT is a C++ program which applies the finite element method
to solve a form of the time-dependent heat equation over an arbitrary
triangulated region.
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.
Normally, the user does not type in this information by hand, but has a
program fill in the nodes, and perhaps another program that constructs
the triangulation. However, in the simplest case, the user might
construct a very crude triangulation by hand, and have
TRIANGULATION_REFINE
refine it to something more reasonable.
For the following ridiculously small example:
10-11-12
|\ |\
| \ | \
6 7 8 9
| \| \
1-2--3--4-5
the node file would be:
0.0 0.0
1.0 0.0
2.0 0.0
3.0 0.0
4.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
3 5 12 4 9 8
12 10 3 11 7 8
The program is set up to handle the time dependent heat equation with a
right hand side function, and nonhomogeneous Dirichlet boundary
conditions. The state variable U(T,X,Y) is then constrained by:
Ut - ( Uxx + Uyy ) + K(x,y,t) * U = F(x,y,t) in the region
U = G(x,y,t) on the boundary
U = H(x,y,t) at initial time TINIT.
To specify the right hand side function F(x,y,t), the linear coefficient
K(x,y,t), the boundary condition function G(x,y,t), and the initial
condition H(x,y,t), the user has to supply a file, perhaps called
myprog.C, containing several functions:
- double rhs ( int node_num, double node_xy[], double time )
evaluates the right hand side forcing term F(x,y,t).
- double k_coef ( int node_num, double node_xy[], double time
) evaluates K(x,y,t);
- double *dirichlet_condition ( int node_num, double
node_xy[], double time ) evaluates G(x,y,t) for all nodes on
the boundary;
- double *initial_condition ( int node_num, double node_xy[],
double time ) evaluates H(x,y,t) for all nodes at the initial
time.
The program is also able to write out a file containing the solution
value at every node. This file may be used to create contour plots of
the solution.
Usage: {#usage align=”center”}
g++ fem2d_heat.cpp myprog.cpp\
mv a.out fem2d_heat\
fem2d_heat prefix
where prefix is the common file prefix:
- “prefix”_nodes.txt, contains the node coordinates.
- “prefix”_elements.txt, contains the indices of nodes that form
elements.
Licensing: {#licensing align=”center”}
The computer code and data files described and made available on this
web page are distributed under the GNU LGPL
license.
Languages: {#languages align=”center”}
FEM2D_HEAT is available in a C++
version and a FORTRAN90
version and a MATLAB
version.
FD2D_HEAT_STEADY,
a C++ program which uses the finite difference method (FDM) to solve the
steady (time independent) heat equation in 2D.
FEM1D_HEAT_STEADY,
a C++ program which uses the finite element method to solve the steady
(time independent) heat equation in 1D.
FEM2D_HEAT_SQUARE,
a C++ library which defines the geometry of a square region, as well as
boundary and initial conditions for a given heat problem, and is called
by FEM2D_HEAT as part of a solution procedure.
STOCHASTIC_HEAT2D,
a C++ program which implements a finite difference method (FDM) for the
steady (time independent) 2D heat equation, with a stochastic heat
diffusivity coefficient, using gnuplot to illustrate the results.
Reference: {#reference align=”center”}
- Hans Rudolf Schwarz,\
Finite Element Methods,\
Academic Press, 1988,\
ISBN: 0126330107,\
LC: TA347.F5.S3313.
- Gilbert Strang, George Fix,\
An Analysis of the Finite Element Method,\
Cambridge, 1973,\
ISBN: 096140888X,\
LC: TA335.S77.
- Olgierd Zienkiewicz,\
The Finite Element Method,\
Sixth Edition,\
Butterworth-Heinemann, 2005,\
ISBN: 0750663200,\
LC: TA640.2.Z54
Source Code: {#source-code align=”center”}
List of Routines: {#list-of-routines align=”center”}
- MAIN is the main routine of FEM2D_HEAT.
- ASSEMBLE_BACKWARD_EULER adjusts the system for the backward
Euler term.
- ASSEMBLE_BOUNDARY modifies the linear system for the boundary
conditions.
- ASSEMBLE_FEM assembles the finite element system for the heat
equation.
- BANDWIDTH determines the bandwidth of the coefficient matrix.
- BASIS_11_T6: one basis at one point for the T6 element.
- CH_CAP capitalizes a single character.
- CH_EQI is true if two characters are equal, disregarding case.
- CH_TO_DIGIT returns the integer value of a base 10 digit.
- DGB_FA performs a LINPACK-style PLU factorization of a DGB
matrix.
- DGB_PRINT_SOME prints some of a DGB matrix.
- DGB_SL solves a system factored by DGB_FA.
- DTABLE_DATA_READ reads the data from a real TABLE file.
- DTABLE_HEADER_READ reads the header from a real TABLE file.
- FILE_COLUMN_COUNT counts the number of columns in the first
line of a file.
- FILE_NAME_INC increments a partially numeric file name.
- FILE_NAME_SPECIFICATION determines the names of the input
files.
- FILE_ROW_COUNT counts the number of row records in a file.
- I4_MAX returns the maximum of two I4’s.
- I4_MIN returns the smaller of two I4’s.
- I4_MODP returns the nonnegative remainder of integer division.
- I4_WRAP forces an integer to lie between given limits by
wrapping.
- I4COL_COMPARE compares columns I and J of an I4COL.
- I4COL_SORT_A ascending sorts the columns of an I4COL.
- I4COL_SWAP swaps two columns of an integer array.
- I4MAT_TRANSPOSE_PRINT_SOME prints some of an I4MAT,
transposed.
- ITABLE_DATA_READ reads data from an integer TABLE file.
- ITABLE_HEADER_READ reads the header from an integer TABLE
file.
- I4VEC_PRINT_SOME prints “some” of an integer vector.
- LVEC_PRINT prints a logical vector.
- POINTS_PLOT plots a pointset.
- QUAD_RULE sets the quadrature rule for assembly.
- R8_ABS returns the absolute value of an R8.
- R8_HUGE returns a “huge” R8.
- R8_NINT returns the nearest integer to an R8.
- R8MAT_TRANSPOSE_PRINT_SOME prints some of an R8MAT,
transposed.
- R8VEC_PRINT_SOME prints “some” of an R8VEC.
- REFERENCE_TO_PHYSICAL_T3 maps reference points to physical
points.
- S_LEN_TRIM returns the length of a string to the last
nonblank.
- S_TO_I4 reads an I4 from a string.
- S_TO_I4VEC reads an I4VEC from a string.
- S_TO_R8 reads an R8 from a string.
- S_TO_R8VEC reads an R8VEC from a string.
- S_WORD_COUNT counts the number of “words” in a string.
- SOLUTION_WRITE writes the solution to a file.
- SORT_HEAP_EXTERNAL externally sorts a list of items into
ascending order.
- TIMESTAMP prints the current YMDHMS date as a time stamp.
- TIMESTRING returns the current YMDHMS date as a string.
- TRIANGLE_AREA_2D computes the area of a triangle in 2D.
- TRIANGULATION_ORDER6_BOUNDARY_NODE indicates nodes on the
boundary.
- TRIANGULATION_ORDER6_PLOT plots a 6-node triangulation of a
set of nodes.
You can go up one level to the C++ source codes.
Last revised on 13 November 2006.