Output of a solution subfield for fault
## Describe the bug

Output of a solution subfield for the fault triggers a PETSc error.

## To Reproduce

**PETSc branch**: `knepley/pylith`
**PyLith branch**: `baagaard/fix-output-fault-solution-field`
```
cd examples/strikeslip-2d
pylith step01_slip.cfg
```

## Error message

```
0 TS dt 0.01 time 0.
    0 SNES Function norm 4.895713226482e-02 
    Linear solve converged due to CONVERGED_ATOL iterations 35
    1 SNES Function norm 2.540378711183e-12 
  Nonlinear solve converged due to CONVERGED_FNORM_ABS iterations 1
[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR: Null argument, when expecting valid pointer
[0]PETSC ERROR: Null Pointer: Parameter # 1
[0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
[0]PETSC ERROR: Petsc Development GIT revision: v3.17.4-1194-g766cdfe255  GIT Date: 2022-09-02 04:19:45 +0000
[0]PETSC ERROR: /Users/baagaard/software/unix/py310-venv/pylith-debug/bin/mpinemesis on a arch-clang-13.1_debug named IGSKCI164LM006 by baagaard Wed Oct 26 10:50:15 2022
[0]PETSC ERROR: Configure options --PETSC_ARCH=arch-clang-13.1_debug --with-debugging=1 --with-clanguage=c --with-mpi-compilers=1 --with-shared-libraries=1 --with-64-bit-points=1 --with-large-file-io=1 --with-lgrind=0 --download-chaco=1 --download-parmetis=1 --download-metis=1 --download-triangle --download-ml=1 --download-superlu=1 --with-fc=0 --download-f2cblaslapack --with-hdf5=1 --with-hdf5-include=/Users/baagaard/software/unix/hdf5-1.12/clang-13.1/include --with-hdf5-lib=/Users/baagaard/software/unix/hdf5-1.12/clang-13.1/lib/libhdf5.dylib --with-zlib=1
[0]PETSC ERROR: #1 PetscQuadratureGetData() at /Users/baagaard/software/unix/petsc-dev/src/dm/dt/interface/dt.c:255
[0]PETSC ERROR: #2 PetscFEGetCellTabulation() at /Users/baagaard/software/unix/petsc-dev/src/dm/dt/fe/interface/fe.c:791
[0]PETSC ERROR: #3 PetscDSSetUp() at /Users/baagaard/software/unix/petsc-dev/src/dm/dt/interface/dtds.c:409
[0]PETSC ERROR: #4 PetscDSGetComponentOffsets() at /Users/baagaard/software/unix/petsc-dev/src/dm/dt/interface/dtds.c:2961
[0]PETSC ERROR: #5 DMProjectPoint_Field_Private() at /Users/baagaard/software/unix/petsc-dev/src/dm/impls/plex/plexproject.c:201
[0]PETSC ERROR: #6 DMProjectPoint_Private() at /Users/baagaard/software/unix/petsc-dev/src/dm/impls/plex/plexproject.c:398
[0]PETSC ERROR: #7 DMProjectLocal_Generic_Plex() at /Users/baagaard/software/unix/petsc-dev/src/dm/impls/plex/plexproject.c:850
[0]PETSC ERROR: #8 DMProjectFieldLabelLocal_Plex() at /Users/baagaard/software/unix/petsc-dev/src/dm/impls/plex/plexproject.c:933
[0]PETSC ERROR: #9 DMProjectFieldLabelLocal() at /Users/baagaard/software/unix/petsc-dev/src/dm/interface/dm.c:8191
[0]PETSC ERROR: #10 DMProjectFieldLabel() at /Users/baagaard/software/unix/petsc-dev/src/dm/interface/dm.c:8267
[0]PETSC ERROR: #11 void pylith::meshio::OutputSubfield::projectWithLabel(const PetscVec &)() at /Users/baagaard/src/cig/pylith/libsrc/pylith/meshio/OutputSubfield.cc:221

```

