Some Notes about Grid Generation

OVERVIEW

Philippe F. has written code to generate a grid. This code runs as part of the
interactive pre-processing part of the thermo-mechanical model.

The grid generation code allows the user to define the different sections of the
grid by defining the characteristics of the "interface" between each part of the
grid.  The code then uses these definitions to interpolate and create a grid.

Figure 1. Example of a desired grid.


(have one or two examples, covering all the types of interface)

INPUT

The grid generation code asks for the following parameters:

hx   width of plotting window
hz   height of plotting window

     hx, hz are the width and height of plotting window (in the same units
     as used for the grid - normally metres)
     (0,0) is the assumed origin at upper left
     NOTE: hx should be the same width that you want your grid to be.
     So if you want a grid 1000 km wide, use hx of 1000000.

igy  number of nodes in the horizontal direction, for the thermal grid
igz  number of nodes in the vertical direction, for the thermal grid

You are then asked for the number of lines.  This is the number of
interfaces in your grid.  In the example below, there are 4 of these
interfaces.

For each interface, you need to enter:

   * The row number in the grid for this interface
   * What type of interface:
       a) line-point
       b) translate previous interface
       c) create a constant thickness slab
       d) line-circle-line interface
       e) circular transition from a horizontal to a horizontal 
	 (e.g subduction) geometry
       f) read-in interface
       g) line segments
     (these types of interfaces are explained below)
   * Parameters needed to create the chosen type of interface (explained
     below)

TYPES OF INTERFACES

a) line-point

   This type of interface is a straight line.  You define it by giving the
   (x,z) coordinates of the first point on the line, followed by the
   (x,z) coordinates of the second point on the line, followed by the number
   of points on the line.
   
   The number of points used determines how well the resulting grid lines
   follow the desired interface.  More points, better fit... but you must
   keep in mind that there is a limit to the total number of points in
   all the interfaces. (Looks like 1000 points at the moment - BL)
   
   In general, you should  put the endpoints of the line slightly outside the
   desired grid area, so the code can extrapolate accurately at the edges of the
   grid.
   
   Example: 
   
 enter reference line npoints and xa za xb zb
0 0 605000 0 201

   This defines a line from (0,0) to (0,605000) with 201 points on the line.

   NOTE for this example:  The prompt is wrong!  The order of the input should
   be: xa za xb zb npoints
   where xa is the x coordinate of the first point of the segment
         za is the z coordinate of the first point of the segment
         xb is the x coordinate of the end point of the segment
         zb is the z coordinate of the end point of the segment
         npoints is the number of points on the line
         
   The "first" point is the one with the smaller x value.
   
       
       
       GENERAL NOTE: Philippe says, be careful that the total number of 
         points on all interface segments is not more than 1000.

     
b) translate previous interface

   This creates an interface which is the same shape as the previous one, but
   translated in the z direction.
   
   [ Add more info here]
   
c) create a constant thickness slab

d) line-circle-line interface

   (diagram)
   
   Use this interface type when one or both of the lines are not horizontal.  If
   both lines are horizontal, use "circular transition from a horizontal to a
   horizontal geometry".
   
   

e) circular transition from a horizontal to a horizontal (e.g. subduction)
   geometry
   
   (diagram)
   
   
   
   NOTE:
   For horizontal-circle-horizontal transition, can only specify horizontal 
   position of singularity, exit angle, and depth of slab.  Can't specify the
   delta x (horizontal distance spanned by the arc); have to let the program
   place the top of the circle.

   Example:
   
 enter horizontal singularity position
300000
 enter zleft zright
55588 20000   [use positive numbers; "z" here is not right; we're looking at
               depth]
 enter slab angle
20
 enter slab length
500000    [use a big number to make sure it goes past edge of grid]
 radius is  590111.1932229433 55588.0 20000.0
 confirm :delta_z=  35588.0
 top circle position  501829.7513867472 -20000.0
 sla1 point  300000.0 -55588.0
 enter angle left
0
 500000.0 300000.0 3.14159 353553.625139613 -200000.0
 -200000.0 -55588.0 300000.0 -55588.0
 enter points on exit-slab ,circle,in-slab
101 31 101
 enter position of top of in-slab
605000 -20000

   
f) read-in interface

   [check the code about this... or ask Philippe]



==========================================================================

INPUT file notes:

in microfem_tuning_i, set ifpa = 1 if we want to generate result_track
(ifpa = i first pass)

[cobra]/disk4/bonny/ptt/c13q_v1_e1-1.75my_rb1 % ../CODE960916/xnew
 1=1/3 formul 2 = 1/2 formul
