# -*- Python -*- (syntax highlighting)
#
# Description:
#
# This journal creates defines the geometry for a mesh with dimensions
# of 2000 km (x), 1000 km (y), and 400 km (z).  The volume is divided
# into 2 horizontal layers:
#   (1) a top layer from 0 to 40 km, marking the surface elastic
#   layer;
#   (2) a bottom layer from 40 to 400 km, which approximates the
# viscoelastic half-space.
#
# The fault plane is also cut horizontally at 20 km depth,
# representing the division between the upper locked section and the
# lower creeping section.  The volume is cut vertically by the fault
# surface, which lies along the x-axis, and by a symmetry plane, which
# lies along the y-axis.  The fault is restricted to the upper two
# (elastic) layers.
#
# ----------------------------------------------------------------------
# Set units to SI.
# ----------------------------------------------------------------------
${Units('si')}
#
# ----------------------------------------------------------------------
# Reset geometry.
# ----------------------------------------------------------------------
reset

# Make sure undo is off to prevent errors in stitching volumes.
undo off

# ----------------------------------------------------------------------
# Create outer block
# ----------------------------------------------------------------------
# Block is 2000 km x 1000 km x 400 km
# -1000 km <= x <= 1000 km
# -500 km <= y <= 500 km
# -400 km <= z <= 0 km
#
${blockWidth=1000.0*km}
${blockLength=2000.0*km}
${blockHeight=400.0*km}
${faultWidth=40.0*km}
${lockingDepth=20.0*km}

brick x {blockLength} y {blockWidth} z {blockHeight}

# Translate block so the top is at z=0
volume 1 move x {-xoffset} y 0 z {-0.5*blockHeight}

# ----------------------------------------------------------------------
# Create interface surfaces
# ----------------------------------------------------------------------
create planar surface with plane xplane offset 0
${sid=Id("surface")}
surface {sid} name "fault_surface"

create planar surface with plane zplane offset {-faultWidth}
${sid=Id("surface")}
surface {sid} name "material_interface"

create planar surface with plane yplane offset 0
${sid=Id("surface")}
surface {sid} name "symmetry_plane"

# ----------------------------------------------------------------------
# Cut volume with interfaces.
# ----------------------------------------------------------------------
# Cut fault with locking transition (not needed in volume)
webcut volume 1 with plane fault_surface
webcut volume 1 5 with plane material_interface
webcut volume 1 5 6 7 with plane symmetry_plane

# Inscribe locking transition on fault plane
brick x {5.0*km} y {blockWidth+10.0*km} z {2.0*lockingDepth}
${vid=Id("volume")}
chop volume 2 with volume {vid}

# Inscribe bottom of fault on fault plane
brick x {5.0*km} y {blockWidth+10.0*km} z {4.0*faultWidth}
${vid=Id("volume")}
chop volume 14 with volume {vid}

# ----------------------------------------------------------------------
# Imprint all volumes, merging surfaces.
# ----------------------------------------------------------------------
imprint all
merge all
delete body 3 4 13 16 17
compress all

