3d~def User's Manual
Introduction
Four files contain source code (Fortran text):
3D_MAIN.F, OKADA_SUB.F, XYZ_OUTPUT.F, VECTOR_OUTPUT.F
In order to compile the source code two include files are required. These are:
SIZES.INC, UNITS.INC
(Include files contain Fortran statements that are inserted into the
source code automatically upon compilation. They simplify making
changes to the program; for example, to change the maximum number of
allowable elements (see Section II.B) one would only modify a single
parameter specified in the file SIZES.INC. All appropriate array
dimensions within the four source code files would then be changed
appropriately without having to edit them.)
SIZES.INC contains program dimensions that depend on the number of
elements. The largest array (computer memory accessed as an indexed
matrix) has dimensions of (MAX3_ELEM)*(MAX3_ELEM+1) where
MAX3_ELEM=3*MAX_ELEM and MAX_ELEM is the maximum number of elements
allowed. You should make this number as large as you will need for your
specific problem by editing the file SIZES.INC before compiling the
program. If the program does not run and the computer returns an
'insufficient memory' error, you will have to reduce this number. The
memory useage could be made much more efficient if dynamic memory
allocation was implemented - this wasn't available for Fortran at the
time the program was written, and we haven't gotten around to updating
it.
UNITS.INC contains all file unit numbers ('tags' used by the computer
to identify a file when reading or writing to it). You can change these
as needed and re-compile the source codes to implement the changes.
The basic parameter specification (input) is done in the code
contained in the file 3D_MAIN.F. All major calculations are done in
routines contained in files 3D_MAIN.F and OKADA_SUB.F. Routines that
control the output grid and information written to column formatted
files are contained in the file XYZ_OUTPUT.F. Each row of these output
files contains the X,Y,Z coordinates of each grid point and specified
calculated quantities. These column formatted output files are read by
routines in the file OUTPUT_VECTOR.F and the information they contain
is reformatted in a vector format (as required by many plotting
programs such as MATLAB, IDL, etc.). Vector reformatting is optional.
The program does not generate any graphical output - displaying the output is left up to you.
II. Definitions and
conventions
A. Green's functions.
3D-DEF uses Green's functions calculated using subroutines provided by
Y. Okada. The Green's functions relate the deformation field
(displacements and displacement gradients) to a rectangular dislocation
in a homogeneous half-space (a semi-infinite medium bounded by a planar
free-surface). Descriptions of the dislocation theory underlying the
routines may be found in Okada, Y., Internal deformation due to shear
and tensile faults in a half-space, Bull. Seismo. Soc. Amer., v.
82, 1018-1040, 1992.
B. Selected definitions.
The meaning of an element and the coordinate systems used in the
computational algorithms are defined in this section. Three coordinate
systems are used. A local coordinate system is defined for each
element because the Green's functions (see Section II.A) are most
efficiently calculated in the local coordinate system of the element.
An in-plane coordinate system is also defined for each element
because the boundary conditions are specified either on or orthogonal
to the element surface. A global coordinate system is required
so that the effect of the deformation field associated with one element
on another element can be calculated. Because the parameters describing
the deformation are calculated in local coordinates, they must be
transformed into a global system that is common to all elements.
Figure 1
C. Boundary condition codes.
The boundary condition codes define the constraints on each
sub-element, one boundary condition for each directional (strike, dip,
normal) component of stress or displacement. Boundary conditions may
include the magnitude of the directional components of stress, of
absolute displacements, or of relative displacements on a sub-element.
For a single sub-element, mixed (stress and displacement) boundary
conditions may be specified.
The stresses are continuous everywhere and displacements are
discontinuous across elements. Understanding the difference between absolute
and relative displacements is very important when
specifying boundary conditions and interpreting the program output.
Once you understand the difference, you are well on your way to mastery
of boundary element modeling! The program calculates (solves for) relative
displacement on all of the sub-elements in the model such that
the strain energy within the half-space is a global minimum. Relative
displacement is simply the net dislocation or slip of one side
of the element with respect to the other side (for
example, the fault offset observed after an earthquake). The relative
displacements solved for by the program (that minimize the strain
energy) are determined by the user specified boundary conditions on the
sub-elements. The relative displacement on a particular sub-element or
group of elements may also be specified by the user, rather than being
treated as an unknown. Absolute displacement refers to
the motion of the footwall side of the element with respect to the
global coordinates (HOWEVER: the displacements are specified in the
planar coordinates of the sub-element; that is: strike direction, dip
direction, fault-normal direction). Mathematically, the relative
displacement is the difference between the absolute
displacements of either side of the element. For example, the normal
component of relative displacement, Dn, is
Dn = un- - un+
where un- is the magnitude of the motion of the footwall (absolute displacement) in the direction normal to the element and is un+ is the absolute displacement of the hanging wall.
The user specifies one boundary condition code (a number) for each
sub-element. Each code actually constrains the magnitude of three
parameters (stresses or displacements) at the center of the
sub-element, one in the strike direction, one in the dip direction, and
one in the direction normal to the sub-element surface. The boundary
condition codes allowed are:
Code#....Direction in the Plane of the
Sub-element...............................................
...........strike (shear)....................dip
(shear).......................normal (tensile)
1 ........ absolute displacement......absolute
displacement......absolute displacement
2 ........
stress...............................stress................................stress
3 ........ stress...............................absolute
displacement......stress
4 ........ absolute
displacement.....stress................................stress
5 ........
stress...............................stress...............................absolute
displacement
10 ...... relative displacement.......relative
displacement.......relative displacement
11 ...... angle (degrees)................stress (at specified
angle).relative displacement
12 ......
stress...............................stress...............................relative
displacement
13 ...... stress...............................relative
displacement.......relative displacement
14 ...... relative
displacement.......stress...............................relative
displacement
Note that Code #11 and Code #12 specify the same stress conditions; for Code #11 the user specifies the direction of shear as an angle measured from the strike direction, which is then resolved into shear stress conditions in the strike and dip directions.
D. Background deformation specification.
It is possible to include a uniform deformation field, referred to as
a background deformation field. (For example, suppose one wanted to
investigate how slip was partitioned among a network of faults driven
by a uniform simple-shear strain field.) The field is specified as a
stress, strain, or displacement gradient tensor. Internally these
tensors are converted to a stress tensor and a rigid body rotation (if
a displacement gradient tensor is specified). The background stresses
are subtracted from the stress boundary conditions before solving for
the relative displacements on the sub-elements. Because the background
deformation is added to that associated with the sub-element
displacements, the resultant stresses on the sub-elements are thus
guaranteed to equal the boundary condition stress value. Use of a
background stress or strain tensor will have no affect on the
sub-element relative displacements if you use Codes 1 or 10 because the
boundary conditions do not involve stress. If a displacement
gradient tensor is specified, then the only affect will be to modify
the output displacements by any rigid body rotations associated with
the background displacement gradient field.
Section III.B explains how to specify the
input to include (or not include) a uniform background deformation
field.
F. Units.
Specify input element (fault) dimensions in kilometers and
displacements in centimeters. Output displacements have units of
centimeters, strains are dimensionless. Stresses (input and output)
have the same units as the Young's modulus specified in the input file.
1 Sxx
2 Sxy
3 Sxz
4 Syy
5 Syz
6 Szz
A. Introduction to the input file.
The program requires one input file and the amount of interactive
input is variable depending on what options the user selects. All input
parameters may be written to a file so that the only interactive step
is to respond to the program prompt for the input file name.
Alternatively, the user can interactively enter the type and
characteristics of the information written to output files.
All lines in the input file containing numbers are free formatted so
the format is not important but items cannot be omitted. Text lines
that begin with an asterisk are provided only for the user's benefit
and are ignored by the program. Note that while you can write anything
in these lines, the number of asterisk lines and blank lines must not
change.
The first part of the input file contains parameters that describe the
material properties of the medium, the element geometry, and the
boundary conditions. The second part of the file controls which
properties of the deformation field are calculated and output and where
within the half space (and in what coordinate system) these properties
are calculated. This second part of the file may be omitted and the
needed input entered interactively.
The next group of numbers (following a blank and three comment lines)
specify the boundary conditions for each sub-element, one row per
sub-element. Figure 2 below illustrates
how the sub-elements are identified. The first two integers are
sub-element indices but are not actually used by the program; they
simply help the user keep track of the sub-elements. The first
sub-element must be the one at the internal reference point or
equivalently, the origin of the in-plane coordinate system (looking at
the plane from the hanging wall the first sub-element is in the lower
(deepest) left corner). The element dip index is the second number of
the pair (right index) and the values increase up-dip. The first number
of the pair is the strike index and these increase with distance from
the origin. For each strike index, the dip index is incremented .
Diagramatically the indices of the sub-elements of the first plane in
the sample file are:
Figure 2
If output of parameters describing the deformation on a plane is
desired, the line must start with the character string
"plane". For a gridded volume, the line must start with
"vol". (For a gridded volume, all output parameters, except
element relative displacements, will be in global coordinates, no
matter what you write in the appropriate line). For output at a set of
user-specified points, the line starts with "user". In
almost all cases, only the first letter of the line is read by the
program.
The next line determines the type of output, "XYZ" or
"vector". For XYZ (column) output, the output files of
deformation parameters contain one line per point, with the X, Y, and Z
coordinates of the point given one on each line. For vector output, the
values within the gridded region (which may be either a plane or a
volume) are output into a file as an array of numbers, without
coordinates. The coordinates of each point are determined by their
position within the array (Seem Obscure? Patience, Grasshopper!
Enlightenment may be obtained by reading more about this issue below in
the following examples!). NOTE: if vector output is requested, both
column output files and vector output files are created because
the column output files are always created, and the program generates
the vector files from the column output files.
The next block (#2) determines the coordinate system for calculating various parameters describing the deformation field. The two choices of coordinate systems are "global" or "in-plane"; the latter is specified with anything that doesn't start with a 'g' or 'G' (in-plane, planar, etc.). This option only applies to output on a gridded plane and affects only output stresses, strains, displacements, and gradients of displacements. Global output coordinates means that stresses, strains, displacements, and gradients of displacement are coordinate stresses, strains, etc. in the same reference frame in which the faults were specified (that is, Sxx, Sxy, Sxz, etc.). In-plane (or planar) coordinates refer to the gridded plane, that is: the output stresses, strains will be referred to the strike, dip, and normal direction of the gridded plane (that is, Sss, Ssd, Ssn, etc).
The next block (#3) determines which parameters describing the deformation field to calculate. The choices are: element displacements (that is, one side relative to the other of each element), and at each grid or user specified point, absolute displacements, gradients of displacements, strains, stresses, strain invariants, strain orientations, rigid body rotations, and failure plane orientations. These are described in detail in Section IV. The list of files containing various characteristics of the calculated deformation field can be in any order. The word on each line sets an appropriate flag so that an output file containing the parameters described by the word is generated. The program only interprets the first letter (the first 4 letters for 'stress' and 'strain') so there is some flexibility in what words are used. For example, you must write 'gradients of displacement' (if fact, simply 'g' will do) to generate an output file containing the components of the displacement gradient tensor; use of 'displacement gradients' would be incorrect and (because the key word starts with a 'd' instead of a 'g') result in an output file containing displacements.
The next block (#4) is a line containing the suffix that will be appended to the output files.
The final block (#5) specifies the points at which output is calculated. If a user array was specified in block #1, then this portion contains the number of points, and their coordinates (see example below). If a gridded plane was specified in block #1, then this section contains the orientation and size of the plane in one line, and the number of points on the plane at which to calculate the output on a second line (see example below). If a volume is specified, the first line of this block contains the coordinates of two of the corners of the rectangular prism of the volume, and the second line contains the number of X, Y, and Z grid points within the volume.
In the examples below, user input that effects the outcome of the model
(that is, which output files are created, what they are called, where
the values are calculated, etc.) is shown in bold. The minimum
input required to obtain the desired result is underlined. Everything
else in the second part of the file is ignored by the program, but is
included to make the input file more intelligible to the user.
(NOTE: lines that begin with '*' must be in the input file
because they are used by the program to separate blocks of differing
types of input. However, they may contain any text following the '*'. )
A. Column output format.
Nine output files are created with column format. Each one has a
similar format. The suffix ??? in the following discussion is replaced
with the suffix that you specified in the input file (or answered
interactively). Each line contains an observation point x,y,z (in
global coordinates) and the parameters describing the various
characteristics of the deformation field at that point. The file names
and their contents are:
The first eight output files may be re-formatted in a vector format (that is, as desired by plotting programs such as MATLAB, IDL, etc.); this re-formatting is of course optional. Remember that if you change the format of these column formatted files, you must also change the routines in OUTPUT_VECTOR.f since they must read these files. The element displacements may also be written out in vector format but this is independent of the column formatted file.
where # is the number of the plane as listed in the input file. Each contains a component of the displacement such that each line of the file corresponds to the sub-elements along strike at a fixed depth (dip index). The first line for each plane is the top (closest to the surface) of the plane and lines increment with increasing depth. A set of lines is written for each plane, and each plane is separated by a blank row.
V. Acknowledgements
The authors would like to thank P. Bodin for his input and discussions
on 3D-DEF and boundary elements in general, bravery in testing the
program, and assistance in writing this manual. We also thank C.
Meertens and B. Anderson for being guinea pigs, and UNAVCO/UCAR and H.
Benz for providing computing facilities. A. Crone and K. Haller
provided extremely useful reviews of this manual.