1
  1e+38=  1.000000000000000E+38
 uncoupled run
 TIME INITIALIZATION :  .0 473041372198.4124
 igridn  0
 201 28 5628 301 28 8428 301 55
 generate (1=interactive 2=file) or read (0) mesh
1
 make thermal grid .choose here your hx hz
 exemple 400000. 100000.
 
        [hx, hz are width and height of plotting window (in the same units
         as used for the grid - normally metres)
         (0,0) is the assumed origin at upper left]
         
600000 300000
 your z is exagerated (or equivalently your x
 is shrunk by a factor  4.0
 enter igy igz

        [igy is number of nodes in horizontal direction in the thermal grid.
         igz is number of nodes in vertical direction in the thermal grid. ]
201 55
io_open_xwindow: Opening X11 window
    [ displays window with title Plot "E mesh 2" with a grid]
 enter number of lines
     [ number of lines = number of depth interfaces ]
     [  "line" (ligne) = "row" ]
4
 line number  1
 enter line number in grid
1
 new menu !!!
 enter 0 for line-point
 enter 1 to translate previous interface
 enter 2 to create a constant thickness slab
 enter 3 for line-circle-line interface
 enter 5 for circular transition from
 a horizontal to a horizontal (e.g subduction)
 geometry
 enter 4 for read-in interface
0
 enter reference line npoints and xa za xb zb
     [ Actually this is not the right order for input.
       Enter xa za first (coordinates of first point of segment)
            in order of increasing x 
       Enter xb zb next (coordinates of end point of segment)
       Enter npoints (number of points on the segment)
       NOTE: be careful that the total number of points on all interface
             segments is not more than 1000.
      ]
0 0 605000 0 201
sub menu along interface ! :      [in plot window, see interface drawn]
 1= get points on screen
 0= translate previous
 -1 =  add no points
     [ these options give the ability to add further points to this interface
       by hand]
-1
    [ ... coordinates of points on interface ...]
 line number  2
 enter line number in grid
28
new menu !!!
 enter 0 for line-point
 enter 1 to translate previous interface
 enter 2 to create a constant thickness slab
 enter 3 for line-circle-line interface
 enter 5 for circular transition from
 a horizontal to a horizontal (e.g subduction)
 geometry
 enter 4 for read-in interface
5
 enter horizontal singularity position
300000
 enter zleft zright
55588 20000
 enter slab angle
20
 enter slab length
500000    [big number]
 radius is  590111.1932229433 55588.0 20000.0
 confirm :delta_z=  35588.0
 top circle position  501829.7513867472 -20000.0
 sla1 point  300000.0 -55588.0
 enter angle left
0
 500000.0 300000.0 3.14159 353553.625139613 -200000.0
 -200000.0 -55588.0 300000.0 -55588.0
 enter points on exit-slab ,circle,in-slab
101 31 101
 enter position of top of in-slab
605000 -20000
 
    [lots of coordinates...]
 line number  3
 enter line number in grid
43
 new menu !!!
 enter 0 for line-point
 enter 1 to translate previous interface
 enter 2 to create a constant thickness slab
 enter 3 for line-circle-line interface
 enter 5 for circular transition from
 a horizontal to a horizontal (e.g subduction)
 geometry
 enter 4 for read-in interface
3
 enter top circle position. radius, angle
500000 -21000 590111.19 20
 enter slab length
500000
 500000.0 298170.2497155633 3.14159 353553.625139613 -171676.111098546
 -171676.111098546 -227597.932937419 298170.2497155633
 -56587.9998056331
 enter points on exit-slab ,circle,in-slab
50 30 50
 enter position of top of in-slab
605000 -21000
    
    [...lots of coordinates ...]

 line number  4
 enter line number in grid
55
 new menu !!!
 enter 0 for line-point
 enter 1 to translate previous interface
 enter 2 to create a constant thickness slab
 enter 3 for line-circle-line interface
 enter 5 for circular transition from
 a horizontal to a horizontal (e.g subduction)
 geometry
 enter 4 for read-in interface
2
 enter slab thickness
60000
 [ points written out ]
 z -interpolate
 now plot
 define line number that you consider as
 the top of the descending slab
43
 done
 E2 grid
 field range  .0 .0
 No variations in field
 velocity range is :  .0
 No variation in velocity field
 
     [ Plot window (xplot) showing thermal grid 
      Click on "change buttons"
      Click on "save grid"
    ]
 saving grid into file save_grid_i
 enter row number of upper slab if relevant
43
 done ! you can now rename the file
 redefine grid E 1
 E1 grid
 field range  .0 .0
 No variations in field
 velocity range is :  .0
 No variation in velocity field
     [ displays E1 grid in xplot window ]

     [ From here on, same as for SUBDUCT960826 instructions ]