[pylithapp]

# ----------------------------------------------------------------------
# PROBLEM DESCRIPTION
# ----------------------------------------------------------------------

#
# This is a purely elastic static problem using Neumann (traction)
# boundary conditions. We apply normal tractions to the top surface
# and roller (fixed normal but free lateral motion) boundary
# conditions on the lateral sides and bottom surfaces.

# ----------------------------------------------------------------------
# RUNNING THE SIMULATON
# ----------------------------------------------------------------------

# This is not a self-contained simulation configuration file. This
# file only specifies parameters specific to tutorial step18.
# The general parameters are specificed in the pylithapp.cfg
# file which PyLith reads by default.
#
# To run the simulation:
# pylith step18.cfg
#
# Output will be directed to directory output.

# ----------------------------------------------------------------------
# problem
# ----------------------------------------------------------------------
[pylithapp.timedependent.formulation.time_step]
# Define the total time for the simulation and the default time step size.
total_time = 0.0*s ; total time of simulation
# Define an appropriat time step for simulations. Important for
# nondimensionalization of velocities and slip rates.
dt = 5.0*year

[pylithapp.timedependent]
# Set bc to an array of 6 boundary conditions (one for each side of the domain).
bc = [x_pos,x_neg,y_pos,y_neg,z_neg,z_pos]

[pylithapp.timedependent.implicit]
# Set the output to an array of 2 output managers.
# We will output the solution over the domain and the ground surface.
output = [domain,subdomain]

# Set subdomain component to OutputSolnSubset (subset of domain).
output.subdomain = pylith.meshio.OutputSolnSubset

# ----------------------------------------------------------------------
# boundary conditions
# ----------------------------------------------------------------------
# Set the parameters for Dirichlet boundary conditions applied on the
# +x, -x, +y, -y, and -z faces of the box, and the Neumann boundary
# conditions applied on the +z face of the box.
#
# We fix the x DOF on the +x and -x faces, the y DOF on the +y and -y
# faces, and the z DOF on the bottom (-z) face. We use the ZeroDispDB
# (default) since we want zero displacements.
#
# We apply axial tractions on the +z face.

# The label corresponds to the name of the nodeset in CUBIT.

# +x face
[pylithapp.timedependent.bc.x_pos]
label = face_xpos
bc_dof = [0]
db_initial.label = Dirichlet BC on +x

# -x face
[pylithapp.timedependent.bc.x_neg]
label = face_xneg
bc_dof = [0]
db_initial.label = Dirichlet BC on -x

# +y face
[pylithapp.timedependent.bc.y_pos]
label = face_ypos
bc_dof = [1]
db_initial.label = Dirichlet BC on +y

# -y face
[pylithapp.timedependent.bc.y_neg]
label = face_yneg
bc_dof = [1]
db_initial.label = Dirichlet BC on -y

# -z face
[pylithapp.timedependent.bc.z_neg]
label = face_zneg
bc_dof = [2]
db_initial.label = Dirichlet BC on -z

# +z face -- change bc type to Neumann
[pylithapp.timedependent.bc]
z_pos = pylith.bc.Neumann

[pylithapp.timedependent.bc.z_pos]
label = face_zpos

db_initial = spatialdata.spatialdb.SimpleDB
db_initial.label = Neumann BC on +z
db_initial.iohandler.filename = spatialdb/tractions_axial_pressure.spatialdb.txt
# Use linear interpolation
db_initial.query_type = linear

output.cell_info_fields = [initial_value]
output.writer.filename = output/step18-traction.vtk
output.cell_filter = pylith.meshio.CellFilterAvgSubMesh

# We must specify quadrature information for the cell faces.
quadrature.cell = pylith.feassemble.FIATLagrange
quadrature.cell.dimension = 2
quadrature.cell.quad_order = 2

# Because normal for +z surface is [0,0,1], the horizontal and
# vertical shear directions are ambiguous. We provide a "fake" up
# direction of [0,1,0] so that the horizontal shear direction ("up" x
# normal) is [1,0,0] and the vertical shear direction (normal x horiz
# shear dir) is [0,1,0].
up_dir = [0,1,0]

# ----------------------------------------------------------------------
# output
# ----------------------------------------------------------------------
# Give basename for VTK domain output of solution over domain.
[pylithapp.problem.formulation.output.domain.writer]
filename = output/step18.vtk

# Give basename for VTK domain output of solution over ground surface.
[pylithapp.problem.formulation.output.subdomain]
label = face_zpos ; Name of nodeset for subdomain
writer.filename = output/step18-groundsurf.vtk

# Give basename for VTK output of upper_crust state variables.
[pylithapp.timedependent.materials.upper_crust.output]
cell_filter = pylith.meshio.CellFilterAvgMesh
writer.filename = output/step18-upper_crust.vtk

# Give basename for VTK output of lower_crust state variables.
[pylithapp.timedependent.materials.lower_crust.output]
cell_filter = pylith.meshio.CellFilterAvgMesh
writer.filename = output/step18-lower_crust.vtk