input.3Dpunch

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!              Input File                                                      |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

This file is read by two subroutines in the code:
- read_controlling_parameters
- read_input_file
The indications between [] indicate the type of the read parameter. It can be
an integer [int], a real*8 [dp], a character chain [char] or a boolean [bool].


CONTROLLING PARAMETERS
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-

[int]debug is a parameter that switches on/off various prints and outputs
(the level of printing for error, warning and messages of the solver id%cntl(4)
is set to the debug value).
if debug is zero, no debugging ; if debug is equal to one, this triggers the terminal 
display of key parameters ; if debug is equal to two, same as debug equal to one, 
and the code produces various vtk files in the DEBUG subdirectoris.

      debug = 1

[bool]doDoRuRe is a flag that triggers the production of output files needed to produce
the DoRuRes. DoRuRe stands for 'Douar Run Report'. 

       doDoRuRe = T

[bool]compute_qp_gram triggers the production of qpgrams for every grid.

      compute_qpgram = F

RESTART
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-

[int]irestart is a restart flag; if irestart is not 0, the run will restart from
an output file given by [char]restartfile and at step.

      irestart = 0

      restartfile=OUT/time_0001.txt


TIMESTEPPING
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-

[dp]dt is the time step length (if dt is negative, courant conidition is used and
automatic time stepping is turned on)

      dt=0.00065d0

[int]nstep is the number of time steps

      nstep = 1

[dp]courant is only used when dt is negative; it determines the size of the time
step from the maximum value of the velocity field amplitude. The time step s
the product of courant by the ratio of the smallest leaf size by the maximum
velocity.

      courant=.5d0

[bool]normaladvect is a flag used to determine which algorithm to use to calculate
the new geoletry of the normals to the surfaces at the nodes on the surfaces
if normaladvect is T, the normals are advected using the velocity gradient
if normaladvect is F, the normals are re-computed from the geometry of the
surface

      normaladvect = T

GRID ITERATIONS
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-

[int]griditer is a flag that allows for nonlinear iterations; when positive, a fixed number
(griditer) of iterations is permitted; when negative, the number of
nonlinear iterations is determined by a convergence criterion.

      griditer = -80

[dp]octree_refine_ratio is the threshold value used to determine whether the octree 
has converged or not. the larger the value, the less stringent the test.

      octree_refine_ratio=.00005d0


NONLINEAR ITERATIONS
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-


[int]nonlinear_iterations is the maximum number of nonlinear iterations (i.e. the 
iterations on a given constant grid)
if nonlinear_iterations is positive, it simply is the number of nonlinear
iterations performed for each grid. When negative it indicates an upper bound
of nonlinear iterations, but the actual number of nonlinear iterations is
determined by a convergence criterion (see the 'tol' parameter)

      nonlinear_iterations = -200

[dp]tol is the relative tolerance used to estimate convergence on the computed
velocity field

      tol=0.0001d0

[bool]adaptive_tol is a flag that toggles on/off the evolution of the tol parameter
with the grid level: when velocity convergence is reached on a grid, the
latest meaure of the velocity difference between the two last obtained
solutions is put in tol, so that on the following generated grid, the solution
reaches at least the same level of convergence. It allows to start with a not
too stringent value of tol at uniform octree level that evolves with every
grid, assuming that increasing the level of refinement of the octree allows to
better capture the solution, hence allowing a tighter convergence. 

      adaptive_tol = F

OCTREES
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-

[int]leveluniform_oct is the level of uniform discretization of space; note that a
level is a power of two used to divide the unit cube

      leveluniform_oct = 5

[int]levelmax_oct is maximum level of octree discretization

      levelmax_oct = 8

[bool]ismooth is a flag to impose an additional level of smoothing after refinement
for the surfaces and strain rate. It ensures that no leaf is flanked by
other leaves differing by more than 1 level of refinement

      ismooth = F

[int]noctreemax is the maximum size of any octree used in all computations

      noctreemax=10000000

[dp]refine_ratio is used to determine octree refinement based on a given criterion. 
All leaves where the criterion is larger than refine_ratio times the maximum 
of this criterion are refined 

      refine_ratio=0.07d0

