diff --git a/configure.ac b/configure.ac
index e2e5645353..2ffe906431 100644
--- a/configure.ac
+++ b/configure.ac
@@ -109,7 +109,7 @@ AC_CHECK_HEADER([mpi.h], [], [AC_MSG_ERROR([header 'mpi.h' not found])])
 
 dnl PETSC
 AC_LANG(C)
-CIT_PATH_PETSC([3.19.4])
+CIT_PATH_PETSC([3.19.5])
 CIT_HEADER_PETSC
 CIT_CHECK_LIB_PETSC
 
@@ -254,6 +254,7 @@ AC_CONFIG_FILES([Makefile
 		tests/fullscale/linearelasticity/nofaults-3d/Makefile
 		tests/fullscale/linearelasticity/faults-2d/Makefile
 		tests/fullscale/linearelasticity/faults-3d/Makefile
+		tests/fullscale/linearelasticity/faults-3d-buried/Makefile
 		tests/fullscale/linearelasticity/greensfns-2d/Makefile
 		tests/fullscale/incompressibleelasticity/Makefile
 		tests/fullscale/incompressibleelasticity/nofaults-2d/Makefile
diff --git a/libsrc/pylith/faults/FaultCohesive.cc b/libsrc/pylith/faults/FaultCohesive.cc
index 81000c572c..fefdd41489 100644
--- a/libsrc/pylith/faults/FaultCohesive.cc
+++ b/libsrc/pylith/faults/FaultCohesive.cc
@@ -375,6 +375,10 @@ pylith::faults::FaultCohesive::createConstraints(const pylith::topology::Field&
             } // if
         } // for
     } // for
+    if (pointIS) {
+        err = ISRestoreIndices(pointIS, &points);PYLITH_CHECK_ERROR(err);
+        err = ISDestroy(&pointIS);PYLITH_CHECK_ERROR(err);
+    } // if
 
     std::vector<pylith::feassemble::Constraint*> constraintArray;
     pylith::feassemble::ConstraintSimple *constraint = new pylith::feassemble::ConstraintSimple(this);assert(constraint);
diff --git a/libsrc/pylith/meshio/DataWriter.cc b/libsrc/pylith/meshio/DataWriter.cc
index 199898dc48..545d0bc37c 100644
--- a/libsrc/pylith/meshio/DataWriter.cc
+++ b/libsrc/pylith/meshio/DataWriter.cc
@@ -248,4 +248,33 @@ pylith::meshio::DataWriter::getCoordsGlobalVec(PetscVec* coordsGlobalVec,
 } // getCoordsGlobalVec
 
 
+// ----------------------------------------------------------------------
+void
+pylith::meshio::DataWriter::_writeVec(PetscVec vector,
+                                      PetscViewer viewer) {
+    // From plexhdf5.c DMPlexGlobalVectorView_HDF5_Internal
+
+    /* To save vec in where we want, we create a new Vec (temp) with           */
+    /* VecCreate(), wrap the vec data in temp, and call VecView(temp, viewer). */
+    PetscVec temp;
+    const PetscScalar *array;
+    const char* vecName;
+    PetscLayout map;
+    PetscErrorCode err = PETSC_SUCCESS;
+
+    err = VecCreate(PetscObjectComm((PetscObject)vector), &temp);PYLITH_CHECK_ERROR(err);
+    err = PetscObjectGetName((PetscObject)vector, &vecName);PYLITH_CHECK_ERROR(err);
+    err = PetscObjectSetName((PetscObject)temp, vecName);PYLITH_CHECK_ERROR(err);
+    err = VecGetLayout(vector, &map);PYLITH_CHECK_ERROR(err);
+    err = VecSetLayout(temp, map);PYLITH_CHECK_ERROR(err);
+    err = VecSetUp(temp);PYLITH_CHECK_ERROR(err);
+    err = VecGetArrayRead(vector, &array);PYLITH_CHECK_ERROR(err);
+    err = VecPlaceArray(temp, array);PYLITH_CHECK_ERROR(err);
+    err = VecView(temp, viewer);PYLITH_CHECK_ERROR(err);
+    err = VecResetArray(temp);PYLITH_CHECK_ERROR(err);
+    err = VecRestoreArrayRead(vector, &array);PYLITH_CHECK_ERROR(err);
+    err = VecDestroy(&temp);PYLITH_CHECK_ERROR(err);
+} // _writeVec
+
+
 // End of file
diff --git a/libsrc/pylith/meshio/DataWriter.hh b/libsrc/pylith/meshio/DataWriter.hh
index 3dc8469c14..a8c95aa6a5 100644
--- a/libsrc/pylith/meshio/DataWriter.hh
+++ b/libsrc/pylith/meshio/DataWriter.hh
@@ -146,6 +146,18 @@ protected:
     void getCoordsGlobalVec(PetscVec* coordinatesVec,
                             const pylith::topology::Mesh& mesh);
 
+    /** Write PETSc Vec using viewer.
+     *
+     * This method allows us to circumvent the usual place a vector is written and write it to
+     * a desired location. See DMPlexGlobalVectorView_HDF5_Internal() in plexhdf5.c in PETSc.
+     *
+     * @param[in] vector PETSc Vec to write.
+     * @param[inout] viewer PETSc viewer to use to write vector.
+     */
+    static
+    void _writeVec(PetscVec vector,
+                   PetscViewer viewer);
+
     // NOT IMPLEMENTED //////////////////////////////////////////////////////
 private:
 
