3D-DEF: A three-dimensional, boundary element program by Joan Gomberg (USGS) 
with Michael Ellis (CERI, Memphis State University).   January 8,1993.


If you have figured out elastic deformation in 2-dimensions and
are game for a new challenge, here is the first iteration of our 3-d 
boundary element (BE) program.  Unlike the 2-d BE algorithm of Crouch and
Starfield implemented by G. King this is not tied to any computer (as long
as whatever computer is used has a Fortran compiler) or graphics system.  
This mini-manual (which will grow to a proper manual as we receive friendly 
suggestions) contains the following information:
	1. Getting the program to run.
	2. Some random comments about the program.
	3. Input file description.
	4. Output file descriptions.
The code does have some comments inside of it for those brave enough to 
look into it.  Please note that as of today, the program has not been
run extensively (an understatement) although many basic tests have been
performed.  We will attempt to notify you herein when you are entering
unchartered waters.  

We are making this code freely available with the stipulation that users will
acknowledge the following reference when necessary:

Seismicity and surface deformation in the New Madrid region: Numerical
experiments using a new three-dimensional boundary element program.  
(J. Gomberg and M. Ellis, In preparation for BSSA/JGR/?)

We will update the reference as updated versions of the program become available. 

1. Getting the program to run.
   There are three files containing source code; 3D_MAIN.F,OKADA_SUB.F,XYZ_OUTPUT.F,
   VECTOR_OUTPUT.F These require two include files; SIZES.INC, UNITS.INC.  
   SIZES.INC contains program dimensions that depend on the number of elements.  
   The largest
   array has dimensions MAX3_ELEM by (MAX3_ELEM+1) where MAX3_ELEM=3*MAX_ELEM
   and MAX_ELEM is the maximum number of elements allowed.  UNITS.INC 
   contains all file unit numbers.  You can change these to suit your
   needs and recompile the three source codes to implement the changes.
   The codes are written in generic fortran.  They were developed using
   Macintosh Language Systems Fortran but should compile with most any
   Fortran compiler (no structures are used, and even do-enddo's are not
   used for those with old systems).
   
   The basic input is done in 3D_MAIN.F and all the serious calculations are
   done in 3D_MAIN.F and OKADA_SUB.F.  All routines that control the output
   grid and choice of information written to X,Y,Z,Data formatted files are 
   contained in XYZ_OUTPUT.F. These X,Y,Z,Data formatted output files are read 
   by routines in the file OUTPUT_VECTOR.F and re-formatted in a vector format 
   (e.g., as required by plotting programs such as MATLAB); this re-formatting
   is of course optional.

   	  Displaying the output is left up to you. See the section on output
   files for further details.
   
   