[int]refine_criterion determines which refinement algorithm is to be used.
Several criteria exist for the refinement of the osolve octree. 1 is the second 
invariant of the deviatoric strain-rate tensor; 2 is the sum of the squares of the
diagonal terms of the deviatoric strain-rate tensor; 3 is the second invariant
of the deviatoric strain rate tensor timses the leaf size. any other value sets
the criterion to zero and leads to no refinement. 

      refine_criterion = 1

[int]initial_refine_level is the initial level at which the refinement of the
octree will be performed. it has to be smaller than levelmax_oct

      initial_refine_level = 5

[bool]renumber_nodes is a flag that can toggle on/off the renumbering of the nodes
by mean of Sloan's algorithm (T/F)

      renumber_nodes = F

PRESSURE
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
[int]smoothing_type is a parameter allows to choose which type of smoothing is to 
be applied to the pressure field: 0 is none, 1 is center->nodes->center, 2 is the 
same, but weighted by neighbouring elemental volumes, 3 is regular grid+SPH, and 4 is SPH.  

      smoothing_type = 1 

CLOUD
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-

[int]npmin and [int]npmax are used to update the 3D volumetric cloud. npmin corresponds
to the minimum number of particles in any leaf; npmax is the maximum allowable
number in any leaf

      npmin=1
      npmax=4


FEM + DIVFEM + MUMPS
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-

[int]levelcut is the number of levels used to estimate the volume integrals in the
divFEM algorithmi for cut cells; by testing, it has been estimated that a value
of 2 is very accurate

      levelcut=2

[int]levelapprox is the number of levels used to estimate the remaining integrals
using an improved version of Marthijn's clever algorithm, usually 3 is plenty

      levelapprox=3

[dp]penalty is a global penalty parameter used to impose the bad faces or
incompatible faces linear constraints

      penalty=1.d8

[bool]excl_vol is a parameter that toggles off the assumption that lsf's are built
on top of one another

      excl_vol= F


TEMPERATURE
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
[bool] calculate_temp is a logical flag that determines whether the temperature
calculations are performed (T) or not (F)
[dp] ztemp is the height interval between which a linear temperature gradient 
is set: the temperature is 1 at the bottom, and 0 at ztemp.
[dp] tempscale

      calculate_temp = F

      ztemp = .713d0

      tempscale=300.d0

MATERIALS
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-

[int]nmat is number of materials

      nmat=1

[int]material0 determines what is the material above the 1st surface (the free
surface) if material0 is 0 then it is the void and the properties "0" are used
for that part of the model; otherwise the material is one of the material,
comprised between 1 and nmat

      material0 =  1

[dp]densityi, [dp]viscosityi and [dp]penaltyi are the density, viscosity and
incompressibility used for material i; there should be nmat sets of material
properties; there should also be a nil material if material0 has been set to 0
[dp]expon is the nonlinear viscosity exponent
[dp]diffusivity is the heat diffusivity
[dp]heat is the heat production

[char]plasticity_type is the type of plasticity
- No    -> no plasticity, purely (nonlinear) viscous material
- vM    -> von Mises yield criterion
        -> [dp]plasticity_1st_param is the yield value
- DPI   -> Drucker-Prager type of yield criterion
        -> the yield locus passes through the outer apices of the Mohr-Coulomb hexagon
        -> [dp]plasticity_1st_param is the angle phi
        -> [dp]plasticity_2nd_param is the cohesion c
- DPII  -> Drucker-Prager type of yield criterion
        -> the yield locus passes through the inner apices of the Mohr-Coulomb hexagon
        -> [dp]plasticity_1st_param is the angle phi
        -> [dp]plasticity_2nd_param is the cohesion c
- DPIII -> Drucker-Prager type of yield criterion
        -> [dp]plasticity_1st_param is alpha
        -> [dp]plasticity_2nd_param is k
- MC    -> Mohr-Coulomb type of yield criterion
        -> [dp]plasticity_1st_param is the angle phi
        -> [dp]plasticity_2nd_param is the cohesion c

      density0              = 0.d0
      viscosity0            = 1.d-8
      penalty0              = 1.d8
      expon0                = 1.d0
      diffusivity0          = 1.d0
      heat0                 = 0.d0
      activationenergy0     = 0.d0
      plasticity_type0      = No

      density1              = 0.d0
      viscosity1            = 10000.d0
      penalty1              = 1.d8
      expon1                = 1.d0
      diffusivity1          = 1.d0
      heat1                 = 0.d0
      activationenergy1     = 0.d0
      plasticity_type1      = vM
      plasticity_1st_param1 = 1.d0

