Geophysical Finite Element Simulation Tool

GeoFEST (Geophysical Finite Element Simulation Tool) is a two- and three-dimensional finite element software package for the modeling of solid stress and strain in geophysical and other continuum domain applications. The program source is written in C, and consists of approximately 7000 lines of source code. The program is targeted to be compiled and run on UNIX systems, and is running on diverse UNIX derivatives including LINUX, HPUX, and SunOS. The present code is descended from earlier generation FORTAN code running under VAX VMS (VISELAS) and an earlier UNIX C code (VISCO). The program uses input and output in the form of formatted plain text files; the data formats can be adapted to accommodate visualization and graphically oriented i/o. The computational engine of the program employs Crout factorization for the direct inversion of the finite element matrices as well as congugate gradient for an iterative solution option. The physics models supported by the code include isotropic linear elasticity and both Newtonian and power-law viscoelasticity, via implicit/explicit quasi-static time stepping. In addition to triangular, quadrilateral, tetrahedral and hexahedral continuum elements, the program supports split-node faulting, body forces and surface tractions. Capabilities under development include frictional faulting and buoyancy driving. Scientific applications of the code include the modeling of static and transient co- and post-seismic earth deformation, Earth response to glacial, atmospheric and hydrological loading, and other scenarios involving the bulk deformation of geologic media.

Download the latest version (4.3)

Click here to download the complete GeoFEST installation package, with validation test problems and makefile.

User's Guide


In order to simulate viscoelastic stress and flow in a realistic model of the earth's crust and upper mantle, the modeling technique must be able to accommodate a vertically layered structure with imbedded faults which may cut ar arbitrary angles. Stress and displacement features will vary most rapidly near the faults and particularly near fault-terminations. These features argue for fully three-dimensional finite element modeling in the time domain. Two-dimensional modeling, semianalytical techniques, finite difference and semispectral methods either cannot model significant features or geometry of interest, or require gross oversampling in regions of little interest, leading to impossible computational requirements.

Finite element modeling in three dimensions allows faithful modeling of complex faulting geometry, inhomogeneous materials, realistic viscous flow, and a wide variety of fault slip models and boundary conditions. While there are particular problems which are more efficiently expressed using alternative approaches, finite elements represent the most generally useful method for inhomogeneous elastostatic and viscoelastic problems. Because finite elements conform to (nearly) any surface geometry and support wide variations in mesh density, solutions may be made arbitrarily accurate with high computational efficiency. This flexibility comes with a price tag for the user or the toolbuilder, that of generating and adapting the mesh of elements upon which the solution is computed. When such generation tools are primitive, researchers may spend substantially more time creating a mesh than solving their problem and interpreting the results. Therefore we describe automated tools for creating and adapting the mesh, including using an initial coarse mesh to generate a solution, whose computable error characteristics inform a further mesh generation cycle and produce an efficient and accurate solution.

We describe the supported features of the GeoFEST code as of this writing in the next section. The precise description of the codes' operation and mathematical specification follows in the Theory of Operation section. The entries in input and output files are listed next, enabling a user to set up an input file with a text editor, and run Visco on that basis. Then the existing tools for processing the input and output are described, allowing a more automated generation of an input mesh, and visualization of the output. Finally some examples are given.


The primary quantity computed by GeoFEST is the displacement at each point in a domain. The stress tensor is also computed as a necessary byproduct. The computational domain represents a region of the earth's crust and possibly underlying mantle. It is typically a square or rectangular domain in map view, with a flat upper free surface and constant depth, but the domain may deviate from this. The only requirement is that it be a bounded 3D domain with appropriate surface boundary conditions to render the problem well defined. These boundary conditions may be specified as surface tractions and/or displacements, which are usually specified on all surfaces and possibly on interior surfaces such as faults. Free surfaces have zero surface traction by definition. Faults are interior surfaces, and may have associated dislocation increments at set times, when a stress threshhold is exceeded, or according to some other fault friction/triggering model. The solid domain may contain layers or other distributions of material with associated rhelogical properites. Currently supported materials are isotropic, newtonian elastic, newtonian viscoelastic, and nonnewtonian power-law viscosity.

Elastostatic solutions are supported, such as computing the displacements and stresses immediately caused by a specified slip distribution on a fault or finding the interior displacement and stress distribution due to a surface traction or displacement. These solutions are not time-dependent.

Viscoelastic solutions are also supported, in which the material flows and relaxes in response to imposed stress, such as an earthquake event. One may compute the viscoelastic reponse to a single event, or to multiple events in a sequence. The sequence may be user- specified, or may entail fault slips that are dynamically determined to occur according to a fault triggering model. Location-specific body forces are supported.

Boundary conditions and solutions apply to a finite-element discretized approximation to this domain. The domain is defined internally as a mesh of space-filling tetrahedral or hexahedral elements, with three components of displacement at each mesh node constituting the solution. Stress is computed for each element, and is element-wise constant for the current linear tetrahedral element type. Surface nodes carry special boundary conditions such as tractions or specified displacements. Nodes on faults are special split-nodes, that define screw or tensile dislocation on the fault without perturbing the mesh geometry. Temporal evolution is by discrete time steps using an implicit solution technique, allowing large time steps without numerical instability.

The code may be used with all nodes, elements, faults, boundary conditions and time history control created or modified by word processor, constituting the only needed input file. For large meshes hand-construction becomes impractical, and we have sought to support tools for automated mesh generation. Initial attempts have focused on having the user specify fault rectangles, materials, and mesh density, and then semiautomated tools produce the ascii input file.

