jburkardt

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.

Related Programs: {#related-programs align=”center”}

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”}

1. Hans Rudolf Schwarz,\ Finite Element Methods,\ Academic Press, 1988,\ ISBN: 0126330107,\ LC: TA347.F5.S3313.
2. Gilbert Strang, George Fix,\ An Analysis of the Finite Element Method,\ Cambridge, 1973,\ ISBN: 096140888X,\ LC: TA335.S77.
3. Olgierd Zienkiewicz,\ The Finite Element Method,\ Sixth Edition,\ Butterworth-Heinemann, 2005,\ ISBN: 0750663200,\ LC: TA640.2.Z54

Source Code: {#source-code align=”center”}

• fem2d_heat.cpp, the source code;

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.