SURFACES
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-

[int]ns is number of surfaces to track

      ns=1

for each surface, one needs to define a levelt, itype, material, surface_type, activation_time 
and surface_params .
- [int]levelt is the inital level for the particles on the surface; to be accurate
and avoid wholes in the surface during definition of the lsf, one should use
levelt eq to levelmax_oct+1 for all surfaces as a minimum value; 
- [int]itype should be 1 for foldable surfaces or 0 for nonfoldable surfaces; 
- [int]material is the material type refering to the table of material available (max nmat); 
- [dp]activation_time is the time the surface becomes active (before that time, it
  is glued to the 0th surface). this parameter is useful when defining stratigraphic horizons; 
  default is -1, ie the surface is not glued to the free surface
- [int]surface_type is comprised between 1 and 8. 
  1 corresponds to  a flat surface, 
    -> surface_param_01 is the z level
  2 to a rectangular emboss, 
    -> surface_param_01 is the z level
    -> surface_param_02 and 03 are x1,x2 
    -> surface_param_04 and 05 are y1,y2 
    -> surface_param_06 is the thickness 
  3 to a convex spherical emboss, 
    -> surface_param_01 is the z level
    -> surface_param_02 and 03 are x0,y0 
    -> surface_param_04 is the radius 
  4 to concave spherical emboss,
    -> surface_param_01 is the z level
    -> surface_param_02 and 03 are x0,y0 
    -> surface_param_04 is the radius 
  5 to a double rectangular emboss, 
    -> surface_param_01 is the z level
    -> surface_param_02 and 03 are x1,x2 
    -> surface_param_04 and 05 are x3,x4 
    -> surface_param_06 and 07 are y1,y2 
    -> surface_param_08 and 09 are y3,y4 
    -> surface_param_10 is the thickness 
  6 to a sinus, 
    -> surface_param_01 is the z level
    -> surface_param_02 is the wavelength 
    -> surface_param_03 is the amplitude 
  7 to a noisy surface, 
    -> surface_param_01 is the z level
    -> surface_param_02 is the noise amplitude
  8 to a double sinus.
    -> surface_param_01 is the z level
    -> surface_param_02 is the x-wavelength 
    -> surface_param_03 is the x-amplitude 
    -> surface_param_04 is the y-wavelength 
    -> surface_param_05 is the y-amplitude 
- [int]leveloct is the level at which the octree will be refined in the vicinity of the surface.

      levelt1             = 7
      itype1              = 0
      surface_type_1      = 1
      rand1               = F 
      surface_param_01_1  = 0.415432
      material1           = 1
      activation_time_1   = -1. 
      leveloct1           = 5

[int]niter_move is the number of iterations used to update particle positions
using an implicit, mid-point algorithm (default is 10)

      niter_move = 10

[dp]stretch is the maximum allowed increase in linear length between two initially
adjacent particles on any surface; when this stretch is achieved, a new
particle is inserted on the surface, half-way along the stretched edge

      stretch = 1.5d0


[int]criterion is criterion used to define the octree in the vicinity of the
sufaces; criterion 1 corresponds to imposing that all leaves that contain at
least one particle of any surface is at levelmax_oct; criterion 2 corresponds
to imposing that discretization is proportional to the curvature of the
surface; curvature is calculated from the local divergence of the normals.

      criterion = 3

[dp]anglemax is only defined for criterion 2; t is the maximum allowable angle
between two adjacent normals; if the angle is greater than nglemax, the local
octree leaves are forced to be at level levelmax_oct; otherwise they are
proportionally larger (smaller levels) (default is 10)

     anglemax = 30.d0


REFINEMENT IN BOXES
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-

[int]nboxes is the number of boxes in which the user imposes a set level of
discretization

      nboxes=0