diff --git a/libsrc/pylith/topology/Field.cc b/libsrc/pylith/topology/Field.cc
index 6786edab0..9abd65d2f 100644
--- a/libsrc/pylith/topology/Field.cc
+++ b/libsrc/pylith/topology/Field.cc
@@ -651,9 +651,14 @@ pylith::topology::Field::createOutputVector(void) {
 
     PetscDM dmOutput = NULL;
     err = DMGetOutputDM(_mesh->getDM(), &dmOutput);PYLITH_CHECK_ERROR(err);
+
     PetscDS dsOutput = NULL;
-    err = DMGetDS(dmOutput, &dsOutput);PYLITH_CHECK_ERROR(err);
-    err = PetscDSSetUp(dsOutput);PYLITH_CHECK_ERROR(err);
+    PetscInt numRegions = 0;
+    err = DMGetNumDS(dmOutput, &numRegions);PYLITH_CHECK_ERROR(err);
+    for (PetscInt i = 0; i < numRegions; ++i) {
+        err = DMGetRegionNumDS(dmOutput, i, NULL, NULL, &dsOutput, NULL);PYLITH_CHECK_ERROR(err);
+        err = PetscDSSetUp(dsOutput);PYLITH_CHECK_ERROR(err);
+    } // for
 
     err = DMCreateGlobalVector(dmOutput, &_outputVec);PYLITH_CHECK_ERROR(err);assert(_outputVec);
     err = PetscObjectSetName((PetscObject) _outputVec, getLabel());PYLITH_CHECK_ERROR(err);
diff --git a/tests/fullscale/linearelasticity/faults-2d/TestShearNoSlip.py b/tests/fullscale/linearelasticity/faults-2d/TestShearNoSlip.py
index 1a4e3af03..07c6e5e14 100644
--- a/tests/fullscale/linearelasticity/faults-2d/TestShearNoSlip.py
+++ b/tests/fullscale/linearelasticity/faults-2d/TestShearNoSlip.py
@@ -68,6 +68,16 @@ class TestCase(FullTestCase):
                 vertex_fields=["displacement"],
                 defaults=defaults,
             ),
+            Check(
+                mesh_entities=["fault_xmid"],
+                vertex_fields=["slip", "lagrange_multiplier_fault"],
+                defaults=defaults,
+            ),
+            Check(
+                mesh_entities=["fault_xneg"],
+                vertex_fields=["slip", "lagrange_multiplier_fault"],
+                defaults=defaults,
+            ),
         ]
 
     def run_pylith(self, testName, args):
diff --git a/tests/fullscale/linearelasticity/faults-2d/TestThreeBlocks.py b/tests/fullscale/linearelasticity/faults-2d/TestThreeBlocks.py
index c65956387..aa6cbb429 100644
--- a/tests/fullscale/linearelasticity/faults-2d/TestThreeBlocks.py
+++ b/tests/fullscale/linearelasticity/faults-2d/TestThreeBlocks.py
@@ -67,6 +67,16 @@ class TestCase(FullTestCase):
                 vertex_fields=["displacement"],
                 defaults=defaults,
             ),
+            Check(
+                mesh_entities=["fault_xmid"],
+                vertex_fields=["slip", "lagrange_multiplier_fault"],
+                defaults=defaults,
+            ),
+            Check(
+                mesh_entities=["fault_xneg"],
+                vertex_fields=["slip", "lagrange_multiplier_fault"],
+                defaults=defaults,
+            ),
         ]
 
     def run_pylith(self, testName, args):
diff --git a/tests/fullscale/linearelasticity/faults-2d/meshes.py b/tests/fullscale/linearelasticity/faults-2d/meshes.py
index c3b70c823..f37d3ca75 100644
--- a/tests/fullscale/linearelasticity/faults-2d/meshes.py
+++ b/tests/fullscale/linearelasticity/faults-2d/meshes.py
@@ -31,8 +31,8 @@ class TriGmsh(object):
         "mat_xposyneg": MeshEntity(ncells=44, ncorners=3, nvertices=31),
 
         # Faults