2. Some random comments about the program.

  A. Green's functions are based on the dislocation code of Y. Okada as
    described in BSSA, Vol. 82, p. 1018, 1992. We are currently looking into
    adding the capability of using triangular elements.
	
  B. Some definitions.
    Plane = Element = Fault or fault segment with constant strike and dip.
 	Sub-element = each plane can be subdivided into smaller sub-elements.
                  For each plane these have the same dimensions equal to
                  the total length/number of sub-elements along strike and
                  the total width/number of sub-elements along the dip direction.
  	Global coordinates = Coordinate system for all planes.  X is E, Y is N, and
  	          +Z is vertical pointing up wrt the free surface.
 	Reference point,Xo,Yo,Zo = For each plane determines the relative position
 		  of all planes in the global coordinate system. For each plane
 		  the user specifies it as the corner of the plane closest
 		  to the surface such that looking away from the reference
 		  point the fault strikes with the hanging wall on the right.
 		  Internally Xo,Yo,Zo is redefined as the same end of the
 		  fault plane but at the lower corner.
    Local coordinates = Coordinate system for each plane as defined by Okada;
                  local XL axis is along strike, YL axis perpendicular to strike,
       		  ZL axis is vertical. The system is right handed such that
 		  looking in the strike direction (+XL) the hanging-wall is on the
 		  right and the +YL axis points into the footwall.  The reference 
                  point,Xo,Yo,Zo (specified in global coordinates) as specified 
                  by the user is the corner of the plane closest to the surface 
                  such that looking away from the reference point the fault 
                  strikes in the +XL direction with the hanging wall on the right.
  	Planar (in-plane) coordinates = Coordinate system for each plane such that XP 
 		  is along the strike direction (=XL), YP is along the dip direction
 		  and is positive going up-dip from the internal Xo,Yo,Zo. +ZP is
 		  normal to the plane with positive defined so that the system
 		  is right handed.
 
  C. Boundary Condition Codes (see the section on input file).
    Boundary conditions currently allowed:     
 	KODE	strike (shear)	dip (shear)		normal(tensile)
 	 1	displ.		displ.			displ.
 	 2	stress		stress			stress
  	 3	stress		displ.			stress
 	 4	displ.		stress			stress
 	 5      stress		stress			disp.
 	10	relative	relative		relative 
 		displ.		displ.			displ.
 	11	angle 		stress (along	        relative 
 		                angle)			displ.
 	12      stress		stress			relative
 	                                                displ.
 	13	stress		relative		relative 
 	                        displ.			displ.
 	14      relative	stress			relative 
 			        displ.			displ.
   Note that 11 and 12 are the same; for KODE=11 the user specifies
   the direction of shear as an angle measured from the
   strike direction which is then resolved into shear stress conditions
   along the strike and dip directions.
 
  D. Directions 
    Fault (element) orientations are specified as follows:
    1. When looking in the strike direction the hanging wall is on the right.
    2. The strike angle is defined clockwise from north.
    3. Dip is defined from the horizontal, from 0 to 90 degrees.
 
    Fault (element) relative displacement conventions are as follows:
    Dislocation Component	Sign		Faulting Mode
 	   Strike	        >0			left-lateral
 	    			<0			right-lateral
 	   Dip    		>0			thrust
 	  	 		<0			normal
 	Tensile 
	(normal to plane)	>0			opening
 		                <0			closing
 
    Displacement boundary conditions refer to displacement of the footwall side
    of the element.
 
  E. Units 
    Specify input element (fault) dimensions in km, displacements in cm.
    Output displacements with have units of cm, strains will be dimensionless.
    Stresses (input and output) have units of whatever the Young's modulus
   is given in.
 
  F. Misc.
    Stresses are written as a vector internally as
    Vector Index	Stress Component
   		1	Sxx
 		2	Sxy
 		3	Sxz
 		4	Syy
 		5	Syz
 		6	Szz


3. Input file description.  
  Presently one input file is required and the amount of interactive
 input is varied. All input can be done without human interaction other
 than responding to the program prompt for the input file name.  
 Alternatively, the user can interactively enter the type and character-
 istics of the information written to output files.  
 	All lines containing numbers are free formatted so that the format
is not important but items cannot be omitted.  Text lines are meant
only for the User's benefit and begin with an asterisk (ignored by
the program).  Note that while you can write anything in these lines,
the number of lines cannot change. 
	The first part of all files is the same, containing parameters
that describe the element geometry and boundary conditions.  The 
second part of the file controls the output and differs depending 
on whether the program is being run interactively, whether gridding
is being done on a plane (of any orientation) or a volume (with
vertical and horizontal sides), or output at a specified list of points is
used.  
	
	The unchanging, first part of the input file is described as follows:
The first line of numbers is basically self-explanatory; recall that
a plane is the largest element and can be sub-divided into sub-
elements.  The coefficient of internal friction is used to determine the 
potential failure planes, assuming a Coulomb failure criterion.  See 
Jaeger and Cook, Fundamentals of Rock Mechanics, 3rd ed., 1979,
section 4.6, p. 95.  The final item in this line is a character
string flag set to 'none' if no background deformation field is 
specified, to 'stress' if a background stress field 
is specified, to 'strain' if a background strain
field is specified, or to 'displacement
gradient' if a background displacement gradient field is specified.

The next line is blank if the above flag is 'none' or contains
the elements of the background deformation field.
The order of stress or strain tensor components is xx,xy,xz,yy,yz,zz.
For a flag set to 'displacement' the order of displacement gradient 
tensor components is dUx/dx,dUx/dy,dUx/dz,dUy/dx,dUy/dy,dUy/dz,
dUz/dx,dUz/dy,dUz/dz.

  The next group of numbers describes each plane, one line