for each box we need to specify the two end corners and the level
the syntax is [dp]boxnx0, [dp]boxnx1, [dp]boxny0, [dp]boxny1, [dp]boxnz0, [dp]boxnz1
and boxnlevel where n is the box number

      box1x0=.34376d0
      box1x1=.65624d0
      box1y0=0.34376d0
      box1y1=.65625d0
      box1z0=.0d0
      box1z1=.078125d0
      box1level=8


REFINEMENT ON CUBE FACES
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
[bool]ref_on_faces toggles on/off the user imposed refinement on faces.
for each of the six faces of the cube, on defines the level at which the
desired area is to be refined. this area is given by bottom, top, left and
right coordinates that are comprised between 0 and 1
faces 1,2,3,4,5,6 respectively correspond to x=0,x=1,y=0,y=1,z=0,z=1

      ref_on_faces = F

      level_face1=5
      b1=.02
      t1=.51
      l1=.11
      r1=.81
      level_face2=5
      b2=.45
      t2=.55
      l2=.46
      r2=.56
      level_face3=5
      b3=.0
      t3=.2
      l3=.3
      r3=.7
      level_face4=5
      b4=.4
      t4=.5
      l4=.14
      r4=.4
      level_face5=8
      b5=0.4
      t5=0.6
      l5=0.4
      r5=0.6
      level_face6=5
      b6=.26
      t6=.56
      l6=.16
      r6=.86


EROSION
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-

[bool]erosion is a flag that toggles on/off the erosion. 

      erosion = F

if erosion is on, one also needs to set the erosion level/height, [dp]zerosion
Note that this is a first attempt at erosion; in future versions, DOUAR
should be easily linked to a surface processes model like CASCADE

      zerosion=.81d0


MATRIX VISUALISATION
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-

[bool]visualise_matrix allows the user to turn on the visual representation of the 
matrices used in the code. Be careful, the generated postscript files are huge !

      visualise_matrix = F

CROSS SECTIONS
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-

[int]nsections is the number of cross-sections to be output. [dp]scale is the scale used to produce the postscript file.
[int]xyz takes values 1, 2 or 3, and corresponds to planes defined by x=constant, y=constant, and z=constant resp.
The [bool]flags are self explanatory. [char]colormap is the chosen colormap (jet or hot). [int]ncolours is the number of
colours used to produce the plot. 

nsections = 3 

xyz_1          = 1 
slice_1        = 0.5001
flag_press_1   = T
flag_spress_1  = T 
flag_e2d_1     = T
flag_e3d_1     = F
flag_lode_1    = F 
flag_crit_1    = F 
flag_grid_1    = T 
flag_mu_1      = T   
flag_u_1       = F
flag_v_1       = F
flag_w_1       = F
flag_q_1       = F 
flag_uvw_1     = F 
flag_lsf_1     = F 
flag_vfield_1  = F 
flag_colour_1  = T
flag_plastic_1 = F
flag_velvect_1 = T
scale_1        = 500.
colormap_1     = jet
ncolours_1     = 256

xyz_2          = 2
slice_2        = 0.5001
flag_press_2   = T
flag_spress_2  = T
flag_e2d_2     = T
flag_e3d_2     = F 
flag_lode_2    = F 
flag_crit_2    = F 
flag_grid_2    = T 
flag_mu_2      = T   
flag_u_2       = F
flag_v_2       = F
flag_w_2       = F
flag_q_2       = F 
flag_uvw_2     = T 
flag_lsf_2     = F 
flag_vfield_2  = F 
flag_colour_2  = T
flag_plastic_2 = F
flag_velvect_2 = F
scale_2        = 500.
colormap_2     = jet 
ncolours_2     = 256


xyz_3          = 3
slice_3        = 0.0001
flag_press_3   = T 
flag_spress_3  = T 
flag_e2d_3     = T
flag_e3d_3     = F
flag_lode_3    = F 
flag_crit_3    = F 
flag_grid_3    = T 
flag_mu_3      = T   
flag_u_3       = F
flag_v_3       = F
flag_w_3       = F
flag_q_3       = F
flag_uvw_3     = T 
flag_lsf_3     = F 
flag_vfield_3  = F 
flag_colour_3  = T
flag_plastic_3 = F
flag_velvect_3 = F
scale_3        = 500.
colormap_3     = jet
ncolours_3     = 256