-        "fault": MeshEntity(ncells=8, ncorners=2, nvertices=9),
-        "fault2": MeshEntity(ncells=8, ncorners=2, nvertices=9),
+        "fault_xmid": MeshEntity(ncells=8, ncorners=2, nvertices=9),
+        "fault_xneg": MeshEntity(ncells=8, ncorners=2, nvertices=9),
 
         # Boundaries
         "bc_xneg": MeshEntity(ncells=8, ncorners=2, nvertices=9),
@@ -57,8 +57,8 @@ class QuadGmsh(object):
         "mat_xposyneg": MeshEntity(ncells=25, ncorners=3, nvertices=36),
 
         # Faults
-        "fault": MeshEntity(ncells=10, ncorners=2, nvertices=11),
-        "fault2": MeshEntity(ncells=10, ncorners=2, nvertices=11),
+        "fault_xmid": MeshEntity(ncells=10, ncorners=2, nvertices=11),
+        "fault_xneg": MeshEntity(ncells=10, ncorners=2, nvertices=11),
 
         # Boundaries
         "bc_xneg": MeshEntity(ncells=10, ncorners=2, nvertices=11),
@@ -83,8 +83,8 @@ class TriCubit(object):
         "mat_xposyneg": MeshEntity(ncells=32, ncorners=3, nvertices=25),
 
         # Faults
-        "fault": MeshEntity(ncells=8, ncorners=2, nvertices=9),
-        "fault2": MeshEntity(ncells=8, ncorners=2, nvertices=9),
+        "fault_xmid": MeshEntity(ncells=8, ncorners=2, nvertices=9),
+        "fault_xneg": MeshEntity(ncells=8, ncorners=2, nvertices=9),
 
         # Boundaries
         "bc_xneg": MeshEntity(ncells=8, ncorners=2, nvertices=9),
@@ -109,8 +109,8 @@ class QuadCubit(object):
         "mat_xposyneg": MeshEntity(ncells=16, ncorners=3, nvertices=25),
 
         # Faults
-        "fault": MeshEntity(ncells=8, ncorners=2, nvertices=9),
-        "fault2": MeshEntity(ncells=8, ncorners=2, nvertices=9),
+        "fault_xmid": MeshEntity(ncells=8, ncorners=2, nvertices=9),
+        "fault_xneg": MeshEntity(ncells=8, ncorners=2, nvertices=9),
 
         # Boundaries
         "bc_xneg": MeshEntity(ncells=8, ncorners=2, nvertices=9),
diff --git a/tests/fullscale/linearelasticity/faults-2d/pylithapp.cfg b/tests/fullscale/linearelasticity/faults-2d/pylithapp.cfg
index 7a45b1518..295967183 100644
--- a/tests/fullscale/linearelasticity/faults-2d/pylithapp.cfg
+++ b/tests/fullscale/linearelasticity/faults-2d/pylithapp.cfg
@@ -142,13 +142,13 @@ interfaces = [fault_xmid, fault_xneg]
 [pylithapp.problem.interfaces.fault_xmid]
 label = fault_xmid
 label_value = 20
-observers.observer.data_fields = [slip]
+observers.observer.data_fields = [slip, lagrange_multiplier_fault]
 
 
 [pylithapp.problem.interfaces.fault_xneg]
 label = fault_xneg
 label_value = 21
-observers.observer.data_fields = [slip]
+observers.observer.data_fields = [slip, lagrange_multiplier_fault]
 
 
 # ----------------------------------------------------------------------
diff --git a/tests/fullscale/linearelasticity/faults-2d/shearnoslip_soln.py b/tests/fullscale/linearelasticity/faults-2d/shearnoslip_soln.py
index 896419040..b393d07ab 100644
--- a/tests/fullscale/linearelasticity/faults-2d/shearnoslip_soln.py
+++ b/tests/fullscale/linearelasticity/faults-2d/shearnoslip_soln.py
@@ -83,7 +83,9 @@ class AnalyticalSoln(object):
                 "bc_xpos": self.displacement,
                 "bc_yneg": self.bc_traction,
                 "bc_ypos": self.bc_traction,