per plane.  The first three numbers are the origin of the plane in
global coordinates (see Section 2B), fourth and fifth numbers are the
length (in the strike direction) and width (dip direction) of the plane,
the sixth and seventh numbers describe how the plane is sub-divided
into sub-elements in the strike and dip directions.  Sub-elements
on a single plane all have the same dimensions; lengths equal the
plane length divided by the sixth number (under STK in the sample file)
and widths equal the plane width divided by the seventh number
(under DIP).  The last two numbers describe the strike and dip of
the plane (see Section 2D).  In this sample file there are two planes 
specified, the first is divided into six sub-elements and the second 
into two sub-elements. 
  The next group of numbers specify the boundary conditions for each
sub-element, one row per sub-element.  The first two integers are sub-
element indices but are not actually used by the program - they are 
meant to help the  user keep track of the sub-elements.  The first
sub-element must be the one touching 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 (right) index and increases up-dip.  The first index is the
strike index and increases in increasing strike direction.  The
dip index is incremented for each strike index.  Diagramatically
the indices of the sub-elements of the first plane in the sample file
are

       +dip direction
                |
		|-----------------------------------------
		|                    |                    |
		|       1,3          |      2,3           |
		|                    |                    |
		 -----------------------------------------
		|                    |                    |
		|       1,2          |      2,2           |
		|                    |                    |
		 -----------------------------------------
		|                    |                    |
		|       1,1          |      2,1           |
		|                    |                    |
		 -------------------------------------------- +strike direction


	The changing, second part of the input file is now described. For
a non-interactive run with a gridded volume the second part
of the file might look like this (ignore dashed lines):
------------------------------------------------------------
* Run everything else non-interactively
Not Interactive  (alternative is 'Interactive' grid specification)
Volume grid (alternative is 'Plane grid')
XYZ Output  (alternative is Vector output - results in both output types)
* Coordinate system for output (global only for volume grids)
global
* Output information (each line results in a separate output file)
stress (in global coordinates)
strain (in global coordinates)
gradients of displacement (in global coordinates)
rigid body rotations (about the global X(E),Y(N),Z(up) axes)
orientation of principal strain axes (w.r.t the global X(E),Y(N),Z(up) coordinates)
displacements (in global coordinates)
invariants of the strain field (independent of coordinate system)
failure planes (in global coordinates)
element relative displacements (strike,dip,normal directions of each element's plane)
* Output file suffix
out1
* XYZ Coordinates of the upper NW and lower SE corners of the volume:
0,0,0,90,90,10
* Grid spacing increments in X (E), Y (N),Z (up) directions:
5,5,1
------------------------------------------------------------
	Note that 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). Remember that lines 
beginning with * are only descriptions and aren't used
by the program (they can say anything).  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 will be generated.  
The program only interprets the first letter (the first 4 letters for
'stress' and 'strain') so that their is some flexibility in what words
are used.  For example, you must write 'gradients of displacement' or
simply 'gradients' to generate an output file containing the components
of the displacement gradient tensor; use of 'displacement gradients'
would be incorrect (the key word must start with a 'g') and result in
an output file containing displacements.


	For a non-interactive run with a gridded plane the second part