diff --git a/libsrc/pylith/meshio/DataWriterHDF5.cc b/libsrc/pylith/meshio/DataWriterHDF5.cc
index 08c92dc811..03210aaf3a 100644
--- a/libsrc/pylith/meshio/DataWriterHDF5.cc
+++ b/libsrc/pylith/meshio/DataWriterHDF5.cc
@@ -40,15 +40,6 @@
 #include <sstream> // USES std::ostringstream
 #include <stdexcept> // USES std::runtime_error
 
-extern "C" {
-    extern PetscErrorCode VecView_Seq(Vec,
-                                      PetscViewer);
-
-    extern PetscErrorCode VecView_MPI(Vec,
-                                      PetscViewer);
-
-}
-
 #if H5_VERS_MAJOR == 1 && H5_VERS_MINOR >= 8
 #define PYLITH_HDF5_USE_API_18
 #endif
@@ -207,13 +198,7 @@ pylith::meshio::DataWriterHDF5::writeVertexField(const PylithScalar t,
         err = PetscViewerHDF5SetTimestep(_viewer, istep);PYLITH_CHECK_ERROR(err);
 
         PetscVec vector = subfield.getVector();assert(vector);
-        PetscBool isseq;
-        err = PetscObjectTypeCompare((PetscObject) vector, VECSEQ, &isseq);PYLITH_CHECK_ERROR(err);
-        if (isseq) {
-            err = VecView_Seq(vector, _viewer);PYLITH_CHECK_ERROR(err);
-        } else {
-            err = VecView_MPI(vector, _viewer);PYLITH_CHECK_ERROR(err);
-        }
+        DataWriter::_writeVec(vector, _viewer);PYLITH_CHECK_ERROR(err);
         err = PetscViewerHDF5PopTimestepping(_viewer);PYLITH_CHECK_ERROR(err);
         err = PetscViewerHDF5PopGroup(_viewer);PYLITH_CHECK_ERROR(err);
 
@@ -276,13 +261,7 @@ pylith::meshio::DataWriterHDF5::writeCellField(const PylithScalar t,
         err = PetscViewerHDF5SetTimestep(_viewer, istep);PYLITH_CHECK_ERROR(err);
 
         PetscVec vector = subfield.getVector();assert(vector);
-        PetscBool isseq;
-        err = PetscObjectTypeCompare((PetscObject) vector, VECSEQ, &isseq);PYLITH_CHECK_ERROR(err);
-        if (isseq) {
-            err = VecView_Seq(vector, _viewer);PYLITH_CHECK_ERROR(err);
-        } else {
-            err = VecView_MPI(vector, _viewer);PYLITH_CHECK_ERROR(err);
-        } // if/else
+        DataWriter::_writeVec(vector, _viewer);
         err = PetscViewerHDF5PopTimestepping(_viewer);PYLITH_CHECK_ERROR(err);
         err = PetscViewerHDF5PopGroup(_viewer);PYLITH_CHECK_ERROR(err);
 
diff --git a/libsrc/pylith/meshio/DataWriterHDF5Ext.cc b/libsrc/pylith/meshio/DataWriterHDF5Ext.cc
index 6fb83543ab..b230417f5f 100644
--- a/libsrc/pylith/meshio/DataWriterHDF5Ext.cc
+++ b/libsrc/pylith/meshio/DataWriterHDF5Ext.cc
@@ -37,15 +37,6 @@
 #include <sstream> // USES std::ostringstream
 #include <stdexcept> // USES std::runtime_error
 
-extern "C" {
-    extern PetscErrorCode VecView_Seq(Vec,
-                                      PetscViewer);
-
-    extern PetscErrorCode VecView_MPI(Vec,
-                                      PetscViewer);
-
-}
-
 // ----------------------------------------------------------------------
 // Constructor
 pylith::meshio::DataWriterHDF5Ext::DataWriterHDF5Ext(void) :
@@ -195,13 +186,7 @@ pylith::meshio::DataWriterHDF5Ext::writeVertexField(const PylithScalar t,
         assert(binaryViewer);
 
         PetscVec vector = subfield.getVector();assert(vector);
-        PetscBool isseq;
-        err = PetscObjectTypeCompare((PetscObject) vector, VECSEQ, &isseq);PYLITH_CHECK_ERROR(err);
-        if (isseq) {
-            err = VecView_Seq(vector, binaryViewer);PYLITH_CHECK_ERROR(err);
-        } else {
-            err = VecView_MPI(vector, binaryViewer);PYLITH_CHECK_ERROR(err);
-        } // if/else
+        DataWriter::_writeVec(vector, binaryViewer);
 
         ExternalDataset& datasetInfo = _datasets[name];
         ++datasetInfo.numTimeSteps;
@@ -339,13 +324,7 @@ pylith::meshio::DataWriterHDF5Ext::writeCellField(const PylithScalar t,
         assert(binaryViewer);
 
         PetscVec vector = subfield.getVector();assert(vector);
-        PetscBool isseq;
-        err = PetscObjectTypeCompare((PetscObject) vector, VECSEQ, &isseq);PYLITH_CHECK_ERROR(err);
-        if (isseq) {
-            err = VecView_Seq(vector, binaryViewer);PYLITH_CHECK_ERROR(err);
-        } else {
-            err = VecView_MPI(vector, binaryViewer);PYLITH_CHECK_ERROR(err);
-        } // if/else
+        DataWriter::_writeVec(vector, binaryViewer);
 
         ExternalDataset& datasetInfo = _datasets[name];
         ++datasetInfo.numTimeSteps;
diff --git a/pylith/faults/KinSrcTimeHistory.py b/pylith/faults/KinSrcTimeHistory.py
index 4a00c7e449..72af5c776f 100644
--- a/pylith/faults/KinSrcTimeHistory.py
+++ b/pylith/faults/KinSrcTimeHistory.py
@@ -61,12 +61,12 @@ def __init__(self, name="kinsrctimehistory"):
         KinSrc.__init__(self, name)
         return
 
-    def preinitialize(self):
+    def preinitialize(self, problem):
         """Do pre-initialization setup.
         """
-        KinSrc.preinitialize(self)
+        KinSrc.preinitialize(self, problem)
 
-        ModuleKinSrc.setTimeHistory(self, self.dbTimeHistory)
+        ModuleKinSrc.setTimeHistoryDB(self, self.dbTimeHistory)
         return
 
     # PRIVATE METHODS ////////////////////////////////////////////////////
diff --git a/tests/fullscale/linearelasticity/Makefile.am b/tests/fullscale/linearelasticity/Makefile.am
index e739ab9874..cc988d7ba7 100644
--- a/tests/fullscale/linearelasticity/Makefile.am
+++ b/tests/fullscale/linearelasticity/Makefile.am
@@ -20,6 +20,7 @@ SUBDIRS = \
 	nofaults-3d \
 	faults-2d \
 	faults-3d \
+	faults-3d-buried \
 	greensfns-2d
 
 
diff --git a/tests/fullscale/linearelasticity/faults-2d/Makefile.am b/tests/fullscale/linearelasticity/faults-2d/Makefile.am
index 2fe40b9b85..ff2a499995 100644
--- a/tests/fullscale/linearelasticity/faults-2d/Makefile.am
+++ b/tests/fullscale/linearelasticity/faults-2d/Makefile.am
@@ -49,7 +49,8 @@ dist_noinst_DATA = \
 	threeblocks_ic_tri.cfg \
 	shearnoslip.cfg \
 	shearnoslip_quad.cfg \
-	shearnoslip_tri.cfg
+	shearnoslip_tri.cfg \
+	zeroslipfn.timedb
 
 
 noinst_TMP =
diff --git a/tests/fullscale/linearelasticity/faults-2d/shearnoslip.cfg b/tests/fullscale/linearelasticity/faults-2d/shearnoslip.cfg
index 2a2ba9fb3f..426e726b0c 100644
--- a/tests/fullscale/linearelasticity/faults-2d/shearnoslip.cfg
+++ b/tests/fullscale/linearelasticity/faults-2d/shearnoslip.cfg
@@ -20,12 +20,18 @@ features = [
 # ----------------------------------------------------------------------
 # faults
 # ----------------------------------------------------------------------
+[pylithapp.problem.interfaces.fault_xmid.eq_ruptures]
+rupture = pylith.faults.KinSrcTimeHistory
+
 [pylithapp.problem.interfaces.fault_xmid.eq_ruptures.rupture]
 db_auxiliary_field = spatialdata.spatialdb.UniformDB
 db_auxiliary_field.description = Fault rupture auxiliary field spatial database
 db_auxiliary_field.values = [initiation_time, final_slip_left_lateral, final_slip_opening]
 db_auxiliary_field.data = [0.0*s, 0.0*m, 0.0*m]
 
+time_history.description = Slip time function time history
+time_history.filename = zeroslipfn.timedb
+
 [pylithapp.problem.interfaces.fault_xneg.eq_ruptures.rupture]
 db_auxiliary_field = spatialdata.spatialdb.UniformDB
 db_auxiliary_field.description = Fault rupture auxiliary field spatial database
diff --git a/tests/fullscale/linearelasticity/faults-2d/zeroslipfn.timedb b/tests/fullscale/linearelasticity/faults-2d/zeroslipfn.timedb
new file mode 100644
index 0000000000..e13154a424
--- /dev/null
+++ b/tests/fullscale/linearelasticity/faults-2d/zeroslipfn.timedb
@@ -0,0 +1,10 @@
+#TIME HISTORY ascii
+TimeHistory {
+  num-points = 5
+  time-units = year
+}
+ 0.0  0.0
+ 1.0  0.0
+ 4.0  0.0
+ 8.0  0.0
+10.0  0.0
diff --git a/tests/fullscale/linearelasticity/faults-3d-buried/Makefile.am b/tests/fullscale/linearelasticity/faults-3d-buried/Makefile.am
new file mode 100644
index 0000000000..629f948d38
--- /dev/null
+++ b/tests/fullscale/linearelasticity/faults-3d-buried/Makefile.am
@@ -0,0 +1,46 @@
+# -*- Makefile -*-
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard, U.S. Geological Survey
+# Charles A. Williams, GNS Science
+# Matthew G. Knepley, University at Buffalo
+#
+# This code was developed as part of the Computational Infrastructure
+# for Geodynamics (http://geodynamics.org).
+#
+# Copyright (c) 2010-2022 University of California, Davis
+#
+# See LICENSE.md for license information.
+#
+# ----------------------------------------------------------------------
+
+TESTS = test_pylith.py
+
+dist_check_SCRIPTS = test_pylith.py
+
+dist_noinst_PYTHON = \
+	generate_gmsh.py \
+	meshes.py \
+	TestUniformSlip.py \
+	uniformslip_soln.py
+
+dist_noinst_DATA = \
+	mesh_tet.msh \
+	pylithapp.cfg \
+	uniformslip.cfg
+
+
+noinst_TMP =
+
+
+export_datadir = $(abs_builddir)
+include $(top_srcdir)/tests/data.am
+
+clean-local: clean-local-tmp clean-data
+.PHONY: clean-local-tmp
+clean-local-tmp:
+	$(RM) $(RM_FLAGS) -r output __pycache__
+
+
+# End of file
diff --git a/tests/fullscale/linearelasticity/faults-3d-buried/TestUniformSlip.py b/tests/fullscale/linearelasticity/faults-3d-buried/TestUniformSlip.py
new file mode 100644
index 0000000000..3ea7bab78e
--- /dev/null
+++ b/tests/fullscale/linearelasticity/faults-3d-buried/TestUniformSlip.py
@@ -0,0 +1,93 @@
+#!/usr/bin/env nemesis
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard, U.S. Geological Survey
+# Charles A. Williams, GNS Science
+# Matthew G. Knepley, University at Buffalo
+#
+# This code was developed as part of the Computational Infrastructure
+# for Geodynamics (http://geodynamics.org).
+#
+# Copyright (c) 2010-2022 University of California, Davis
+#
+# See LICENSE.md for license information.
+#
+# ----------------------------------------------------------------------
+
+import unittest
+
+from pylith.testing.FullTestApp import (FullTestCase, Check)
+
+import meshes
+import uniformslip_soln
+
+
+# -------------------------------------------------------------------------------------------------
+class TestCase(FullTestCase):
+
+    def setUp(self):
+        defaults = {
+            "filename": "output/{name}-{mesh_entity}.h5",
+            "exact_soln": uniformslip_soln.AnalyticalSoln(),
+            "mesh": self.mesh,
+        }
+        self.checks = [
+            Check(
+                mesh_entities=["domain"],
+                vertex_fields=[],
+                defaults=defaults,
+            ),
+            Check(
+                mesh_entities=["mat_elastic"],
+                filename="output/{name}-{mesh_entity}_info.h5",
+                cell_fields = ["density", "bulk_modulus", "shear_modulus"],
+                defaults=defaults,
+            ),
+            Check(
+                mesh_entities=["bc_xneg", "bc_xpos", "bc_yneg", "bc_ypos", "bc_zneg"],
+                filename="output/{name}-{mesh_entity}_info.h5",
+                vertex_fields=["initial_amplitude"],
+                defaults=defaults,
+            ),
+            Check(
+                mesh_entities=["fault"],
+                vertex_fields=["slip"],
+                defaults=defaults,
+            ),
+        ]
+
+    def run_pylith(self, testName, args, nprocs=1):
+        FullTestCase.run_pylith(self, testName, args, nprocs=nprocs)
+
+
+# -------------------------------------------------------------------------------------------------
+class TestTetGmsh(TestCase):
+
+    def setUp(self):
+        self.name = "uniformslip"
+        self.mesh = meshes.TetGmsh()
+        super().setUp()
+
+        TestCase.run_pylith(self, self.name, ["uniformslip.cfg"], nprocs=1)
+        return
+
+
+# -------------------------------------------------------------------------------------------------
+def test_cases():
+    return [
+        TestTetGmsh,
+    ]
+
+
+# -------------------------------------------------------------------------------------------------
+if __name__ == '__main__':
+    FullTestCase.parse_args()
+
+    suite = unittest.TestSuite()
+    for test in test_cases():
+        suite.addTest(unittest.makeSuite(test))
+    unittest.TextTestRunner(verbosity=2).run(suite)
+
+
+# End of file
diff --git a/tests/fullscale/linearelasticity/faults-3d-buried/generate_gmsh.py b/tests/fullscale/linearelasticity/faults-3d-buried/generate_gmsh.py
new file mode 100755
index 0000000000..2df4ec47c1
--- /dev/null
+++ b/tests/fullscale/linearelasticity/faults-3d-buried/generate_gmsh.py
@@ -0,0 +1,103 @@
+#!/usr/bin/env nemesis
+
+import numpy
+import gmsh
+from pylith.meshio.gmsh_utils import (VertexGroup, MaterialGroup, GenerateMesh)
+
+class App(GenerateMesh):
+    """
+    Block is DOMAIN_X by DOMAIN_Y x DOMAN_Z with discretization size DX.
+
+    -30km <= x <= +30km
+    -30km <= y <= +30km
+    -30km <= z <= 0
+    """
+    DOMAIN_X = DOMAIN_Y = 40.0e+3
+    DOMAIN_Z = 20.0e+3
+    FAULT_LENGTH = 20.0e+3
+    FAULT_WIDTH = 10.0e+3
+    DX = 5.0e+3
+
+    def __init__(self):
+        self.cell_choices = {
+            "required": False,
+            "default": "tet",
+            "choices": ["tet"],
+            }
+        self.filename = "mesh_tet.msh"
+
+    def create_geometry(self):
+        """Create geometry.
+        """
+        gmsh.initialize()
+        self.box = gmsh.model.occ.add_box(-0.5*self.DOMAIN_X, -0.5*self.DOMAIN_Y, -self.DOMAIN_Z,
+                                    self.DOMAIN_X, self.DOMAIN_Y, self.DOMAIN_Z)
+
+        self.s_xneg = 1
+        self.s_xpos = 2
+        self.s_yneg = 3
+        self.s_ypos = 4
+        self.s_zneg = 5
+        self.s_zpos = 6
+        
+        p_faulttrace_yneg = gmsh.model.occ.add_point(0.0, -0.5*self.FAULT_LENGTH, 0.0)
+        p_faulttrace_ypos = gmsh.model.occ.add_point(0.0, +0.5*self.FAULT_LENGTH, 0.0)
+        p_faultbot_yneg = gmsh.model.occ.add_point(0.0, -0.5*self.FAULT_LENGTH, -self.FAULT_WIDTH)
+        p_faultbot_ypos = gmsh.model.occ.add_point(0.0, +0.5*self.FAULT_LENGTH, -self.FAULT_WIDTH)
+
+        self.c_faulttrace = gmsh.model.occ.add_line(p_faulttrace_yneg, p_faulttrace_ypos)
+        self.c_fault_yneg = gmsh.model.occ.add_line(p_faulttrace_yneg, p_faultbot_yneg)
+        self.c_faultbot = gmsh.model.occ.add_line(p_faultbot_yneg, p_faultbot_ypos)
+        self.c_fault_ypos = gmsh.model.occ.add_line(p_faulttrace_ypos, p_faultbot_ypos)
+
+        loop_fault = gmsh.model.occ.add_curve_loop([-self.c_faulttrace, self.c_fault_yneg, self.c_faultbot, -self.c_fault_ypos])
+        self.s_fault = gmsh.model.occ.add_plane_surface([loop_fault])
+        
+        gmsh.model.occ.synchronize()
+        gmsh.model.mesh.embed(2, [self.s_fault], 3, self.box)
+        gmsh.model.mesh.embed(1, [self.c_faulttrace], 2, self.s_zpos)
+
+    def mark(self):
+        """Mark geometry for materials, boundary conditions, faults, etc.
+        """
+        materials = (
+            MaterialGroup(tag=1, entities=[self.box]),
+        )
+        for material in materials:
+            material.create_physical_group()
+
+        vertex_groups = (
+            VertexGroup(name="boundary_xneg", tag=10, dim=2, entities=[self.s_xneg]),
+            VertexGroup(name="boundary_xpos", tag=11, dim=2, entities=[self.s_xpos]),
+            VertexGroup(name="boundary_yneg", tag=12, dim=2, entities=[self.s_yneg]),
+            VertexGroup(name="boundary_ypos", tag=13, dim=2, entities=[self.s_ypos]),
+            VertexGroup(name="boundary_zneg", tag=14, dim=2, entities=[self.s_zneg]),
+            VertexGroup(name="boundary_zpos", tag=15, dim=2, entities=[self.s_zpos]),
+            VertexGroup(name="fault", tag=20, dim=2, entities=[self.s_fault]),
+            VertexGroup(name="fault_edges", tag=21, dim=1, entities=[self.c_fault_yneg, self.c_faultbot, self.c_fault_ypos]),
+        )
+        for group in vertex_groups:
+            group.create_physical_group()
+
+    def generate_mesh(self, cell):
+        """Generate the mesh. Should also include optimizing the mesh quality.
+        """
+        gmsh.option.setNumber("Mesh.MeshSizeMin", self.DX)
+        gmsh.option.setNumber("Mesh.MeshSizeMax", self.DX)
+        if cell == "hex":
+            gmsh.model.mesh.set_transfinite_automatic(recombine=True)
+            #nnodes = 3 + int(self.DOMAIN_Y/self.DX + 0.5)
+            #gmsh.model.mesh.set_transfinite_curve(self.l_xneg, nnodes)
+            #gmsh.model.mesh.set_transfinite_curve(self.l_split0, nnodes)
+            #gmsh.model.mesh.set_transfinite_surface(self.s_xmid, cornerTags=self.xmid_corners)
+            #gmsh.model.mesh.set_recombine(3, self.s_xmid)
+
+        gmsh.model.mesh.generate(3)
+        gmsh.model.mesh.optimize("Laplace2D")
+
+
+if __name__ == "__main__":
+    App().main()
+
+
+# End of file
diff --git a/tests/fullscale/linearelasticity/faults-3d-buried/mesh_tet.msh b/tests/fullscale/linearelasticity/faults-3d-buried/mesh_tet.msh
new file mode 100644
index 0000000000..fdfd2d31bb
Binary files /dev/null and b/tests/fullscale/linearelasticity/faults-3d-buried/mesh_tet.msh differ
diff --git a/tests/fullscale/linearelasticity/faults-3d-buried/meshes.py b/tests/fullscale/linearelasticity/faults-3d-buried/meshes.py
new file mode 100644
index 0000000000..46cb7ee9b8
--- /dev/null
+++ b/tests/fullscale/linearelasticity/faults-3d-buried/meshes.py
@@ -0,0 +1,42 @@
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard, U.S. Geological Survey
+# Charles A. Williams, GNS Science
+# Matthew G. Knepley, University at Buffalo
+#
+# This code was developed as part of the Computational Infrastructure
+# for Geodynamics (http://geodynamics.org).
+#
+# Copyright (c) 2010-2022 University of California, Davis
+#
+# See LICENSE.md for license information.
+#
+# ----------------------------------------------------------------------
+
+from pylith.testing.FullTestApp import MeshEntity
+
+
+class TetGmsh(object):
+    """Mesh information for tet mesh using Gmsh.
+    """
+    ENTITIES = {
+        "domain": MeshEntity(ncells=1450, ncorners=4, nvertices=425+9),
+
+        # Materials
+        "mat_elastic": MeshEntity(ncells=1450, ncorners=4, nvertices=425+9),
+
+        # Faults
+        "fault": MeshEntity(ncells=22, ncorners=3, nvertices=18),
+
+        # Boundaries
+        "bc_xneg": MeshEntity(ncells=84, ncorners=3, nvertices=55),
+        "bc_xpos": MeshEntity(ncells=84, ncorners=3, nvertices=55),
+        "bc_yneg": MeshEntity(ncells=84, ncorners=3, nvertices=55),
+        "bc_ypos": MeshEntity(ncells=84, ncorners=3, nvertices=55),
+        "bc_zneg": MeshEntity(ncells=162, ncorners=3, nvertices=98),
+        "bc_zpos": MeshEntity(ncells=164, ncorners=3, nvertices=94),
+    }
+
+
+# End of file
diff --git a/tests/fullscale/linearelasticity/faults-3d-buried/pylithapp.cfg b/tests/fullscale/linearelasticity/faults-3d-buried/pylithapp.cfg
new file mode 100644
index 0000000000..40cd56965f
--- /dev/null
+++ b/tests/fullscale/linearelasticity/faults-3d-buried/pylithapp.cfg
@@ -0,0 +1,121 @@
+[pylithapp.metadata]
+keywords = [full-scale test, 3D, box, fault, buried edges, field split conditioner, schur complement]
+features = [
+    Field split preconditioner,
+    Schur preconditioner,
+    pylith.meshio.MeshIOPetsc,
+    pylith.problems.SolnDispLagrange,
+    pylith.problems.TimeDependent,
+    pylith.meshio.DataWriterHDF5,
+    pylith.meshio.OutputSolnBoundary,
+    spatialdata.spatialdb.UniformDB
+    ]
+
+[pylithapp]
+# dump_parameters.filename = output/SIMULATION-parameters.json
+# problem.progress_monitor.filename = output/SIMULATION-progress.txt
+
+[pylithapp.launcher] # WARNING: THIS IS NOT PORTABLE
+command = mpiexec -np ${nodes}
+
+# ----------------------------------------------------------------------
+# journal
+# ----------------------------------------------------------------------
+[pylithapp.journal.info]
+timedependent = 1
+solution = 1
+petsc = 1
+meshiopetsc = 1
+isotropiclinearelasticity = 1
+dirichlettimedependent = 1
+faultcohesivekin = 1
+
+[pylithapp.journal.debug]
+#timedependent = 1
+#solution = 1
+#isotropiclinearelasticity = 1
+#dirichlettimedependent = 1
+#constraintspatialdb = 1
+#faultcohesivekin = 1
+#integratorinterface = 1
+#kinsrcstep = 1
+#outputphysics = 1
+#outputsolndomain = 1
+
+# ----------------------------------------------------------------------
+# mesh_generator
+# ----------------------------------------------------------------------
+[pylithapp.mesh_generator]
+reader = pylith.meshio.MeshIOPetsc
+
+[pylithapp.mesh_generator.reader]
+filename = mesh_tet.msh
+coordsys.space_dim = 3
+
+# ----------------------------------------------------------------------
+# problem
+# ----------------------------------------------------------------------
+[pylithapp.problem]
+defaults.quadrature_order = 1
+
+# Use nonlinear solver to ensure residual and Jacobian are consistent.
+solver = nonlinear
+
+solution = pylith.problems.SolnDispLagrange
+
+[pylithapp.problem.solution.subfields]
+displacement.basis_order = 1
+lagrange_fault.basis_order = 1
+
+[pylithapp.problem]
+solution_observers = [domain, bc_ypos]
+solution_observers.bc_ypos = pylith.meshio.OutputSolnBoundary
+
+[pylithapp.problem.solution_observers.bc_ypos]
+label = boundary_ypos
+label_value = 13
+
+# ----------------------------------------------------------------------
+# materials
+# ----------------------------------------------------------------------
+[pylithapp.problem]
+materials = [mat_elastic]
+
+[pylithapp.problem.materials.mat_elastic]
+description = Elastic material
+label_value = 1
+
+db_auxiliary_field = spatialdata.spatialdb.UniformDB
+db_auxiliary_field.description = Elastic properties
+db_auxiliary_field.values = [density, vs, vp]
+db_auxiliary_field.data = [2500*kg/m**3, 3.0*km/s, 5.2915026*km/s]
+
+
+auxiliary_subfields.density.basis_order = 0
+bulk_rheology.auxiliary_subfields.bulk_modulus.basis_order = 0
+bulk_rheology.auxiliary_subfields.shear_modulus.basis_order = 0
+
+derived_subfields.cauchy_stress.basis_order = 0
+derived_subfields.cauchy_strain.basis_order = 0
+
+
+# ----------------------------------------------------------------------
+# PETSc
+# ----------------------------------------------------------------------
+[pylithapp.problem.petsc_defaults]
+solver = True
+testing = True
+monitors = True
+
+[pylithapp.petsc]
+ksp_rtol = 1.0e-15
+ksp_atol = 1.0e-15
+ksp_max_it = 100
+ksp_gmres_restart = 50
+
+snes_atol = 1.0e-12
+
+snes_max_it = 1
+
+
+# End of file
diff --git a/tests/fullscale/linearelasticity/faults-3d-buried/test_pylith.py b/tests/fullscale/linearelasticity/faults-3d-buried/test_pylith.py
new file mode 100755
index 0000000000..73621bbc4a
--- /dev/null
+++ b/tests/fullscale/linearelasticity/faults-3d-buried/test_pylith.py
@@ -0,0 +1,51 @@
+#!/usr/bin/env nemesis
+#
+# ======================================================================
+#
+# Brad T. Aagaard, U.S. Geological Survey
+# Charles A. Williams, GNS Science
+# Matthew G. Knepley, University at Buffalo
+#
+# This code was developed as part of the Computational Infrastructure
+# for Geodynamics (http://geodynamics.org).
+#
+# Copyright (c) 2010-2022 University of California, Davis
+#
+# See LICENSE.md for license information.
+#
+# ======================================================================
+
+from pylith.testing.FullTestApp import TestDriver, FullTestCase
+
+import unittest
+
+
+class TestApp(TestDriver):
+    """Driver application for full-scale tests.
+    """
+
+    def __init__(self):
+        """Constructor.
+        """
+        TestDriver.__init__(self)
+        return
+
+    def _suite(self):
+        """Create test suite.
+        """
+        suite = unittest.TestSuite()
+
+        import TestUniformSlip
+        for test in TestUniformSlip.test_cases():
+            suite.addTest(unittest.makeSuite(test))
+
+        return suite
+
+
+# ----------------------------------------------------------------------
+if __name__ == '__main__':
+    FullTestCase.parse_args()
+    TestApp().main()
+
+
+# End of file
diff --git a/tests/fullscale/linearelasticity/faults-3d-buried/uniformslip.cfg b/tests/fullscale/linearelasticity/faults-3d-buried/uniformslip.cfg
new file mode 100644
index 0000000000..64c9ffca5f
--- /dev/null
+++ b/tests/fullscale/linearelasticity/faults-3d-buried/uniformslip.cfg
@@ -0,0 +1,108 @@
+[pylithapp.metadata]
+# There is no analytical solution for this problem. We just check ranges of values in the output.
+#
+description = Prescribed uniform left-lateral fault slip with Dirichlet boundary conditions.
+authors = [Brad Aagaard]
+keywords = [fault, prescribed slip, uniform slip, tetrahedral cells]
+version = 1.0.0
+pylith_version = [>=3.0, <4.0]
+base = [pylithapp.cfg]
+arguments = [uniformslip.cfg]
+
+features = [
+    Static simulation,
+    Field split preconditioner,
+    pylith.faults.FaultCohesiveKin,
+    pylith.faults.KinSrcStep,
+    pylith.materials.Elasticity,
+    pylith.materials.IsotropicLinearElasticity,
+    pylith.bc.DirichletTimeDependent,
+    spatialdata.spatialdb.UniformDB
+    ]
+
+# ----------------------------------------------------------------------
+# output
+# ----------------------------------------------------------------------
+[pylithapp]
+dump_parameters.filename = output/uniformslip-parameters.json
+problem.progress_monitor.filename = output/uniformslip-progress.txt
+
+problem.defaults.name = uniformslip
+
+# ----------------------------------------------------------------------
+# solution
+# ----------------------------------------------------------------------
+[pylithapp.problem.solution.subfields.displacement]
+basis_order = 1
+
+# ----------------------------------------------------------------------
+# fault
+# ----------------------------------------------------------------------
+[pylithapp.problem]
+interfaces = [fault]
+
+[pylithapp.problem.interfaces.fault]
+label = fault
+label_value = 20
+
+edge = fault_edges
+edge_value = 21
+
+observers.observer.data_fields = [slip, lagrange_multiplier_fault]
+
+[pylithapp.problem.interfaces.fault.eq_ruptures.rupture]
+db_auxiliary_field = spatialdata.spatialdb.UniformDB
+db_auxiliary_field.description = Fault rupture auxiliary field spatial database
+db_auxiliary_field.values = [initiation_time, final_slip_left_lateral, final_slip_reverse, final_slip_opening]
+db_auxiliary_field.data = [0.0*s, +2.0*m, 0.0*m, 0.0*m]
+
+
+# ----------------------------------------------------------------------
+# boundary conditions
+# ----------------------------------------------------------------------
+[pylithapp.problem]
+bc = [bc_xneg, bc_xpos, bc_yneg, bc_ypos, bc_zneg]
+bc.bc_xneg = pylith.bc.DirichletTimeDependent
+bc.bc_xpos = pylith.bc.DirichletTimeDependent
+bc.bc_yneg = pylith.bc.DirichletTimeDependent
+bc.bc_ypos = pylith.bc.DirichletTimeDependent
+bc.bc_zneg = pylith.bc.DirichletTimeDependent
+
+[pylithapp.problem.bc.bc_xneg]
+label = boundary_xneg
+label_value = 10
+constrained_dof = [0]
+db_auxiliary_field = pylith.bc.ZeroDB
+db_auxiliary_field.description = Dirichlet BC -x edge
+
+[pylithapp.problem.bc.bc_xpos]
+label = boundary_xpos
+label_value = 11
+constrained_dof = [0]
+db_auxiliary_field = pylith.bc.ZeroDB
+db_auxiliary_field.description = Dirichlet BC +x edge
+
+[pylithapp.problem.bc.bc_yneg]
+label = boundary_yneg
+label_value = 12
+constrained_dof = [1]
+db_auxiliary_field = pylith.bc.ZeroDB
+db_auxiliary_field.description = Dirichlet BC -y edge
+
+[pylithapp.problem.bc.bc_ypos]
+label = boundary_ypos
+label_value = 13
+constrained_dof = [1]
+db_auxiliary_field = pylith.bc.ZeroDB
+db_auxiliary_field.description = Dirichlet BC +y edge
+
+[pylithapp.problem.bc.bc_zneg]
+label = boundary_zneg
+label_value = 14
+constrained_dof = [2]
+db_auxiliary_field = pylith.bc.ZeroDB
+db_auxiliary_field.description = Dirichlet BC -z edge
+
+
+
+# End of file
diff --git a/tests/fullscale/linearelasticity/faults-3d-buried/uniformslip_soln.py b/tests/fullscale/linearelasticity/faults-3d-buried/uniformslip_soln.py
new file mode 100644
index 0000000000..1a9d30c93e
--- /dev/null
+++ b/tests/fullscale/linearelasticity/faults-3d-buried/uniformslip_soln.py
@@ -0,0 +1,96 @@
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard, U.S. Geological Survey
+# Charles A. Williams, GNS Science
+# Matthew G. Knepley, University at Buffalo
+#
+# This code was developed as part of the Computational Infrastructure
+# for Geodynamics (http://geodynamics.org).
+#
+# Copyright (c) 2010-2022 University of California, Davis
+#
+# See LICENSE.md for license information.
+#
+# ----------------------------------------------------------------------
+#
+# @brief Uniform slip on a buried fault (no analytical solution).
+#
+# We just check info fields.
+
+import numpy
+
+# Physical properties
+p_density = 2500.0
+p_vs = 3000.0
+p_vp = 5291.502622129181
+
+p_mu = p_density * p_vs**2
+p_lambda = p_density * p_vp**2 - 2 * p_mu
+
+# ----------------------------------------------------------------------
+class AnalyticalSoln(object):
+    """Analytical solution to axial/shear displacement problem.
+    """
+    SPACE_DIM = 3
+    TENSOR_SIZE = 6
+
+    def __init__(self):
+        self.fields = {
+            "density": self.density,
+            "shear_modulus": self.shear_modulus,
+            "bulk_modulus": self.bulk_modulus,
+            "initial_amplitude": self.zero_displacement,
+            "slip": self.slip,
+        }
+        return
+
+    def getField(self, name, mesh_entity, pts):
+        return self.fields[name](pts)
+
+    def getMask(self, name, mesh_entity, pts):
+        mask = None
+        if name == "displacement":
+            fdist = numpy.dot(pts-fault_center, fault_normal)
+            index = numpy.argmin(numpy.abs(numpy.min(fdist)))
+            mask = numpy.abs(fdist) < 1.0
+        return mask
+    
+    def zero_displacement(self, locs):
+        """Compute displacement field at locations.
+        """
+        (npts, dim) = locs.shape
+
+        disp = numpy.zeros((1, npts, self.SPACE_DIM), dtype=numpy.float64)
+        return disp
+
+    def density(self, locs):
+        """Compute density field at locations.
+        """
+        (npts, dim) = locs.shape
+        density = p_density * numpy.ones((1, npts, 1), dtype=numpy.float64)
+        return density
+
+    def shear_modulus(self, locs):
+        """Compute shear modulus field at locations.
+        """
+        (npts, dim) = locs.shape
+        shear_modulus = p_mu * numpy.ones((1, npts, 1), dtype=numpy.float64)
+        return shear_modulus
+
+    def bulk_modulus(self, locs):
+        """Compute bulk modulus field at locations.
+        """
+        (npts, dim) = locs.shape
+        bulk_modulus = (p_lambda + 2.0 / 3.0 * p_mu) * numpy.ones((1, npts, 1), dtype=numpy.float64)
+        return bulk_modulus
+
+    def slip(self, locs):
+        """Compute slip field at locations.
+        """
+        (npts, dim) = locs.shape
+        slip = numpy.zeros((1, npts, 3), dtype=numpy.float64)
+        slip[:,:,1] = +2.0
+        return slip
+
+
+# End of file
diff --git a/tests/mmstests/linearelasticity/faults-2d/ThreeBlocksStatic.cc b/tests/mmstests/linearelasticity/faults-2d/ThreeBlocksStatic.cc
index f1351501ef..34fa9dd9fd 100644
--- a/tests/mmstests/linearelasticity/faults-2d/ThreeBlocksStatic.cc
+++ b/tests/mmstests/linearelasticity/faults-2d/ThreeBlocksStatic.cc
@@ -213,6 +213,7 @@ class pylith::_ThreeBlocksStatic {
         data->journalName = "ThreeBlocksStatic";
 
         data->isJacobianLinear = true;
+        data->allowZeroResidual = true;
 
         data->meshFilename = ":UNKNOWN:"; // Set in child class.