-            }
+            },
+            "slip": self.slip,
+            "lagrange_multiplier_fault": self.lagrange_multiplier_fault,
         }
         return
 
@@ -155,5 +157,20 @@ class AnalyticalSoln(object):
         traction[0,:, 1] = 0.0
         return traction
 
+    def slip(self, locs):
+        """Compute slip field.
+        """
+        (npts, dim) = locs.shape
+        slip = numpy.zeros((1, npts, 2), dtype=numpy.float64)
+        return slip
+
+    def lagrange_multiplier_fault(self, locs):
+        """Compute Lagrange multiplier on faults.
+        """
+        (npts, dim) = locs.shape
+        traction = numpy.zeros((1, npts, 2), dtype=numpy.float64)
+        traction[:,:,1] = sxy
+        return traction
+
 
 # End of file
diff --git a/tests/fullscale/linearelasticity/faults-2d/threeblocks_soln.py b/tests/fullscale/linearelasticity/faults-2d/threeblocks_soln.py
index 0322ce459..1160114a1 100644
--- a/tests/fullscale/linearelasticity/faults-2d/threeblocks_soln.py
+++ b/tests/fullscale/linearelasticity/faults-2d/threeblocks_soln.py
@@ -73,11 +73,21 @@ class AnalyticalSoln(object):
             "cauchy_strain": self.strain,
             "cauchy_stress": self.stress,
             "initial_amplitude": self.displacement,
+            "slip": {
+                "fault_xmid": self.slip_xmid,
+                "fault_xneg": self.slip_xneg,
+            },
+            "lagrange_multiplier_fault": self.lagrange_multiplier_fault,
         }
         return
 
     def getField(self, name, mesh_entity, pts):
-        return self.fields[name](pts)
+        if name == "slip":
+            field = self.fields[name][mesh_entity](pts)
+        else:
+            field = self.fields[name](pts)
+
+        return field
 
     def getMask(self, name, mesh_entity, pts):
         mask = None
@@ -137,5 +147,29 @@ class AnalyticalSoln(object):
         stress[0,:, 3] = sxy
         return stress
 
+    def slip_xmid(self, locs):
+        """Compute slip field on fault xmid.
+        """
+        (npts, dim) = locs.shape
+        slip = numpy.zeros((1, npts, 2), dtype=numpy.float64)
+        slip[0,:, 1] = -2.0
+        return slip
+
+    def slip_xneg(self, locs):
+        """Compute slip field on fault xneg.
+        """
+        (npts, dim) = locs.shape
+        slip = numpy.zeros((1, npts, 2), dtype=numpy.float64)
+        slip[0,:, 1] = +2.0
+        return slip
+
+    def lagrange_multiplier_fault(self, locs):
+        """Compute Lagrange multiplier on faults.
+        """
+        (npts, dim) = locs.shape
+        traction = numpy.zeros((1, npts, 2), dtype=numpy.float64)
+        return traction
+
+
 
 # End of file
diff --git a/tests/fullscale/linearelasticity/faults-3d/TestTwoBlocks.py b/tests/fullscale/linearelasticity/faults-3d/TestTwoBlocks.py
index 358904c6b..dc6f19287 100644
--- a/tests/fullscale/linearelasticity/faults-3d/TestTwoBlocks.py
+++ b/tests/fullscale/linearelasticity/faults-3d/TestTwoBlocks.py
@@ -67,6 +67,17 @@ class TestCase(FullTestCase):
                 vertex_fields=["displacement"],
                 defaults=defaults,
             ),