of the file might look like this (ignore dashed lines):
------------------------------------------------------------
* Run everything else non-interactively
Not Interactive  (alternative is 'Interactive' grid specification)
Plane grid (alternative is 'Volume grid')
* Coordinate system for output (global only for volume grids)
local (local coords of gridded plane)
XYZ Output  (alternative is Vector output - results in both output types)
* Output information (each line results in a separate output file)
stress (coordinate system of the gridded plane;X=strike,Y=dip,Z=normal directions)
strain (coordinate system of the gridded plane;X=strike,Y=dip,Z=normal directions)
gradients of displacement (coord. system as for stress and strain)
rigid body rotations (about the global X(E),Y(N),Z(up) axes)
failure planes (about the global X(E),Y(N),Z(up) axes)
orientation of principal strain axes (w.r.t the global X(E),Y(N),Z(up) coordinates)
displacements (coord. system of the gridded plane;X=strike,Y=dip,Z=normal directions)
invariants of the strain field (independent of coordinate system)
element relative displacements (strike,dip,normal directions of each element's plane)
* Output file suffix
out2
* Xo,Yo,Zo (Ref. corner),strike (cw wrt N),dip (0-90),length (strike), width (dip):
0,0,0,90,90,10,5 
* Number of grid points in the strike and dip directions:
5,5
------------------------------------------------------------
	Note that for a gridded plane, if the coordinate system requested
is anything but 'global'  output parameters that are referenced to the 
global coordinate system include: 
		rigid body rotations
 		orientation of principal strain axes
		failure plane descriptions
output parameters that are in the coordinate system of the gridded plane 
(X=strike,Y=dip,Z=normal directions) include:
   		stress tensor elements
   		strain tensor elements
   		displacement gradients tensor elements
 		displacements
Element relative displacements correspond to the strike,dip,normal directions 
of each element's plane. As for the volume, remember that lines 
beginning with * are only descriptions and aren't used
by the program (they can say anything). The list of files containing various
characteristics of the calculated deformation field can be in any order but
the first character must be correct (see gridded volume description above).

      For a non-interactive run on a user-specified array of points
------------------------------------------------------------
* Run everything else non-interactively
Not Interactive  (alternative is 'Interactive' grid specification)
User array (alternative is 'Plane grid' or 'Volume grid')
XYZ Output  (alternative is Vector output - results in both output types)
* Coordinate system for output (global only for volume grids)
global
* Output information (each line results in a separate output file)
stress (in global coordinates)
strain (in global coordinates)
gradients of displacement (in global coordinates)
rigid body rotations (about the global X(E),Y(N),Z(up) axes)
orientation of principal strain axes (w.r.t the global X(E),Y(N),Z(up) coordinates)
displacements (in global coordinates)
invariants of the strain field (independent of coordinate system)
failure planes (in global coordinates)
element relative displacements (strike,dip,normal directions of each element's plane)
* Output file suffix
out1
* Number of points followed by x,y,z
5
10 5 -5
12 6 -5
14 5 -5
16 6 -5

	 For an interactive run the second part
of the file might look like this (ignore dashed lines):
------------------------------------------------------------
Run everything else interactively
Interactive grid specification
___________________________________________________________________
The User will be asked by the program to answer lots of questions.
 	
	

4. Output file descriptions.

   Presently nine output files may be created with X,Y,Z,Data format. 
   Each one has a similar format.  The suffix (??? below) is whatever you specified
   in the input file (or answered to the prompt).
   - each line contains an observation point x,y,z (in global coordinates)
   and parameters describing various characteristics of the deformation field
   at that point.  The characteristics in each of these files are as follows:
     File 		         Contents
     stressten.??? - Six independent elements of the stress tensor, 
 	          Sxx,Sxy,Sxz,Syy,Syz,Szz, and the magnitude and orientation
		  of the maximum shear stress on the inspectoin plane, if specified.
		  If output on a plane is 
                  requested these correspond to the in-plane coordinates 
		  (e.g., Sss,Ssd,Ssn,Sdd,Sdn,Snn where 's' refers to
		  the strike direction, 'd' to the dip direction, and
		  'n' to the normal direction.  Thus Szz is really
		  Snn and corresponds to the stress normal (tensile)
		  to the fault plane, Sxz corresponds to Ssn and is
		  the traction in the strike direction.)
     strainten.??? - Six independent elements of the strain tensor, 
 		  Exx,Exy,Exz,Eyy,Eyz,Ezz, and the volumetric strain 
 		  (Exx+Eyy+Ezz,independent of the coordinate system). If  
 		  output on a plane is requested these correspond to  
 		  the in-plane coordinates as in stressten.out.
     displacements.??? - Displacements in the x,y,z directions. If output 
 		  on a plane is requested these correspond to the 
 		  in-plane coordinates; Ux is in the strike direction
 		  Uy in the up-dip direction, Uz normal to the plane.
     dgten.??? -  Displacement gradient tensor; dUx/dx,dUx/dy,dUx/dz,
 		  dUy/dx,dUy/dy,dUy/dz, dUz/dx,dUz/dy,dUz/dz.  As for
 		  the above, these correspond to the in-plane coordinates
 		  if the user has selected the plane option.
     strainorient.??? - Principal strains,their plunge and trend with respect
 		  to the global axes.  They are listed from maximum to
 		  minimum.
     invariants.??? - Strain & stress invariants: the volume change, failure 
 		  stresses,octahedral shear stress,strain energy density.
     rbr.- Rotations about the global x,y,and z axes
                  Output is in degrees; positive rotations are counter-clockwise                  negative are clockwise.
     failure.??? -  For each of the two failure planes:strike and
	          dip of the plane and the direction cosines of
		  the slip vector.  The strike is cw from north,
		  the dip is from the vertical, and the direction
		  cosines are wrt X(E),Y(N),Z(up).
     elements.??? - Each line of the output file will contain the
		  coordinates of the center of the element in global
		  coordinates, the center in in-plane coordinates,
		  and the relative displacements in the strike,
		  dip, and tensile (normal) directions.
	
	The first eight files are opened in open_opts, closed in close_opts,
and written to in XYZ_Results. The nineth file, elements.???, is opened and 
closed in subroutine final_write and written to in XYZ_elements.
 	The first eight output files may be re-formatted in a vector format (e.g.as 
required by plotting programs such as MATLAB); this re-formatting is of course optional.   
Remember that if you change the format of these X,Y,Z,Data 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 X,Y,Z,Data file.


	If the User requests Vector output deformation parameters in 
will  be re-written in vector format (in new files). Each parameter
will be written to its own file.  If the program is being
run non-interactively, all parameters of a particular selected
type will be written out.  For example, if the User has
specified 'gradients of displacement' in the input file, 
9 output file will be generated - one for each element of
the displacement gradient tensor. If running in an interactive 
mode, the User can select which parameters to be written out 
(e.g., only the dUx/dz, dUy/dz,dUz/dz components of the
displacement gradient tensor or whatever subset is desired).
See the description of the X,Y,Z,Data formatted files to 
find out what coordinate systems the output parameters are
calculated for. Note that for in-plane coordinates,
x is the strike direction, y is the dip direction, and z is
normal to the plane.  For a gridded plane each line of the 
output file corresponds to a row (constant depth) on the plane; 
going from left to right on each line corresponds to the increasing
strike direction.  For a gridded volume consider the volume to
be a stack of horizontal planes, the first plane being the one
closest to the surface.  The output is the same as for the
gridded plane for each of the horizontal planes.  The first
value of each plane corresponds to the northwest corner.
	The following describes the output file possibilities for each
output type requested.  
    Output type		Files		Contents
  	element	        elemds#.???      strike component of displacements 
   	 displacements	elemdd#.???      dip component of displacements     	
      		        elemdn#.???      tensile component of displacements
 				        #=index number, 1... no. of planes

 	stress		Sxx.???        xx component of the stress tensor
 		  	Sxy.???        xy component of the stress tensor
 			Sxz.???        xz component of the stress tensor
 			Syy.???        yy component of the stress tensor
 			Syz.???        yz component of the stress tensor
  			Szz.???        zz component of the stress tensor
			tmax	       max shear stress magnitude on the plane
			tmaxa          orientation (rake) of tmax
			tmaxo(2)       vector corresponding to tmaxa (for MATLAB users)
 	strain		Exx.???        xx component of the strain tensor
 		  	Exy.???        xy component of the strain tensor
 			Exz.???        xz component of the strain tensor
 			Eyy.???        yy component of the strain tensor
 			Eyz.???        yz component of the strain tensor
 			Ezz.???        zz component of the strain tensor
 	gradients in
 	 displacements	dUxdx.???       gradient of x-displacement in x-direction
 			dUxdy.???       gradient of x-displacement in y-direction
 			dUxdz.???       gradient of x-displacement in z-direction
 			dUydx.???       gradient of y-displacement in x-direction
 			dUydy.???       gradient of y-displacement in y-direction
 			dUydz.???       gradient of y-displacement in z-direction
 			dUzdx.???       gradient of z-displacement in x-direction
 			dUzdy.???       gradient of z-displacement in y-direction
 			dUzdz.???       gradient of z-displacement in z-direction
        rigid body 
 	 rotations      rbr_deg_x.???	rotation about the X (E) axis
 			rbr_deg_y.???	rotation about the Y (N) axis
 			rbr_deg_z.???	rotation about the Z (vertical) axis
 	orientation of
     principal axes	prince_max.???  magnitude of the max. principal strain
 			plunge_max.???  plunge wrt horiz. of the max. principal strain 
 			trend_max.???	trend (cw wrt N) of the max. principal strain
 			prince_int.???	magnitude of the intermed. principal strain
                        plunge_int.???	plunge wrt horiz. of the intermed. princ. strain
 			trend_int.???	trend (cw wrt N) of intermed. principal strain
 			prince_min.???	magnitude of the min. principal strain
 			plunge_min.???	plunge wrt horiz. of the min. principal strain
 		        trend_min.???	trend (cw wrt N) of the min. principal strain
     displacements	disp1.???	displacements in the x-direction.
 		        disp2.???	displacements in the y-direction.
 		        disp3.???	displacements in the z-direction.
 	invariants      delv.???	volumetric strain
 		        aftershocks.???	failure stress
 			oss.???		octahedral shear stress
 			work.???	strain energy density
        failure planes  strike1.???	strike of first failure plane           
                        dip1.???	dip of first failure plane           
                        ce1.???	        slip vector dir. cos. wrt X(E) 
                        cn1.???	        slip vector dir. cos. wrt Y(N) 
                        cz1.???	        slip vector dir. cos. wrt Z(up) 
                        strike2.???	strike of second failure plane          
                        dip2.???	dip of second failure plane           
                        ce2.???	        slip vector dir. cos. wrt X(E) 
                        cn2.???	        slip vector dir. cos. wrt Y(N) 
                        cz2.???	        slip vector dir. cos. wrt Z(up) 


	For the element displacement output on each plane 3 file types are 
created for each plane -
      elemds#.???	   strike component of displacements 
      elemdd#.???	   dip component of displacements 
      elemdn#.???	   tensile component of displacements 
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 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 are written for each plane, each plane is separated
by a blank row.


Supplemental Documentation

15-Feb-1994

Calculation of maximum shear stress on an inspection plane.

The magnitude and orientation of the max shear stress on a 
specified inspection plane is calculated in the subroutine
TOPLANE.  Tau max is obtained in the following way: the in-plane
stress tensor, S, contains two components that represent shear 
stress across the inspection plane: tau(zx) and tau(zy), where -
recall - x is along strike, y is up dip, and z is normal to the 
plane.  In order to find an expression for tau max, we first express
tau(zx) as a function of an arbitrary rotation about the z-axis.
Let the rotation matrix be r= [cos(a) sin(a) 0; -sin(a) cos(a) 0; 0 0 1],
where "a" is the angle of rotation.  Thus, from: S'=r.S.r' we obtain an
expression for tau(zx)'=tau(zx).cos(a)+tau(zy).sin(a).  The maximum occurs
when the dtau(zx)'/da = 0, i.e. when a=atan(tau(zy)/tau(zx)), and its 
orientation is given by "a".  

This is a useful set of results since the distribution of tau max can be 
compared to the slip direction.  In stress inversion techniques, these two
directions are assumed to be the same.

The results are written to stressten.suffix and to separate vector files,  
tmax.suffix (magnitude), tmx. and tmy. (vector components of the orientation).
The latter files are good for the QUIVER function inside MATLAB.

*******************************************************************

18-Feb-1994


Modifications to failure.sub
by 
Mike Ellis and Debi Kilb

This new version of 3d_def has the coding for failure.f to return
values strike, dip and rake at the user specified points.

Routines that were changed:
	failure.f    (Complete overhaul)
	xyz_output.f (calling sequence and what is printed)

Routines that were added:
	matmult

Files that were changed:
	3dmain.f
	xyz_output.f

Files that were changed to get rid of common block errors.
 
	3dmain.f
	okada_sub.f
	xyz_output.f
*******************************************************************
