# -*- Python -*- (syntax highlighting)
# ----------------------------------------------------------------------
#
# Brad T. Aagaard, U.S. Geological Survey
#
# ----------------------------------------------------------------------
#
# CUBIT journal file with geometry for TPV210.
#
# ----------------------------------------------------------------------
# Set units to SI.
# ----------------------------------------------------------------------
#{Units('si')}
#
# ----------------------------------------------------------------------
# Reset geometry.
# ----------------------------------------------------------------------
reset

# ----------------------------------------------------------------------
# Create block
# ----------------------------------------------------------------------
# Block is 36.0 km x 66.0 km x 33.0 km
# -18.0 km <= x <= 18.0 km
# -33.0 km <= y <= 33.0 km
# -33.0 km <= z <= 0.0 km
#
#{blockWidth=36.0*km}
#{blockLength=66.0*km}
#{blockHeight=33.0*km}
#
#{faultLength=30.0*km}
#{faultWidth=15.0*km}

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

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

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

# ----------------------------------------------------------------------
# Divide volumes using interface surfaces
# ----------------------------------------------------------------------
webcut volume 1 with plane surface {sI}

# ----------------------------------------------------------------------
# Split fault surface
# ----------------------------------------------------------------------
create vertex 0 {-0.5*faultLength} 0
#{vAn=Id("vertex")}
create vertex 0 {+0.5*faultLength} 0
#{vAp=Id("vertex")}
create vertex 0 {-0.5*faultLength} {-blockHeight}
#{vBn=Id("vertex")}
create vertex 0 {+0.5*faultLength} {-blockHeight}
#{vBp=Id("vertex")}
create vertex 0 {-0.5*faultLength} {-faultWidth}
#{vCn=Id("vertex")}
create vertex 0 {+0.5*faultLength} {-faultWidth}
#{vCp=Id("vertex")}
split surface 7 across location vertex {vAn} location vertex {vBn}
split surface 18 across location vertex {vAp} location vertex {vBp}
split surface 21 across location vertex {vCn} location vertex {vCp}
surface 22 name "fault_surface"
delete vertex all

#{ifdef(chopfault)}
# asperities (to ensure exact alignment with mesh)
brick x {3.0*km} y {3.0*km} z {3.0*km}
{idChopper=Id("volume")}
volume {idChopper} move x 0 y 0 z {-7.5*km}
chop volume 2 with volume {idChopper} keep
volume {idChopper} move x 0 y {7.5*km} z 0
chop volume 2 with volume {idChopper} keep
volume {idChopper} move x 0 y {-15.0*km} z 0
chop volume 2 with volume {idChopper} keep
delete volume {idChopper}
#{endif}

# ----------------------------------------------------------------------
# Imprint all volumes, merging surfaces
# ----------------------------------------------------------------------
imprint all
merge all
delete body 2
#{ifdef(chopfault)}
delete body 5 6 7 8 9 10
#{endif}