+            Check(
+                mesh_entities=["fault"],
+                vertex_fields=["slip"],
+                defaults=defaults,
+            ),
+            Check(
+                mesh_entities=["fault"],
+                vertex_fields=["lagrange_multiplier_fault"],
+                scale = 1.0e+6,
+                defaults=defaults,
+            ),
         ]
 
     def run_pylith(self, testName, args, nprocs=1):
diff --git a/tests/fullscale/linearelasticity/faults-3d/meshes.py b/tests/fullscale/linearelasticity/faults-3d/meshes.py
index acebac22f..cc24d760f 100644
--- a/tests/fullscale/linearelasticity/faults-3d/meshes.py
+++ b/tests/fullscale/linearelasticity/faults-3d/meshes.py
@@ -51,7 +51,7 @@ class HexGmsh(object):
         "mat_elastic": MeshEntity(ncells=150, ncorners=8, nvertices=252+36),
 
         # Faults
-        "fault": MeshEntity(ncells=10, ncorners=4, nvertices=11),
+        "fault": MeshEntity(ncells=25, ncorners=4, nvertices=36),
 
         # Boundaries
         "bc_xneg": MeshEntity(ncells=25, ncorners=4, nvertices=36),
diff --git a/tests/fullscale/linearelasticity/faults-3d/twoblocks.cfg b/tests/fullscale/linearelasticity/faults-3d/twoblocks.cfg
index 99fe268bf..a4cf512bc 100644
--- a/tests/fullscale/linearelasticity/faults-3d/twoblocks.cfg
+++ b/tests/fullscale/linearelasticity/faults-3d/twoblocks.cfg
@@ -34,7 +34,7 @@ interfaces = [fault]
 label = fault
 label_value = 20
 
-observers.observer.data_fields = [slip]
+observers.observer.data_fields = [slip, lagrange_multiplier_fault]
 
 [pylithapp.problem.interfaces.fault.eq_ruptures.rupture]
 db_auxiliary_field = spatialdata.spatialdb.UniformDB
diff --git a/tests/fullscale/linearelasticity/faults-3d/twoblocks_soln.py b/tests/fullscale/linearelasticity/faults-3d/twoblocks_soln.py
index c315424d8..7e942c45f 100644
--- a/tests/fullscale/linearelasticity/faults-3d/twoblocks_soln.py
+++ b/tests/fullscale/linearelasticity/faults-3d/twoblocks_soln.py
@@ -76,6 +76,8 @@ class AnalyticalSoln(object):
             "cauchy_strain": self.strain,
             "cauchy_stress": self.stress,
             "initial_amplitude": self.displacement,
+            "slip": self.slip,
+            "lagrange_multiplier_fault": self.lagrange_multiplier_fault,
         }
         return
 
@@ -149,5 +151,21 @@ class AnalyticalSoln(object):
         stress[0,:, 3] = sxy
         return stress
 
+    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
+        slip[:,:,2] = +1.0
+        return slip
+
+    def lagrange_multiplier_fault(self, locs):
+        """Compute change in fault traction field at locations.
+        """
+        (npts, dim) = locs.shape
+        traction_change = numpy.zeros((1, npts, 3), dtype=numpy.float64)
+        return traction_change
+
 
 # End of file
['tests/fullscale/linearelasticity/faults-2d/TestShearNoSlip.py', 'tests/fullscale/linearelasticity/faults-2d/threeblocks_soln.py', 'tests/fullscale/linearelasticity/faults-3d/TestTwoBlocks.py', 'tests/fullscale/linearelasticity/faults-2d/shearnoslip_soln.py', 'tests/fullscale/linearelasticity/faults-3d/twoblocks_soln.py', 'libsrc/pylith/topology/Field.cc', 'tests/fullscale/linearelasticity/faults-2d/meshes.py', 'tests/fullscale/linearelasticity/faults-2d/TestThreeBlocks.py', 'tests/fullscale/linearelasticity/faults-2d/pylithapp.cfg', 'tests/fullscale/linearelasticity/faults-3d/twoblocks.cfg', 'tests/fullscale/linearelasticity/faults-3d/meshes.py']