Other supported features include: Specification of temporal epochs, each with differing steady boundary conditions Boundary velocity condition (steady change in the displacement, imposed) Controls for generating output on a subset of nodes/elements Control of implicit integration parameter Ability to shortcut temporal advance by fixing the sparse system for several time steps Control of checkpointing, saving state and allowing restart Checks a control file at each iteration, allowing clean interruption by the user.

Theory of operation

It is worth noting that a central theme to finite elements is the set of shape functions for individual elements. For each linear tetrahedron these are quite simple: for each node, define the function that is unity at that node and linearly declines to zero at the entire opposite face. We may also speak of global shape functions: for each element that shares a particular node, define the function that is unity at that node and declines linearly to zero at the opposite faces of each of the sharing tetrahedra. Shape functions so defined are called interpolatory: if we determine a list of displacements at each node, we may interpolate those values to any point interior to any element. Simply multiply each node displacement by that node's shape function, and sum over those product terms that are nonzero (these will be terms from the nodes defining the enclosing tetrahedron of the interpolation point). Such interpolated displacements will be continuous across all element boundaries (tetrahedral facets).

Within the volume of an isolated finite element, appropriate derivatives of the element shape functions weight the spatial derivatives of the stress tensor (defined from the strain and generalized Hooke's law, relying on locally defined element rheology) to enforce local equilibrium. The shape functions constitute a finite degrees-of-freedom approximation to the continuous system, and a volume integral enforces the equilibrium in a weighted average sense. Elements sharing a node contribute such weighted average terms to an equation for a single displacement. Away from boundary conditions and absent body forces, this set of displacement equation terms is set to zero. Body forces and surface tractions add forcing terms to the right-hand side.

The ensemble system of equations so defined is sparse, and if boundary conditions are sufficient the system is closed and solvable by standard sparse matrix techniques (currently by direct Crout factorization, but also possible to do by iterative techniques). Closed boundary conditions are usually easy to obtain; the user does have to ensure that enough of the boundary is constrained that there are no unconstrained body translations or rotations permitted to the domain.

The solution is the elastostatic solution for the posed problem in terms of displacements at each node. Local linear combinations of these displacements with shape function derivatives yield the stress tensor in each element.

The elastostatic solution is required for any viscous relaxation computation. Once the static step is complete, the time evolution of quasi-static viscous relaxation may begin by computing the viscoplastic strain rate, which is directly determined by the stress and the viscosity parameters. Conceptually this rate adds a force term to the right-hand-side of a sparse system similar to the elastostatic equilibrium. In practice, to obtain the advantages of implicit temporal development, terms are rearranged to modify the sparse equation coefficients as well. Each time step involves a solution to a sparse equation system, of similar cost as the elastostatic solution.

Faults are specified as either fault elements or as split nodes. With fault elements the user specifies fault properties such as failure criterion in an input record similar to that for a finite element. With split nodes the user specifies for each node on the fault its direction and amount of slip. The fault will slip at the initial time step, and also add the same increment of slip at each fault failure event, specified by the user. The split node formalism may be considered to represent the screw dislocation at a node as a separate entity from the displacement at that node, but is implemented here as an equivalent increment in the stress affecting the nodes immediately adjacent.


The user (perhaps via higher-level tools) is required to specify the model domain, constraints, temporal effects, and other solution controls in a single input file. Note that the model domain is a discretized version of the continuous domain described above. Therefore the volume is defined by its tetrahedral/hexahedral elements (or quad/triangular in 2-D), faults and other surfaces are sets of tetrahedral facets (triangles), and many constraints are specified on nodes. Higher-level constructs such as lines or surfaces may not have explicit descriptive records (although there is a construct for element groups sharing material properties; read on). This input file may be roughly characterized by the following brief list of items, which is in file order:
An output filename
Arbitrary comments
Number of processors and node decomposition (unsupported)
Number of free displacement dimensions
Number of time periods (0 for pure elastic; else number of steady epochs)
Flags: "Save shape functions", "Solver option"
Node free/fixed axis flags
Node coordinates
Node imposed conditions: displacement, velocity
Properties for clusters of elements: type, integration rule, material
Element records (nodes, property ID) (may be a fault element)
Split nodes
Surface tractions
List of nodes and elements for results printing
Number of time groups, plus reform interval, save time, quake time
Time group start times, implicitness parameter, time step
Output times (list, or flag for every earthquake)
Restart file name (optional) (i.e., start where another run left off, using this file)
Save-state file name (optional) (i.e., save this run so another run can start where we left off).

Outputs come in three flavors:
Standard Out: comments, summary info, progress messages.
Named output file (the name is in the first line of the input): the requested displacement and stress information at the requested times.
Restart file: contains simulation state. This is in ascii, but is intended for machine reading.

The input file affords considerable power to the user, but is very terse and cryptic. Without a full description of its options and features mistakes and misunderstandings inevitably arise. As a guide to the proper construction and formatting of a GeoFEST input file, we provide below a link to a hypertext-annotated template. The example is for a two dimensional problem, in order to facilitate easy visualization of the domain, however the generalization to 3-D is straightforward. Note the convention used in the input file that many vertical lists use the special flag "0 0" to indicate termination.

Example Format of 2-d Input File For GeoFEST Finite Element Program

Example Format of 3-d Input File For GeoFEST Finite Element Program

Reference Documentation

The following links lead to more detailed references describing the physical models and code structure underlying GeoFEST. These documents will assist investigators in understanding and modifying the source code.

Physics Technical Design Document
GeoFEST Source Structure and Subroutines

Contact Information

Please direct inquiries, requests, suggestions and bug reports to: