[cig-commits] r7117 - in short/3D/PyLith/trunk: . libsrc libsrc/faults libsrc/feassemble modulesrc/feassemble pylith pylith/feassemble/quadrature unittests/libtests/faults unittests/libtests/faults/data unittests/libtests/feassemble unittests/pytests/feassemble

[brad at geodynamics.org]

Author: brad
Date: 2007-06-08 18:30:56 -0700 (Fri, 08 Jun 2007)
New Revision: 7117

Added:
   short/3D/PyLith/trunk/libsrc/feassemble/Quadrature0D.cc
   short/3D/PyLith/trunk/libsrc/feassemble/Quadrature0D.hh
   short/3D/PyLith/trunk/libsrc/feassemble/Quadrature0D.icc
   short/3D/PyLith/trunk/pylith/feassemble/quadrature/Quadrature0D.py
   short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKinTri3.cc
   short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKinTri3.hh
   short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinData.cc
   short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinData.hh
   short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataLine2.cc
   short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataLine2.hh
   short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTri3.cc
   short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTri3.hh
   short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature0D.cc
   short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature0D.hh
Modified:
   short/3D/PyLith/trunk/TODO
   short/3D/PyLith/trunk/libsrc/Makefile.am
   short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
   short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am
   short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc
   short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe.src
   short/3D/PyLith/trunk/pylith/Makefile.am
   short/3D/PyLith/trunk/pylith/feassemble/quadrature/__init__.py
   short/3D/PyLith/trunk/unittests/libtests/faults/Makefile.am
   short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
   short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.hh
   short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKinLine2.cc
   short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKinLine2.hh
   short/3D/PyLith/trunk/unittests/libtests/faults/test_faults.cc
   short/3D/PyLith/trunk/unittests/libtests/feassemble/Makefile.am
   short/3D/PyLith/trunk/unittests/pytests/feassemble/TestQuadrature.py
Log:
Worked on C++ implementation and unit tests for kinematic faults with cohesive elements. Fault mesh created in generating cohesive cells doesn't seem useful.

Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO	2007-06-08 23:11:50 UTC (rev 7116)
+++ short/3D/PyLith/trunk/TODO	2007-06-09 01:30:56 UTC (rev 7117)
@@ -2,6 +2,19 @@
 MAIN PRIORITIES (Brad)
 ======================================================================
 
+ISSUES: 
+
+  * What is the fault mesh supposed to be??????
+
+  * BruneSlipFn assumes fault mesh is associated with constraint
+    vertices (THIS IS NOT THE CASE).
+
+    The SlipTimeFn stuff doesn't care about topology, only
+    points. Need to eliminate passing both the constraint points and
+    the fault mesh; use one or the other but not both. Perhaps the
+    fault mesh should be composed of the constraint vertices; this
+    would eliminate the need to store the set of constraint vertices.
+
 1. Simple tests with analytical solutions
 
    a. 2-D

Modified: short/3D/PyLith/trunk/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/Makefile.am	2007-06-08 23:11:50 UTC (rev 7116)
+++ short/3D/PyLith/trunk/libsrc/Makefile.am	2007-06-09 01:30:56 UTC (rev 7117)
@@ -51,6 +51,7 @@
 	feassemble/GeometryHex3D.cc \
 	feassemble/Integrator.cc \
 	feassemble/Quadrature.cc \
+	feassemble/Quadrature0D.cc \
 	feassemble/Quadrature1D.cc \
 	feassemble/Quadrature1Din2D.cc \
 	feassemble/Quadrature1Din3D.cc \

Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc	2007-06-08 23:11:50 UTC (rev 7116)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc	2007-06-09 01:30:56 UTC (rev 7117)
@@ -58,78 +58,115 @@
 { // initialize
   assert(0 != _quadrature);
   assert(0 != _eqsrc);
-  assert(0 != _faultMesh);
-  assert(!_faultMesh->isNull());
   
   if (3 != upDir.size())
     throw std::runtime_error("Up direction for fault orientation must be "
 			     "a vector with 3 components.");
 
-  /* First use fault mesh to get orientation of vertices in fault mesh
-   * We will transfer the orientation from the fault mesh vertices to
-   * the Lagrange constraint vertices.
+  /* First find vertices associated with Lagrange multiplier
+   * constraints in cohesive cells and compute the orientation of the
+   * fault at these locations.
    */
 
-  // Allocate section for orientation of vertices in fault mesh
-  ALE::Obj<real_section_type> orientation = 
-    new real_section_type((*_faultMesh)->comm(), (*_faultMesh)->debug());
-  assert(!orientation.isNull());
-  const int cellDim = (*_faultMesh)->getDimension();
+  const ALE::Obj<sieve_type>& sieve = mesh->getSieve();
+  assert(!sieve.isNull());
+  const ALE::Obj<ALE::Mesh::label_sequence>& cellsCohesive = 
+    mesh->getLabelStratum("material-id", id());
+  assert(!cellsCohesive.isNull());
+  const Mesh::label_sequence::iterator cellsCohesiveBegin =
+    cellsCohesive->begin();
+  const Mesh::label_sequence::iterator cellsCohesiveEnd =
+    cellsCohesive->end();
+
+  // Create set of vertices associated with Lagrange multiplier constraints
+  _constraintVert.clear();
+  for (Mesh::label_sequence::iterator c_iter=cellsCohesiveBegin;
+       c_iter != cellsCohesiveEnd;
+       ++c_iter) {
+    // Vertices for each cohesive cell are in three groups of N.
+    // 0 to N-1: vertices on negative side of the fault
+    // N-1 to 2N-1: vertices on positive side of the fault
+    // 2N to 3N-1: vertices associated with constraint forces
+    const ALE::Obj<sieve_type::traits::coneSequence>& cone = 
+      sieve->cone(*c_iter);
+    assert(!cone.isNull());
+    const sieve_type::traits::coneSequence::iterator vBegin = cone->begin();
+    const sieve_type::traits::coneSequence::iterator vEnd = cone->end();
+    const int coneSize = cone->size();
+    assert(coneSize % 3 == 0);
+    sieve_type::traits::coneSequence::iterator v_iter = vBegin;
+    // Skip over non-constraint vertices
+    for (int i=0, numSkip=2*coneSize/3; i < numSkip; ++i)
+      ++v_iter;
+    // Add constraint vertices to set
+    for(int i=0, numConstraintVert=coneSize/3; 
+	i < numConstraintVert; 
+	++i, ++v_iter)
+      _constraintVert.insert(*v_iter);
+  } // for
+
+  // Create orientation section for constraint vertices
+  const int cohesiveDim = mesh->getDimension()-1;
   const int spaceDim = cs->spaceDim();
-  const int orientationSize = (cellDim > 0) ? cellDim*spaceDim : 1;
-  orientation->setFiberDimension((*_faultMesh)->depthStratum(0), 
-				 orientationSize);
-  (*_faultMesh)->allocate(orientation);
-  
+  const int orientationSize = (cohesiveDim > 0) ? cohesiveDim*spaceDim : 1;
+  _orientation = new real_section_type(mesh->comm(), mesh->debug());
+  assert(!_orientation.isNull());
+  const std::set<Mesh::point_type>::const_iterator vertCohesiveBegin = 
+    _constraintVert.begin();
+  const std::set<Mesh::point_type>::const_iterator vertCohesiveEnd = 
+    _constraintVert.end();
+  for (std::set<Mesh::point_type>::const_iterator v_iter=vertCohesiveBegin;
+       v_iter != vertCohesiveEnd;
+       ++v_iter)
+    _orientation->setFiberDimension(*v_iter, orientationSize);
+  mesh->allocate(_orientation);
+
+  // Compute orientation of fault at constraint vertices
+
   // Get section containing coordinates of vertices
   const ALE::Obj<real_section_type>& coordinates = 
     mesh->getRealSection("coordinates");
   assert(!coordinates.isNull());
 
   // Set orientation function
-  assert(cellDim == _quadrature->cellDim());
+  assert(cohesiveDim == _quadrature->cellDim());
   assert(spaceDim == _quadrature->spaceDim());
   orient_fn_type orientFn;
-  switch (cellDim)
+  switch (cohesiveDim)
     { // switch
-    case 1 :
+    case 0 :
       orientFn = _orient1D;
       break;
-    case 2 :
+    case 1 :
       orientFn = _orient2D;
       break;
-    case 3 :
+    case 2 :
       orientFn = _orient3D;
       break;
     default :
       assert(0);
     } // switch
 
-  // Loop over cells in fault mesh, computing orientation at each vertex
-  const ALE::Obj<sieve_type>& sieve = (*_faultMesh)->getSieve();
-  assert(!sieve.isNull());
+  // Loop over cohesive cells, computing orientation at each vertex
 
   const int numBasis = _quadrature->numBasis();
   const feassemble::CellGeometry& cellGeometry = _quadrature->refGeometry();
   const double_array& verticesRef = _quadrature->vertices();
-  const int jacobianSize = spaceDim * cellDim;
+  const int jacobianSize = spaceDim * cohesiveDim;
   double_array jacobian(jacobianSize);
   double jacobianDet = 0;
   double_array vertexOrientation(orientationSize);
-  double_array cellVertices(numBasis*spaceDim);
+  double_array cohesiveVertices(3*numBasis*spaceDim);
+  double_array faceVertices(numBasis*spaceDim);
 
-  const ALE::Obj<Mesh::label_sequence>& cellsFault = 
-    (*_faultMesh)->heightStratum(0);
-  assert(!cellsFault.isNull());
-  const Mesh::label_sequence::iterator cellsFaultBegin = cellsFault->begin();
-  const Mesh::label_sequence::iterator cellsFaultEnd = cellsFault->end();
-  for (Mesh::label_sequence::iterator c_iter=cellsFaultBegin;
-       c_iter != cellsFaultEnd;
+  for (Mesh::label_sequence::iterator c_iter=cellsCohesiveBegin;
+       c_iter != cellsCohesiveEnd;
        ++c_iter) {
     mesh->restrict(coordinates, *c_iter, 
-		   &cellVertices[0], cellVertices.size());
+		   &cohesiveVertices[0], cohesiveVertices.size());
+    for (int i=0, offset=2*numBasis; i < numBasis; ++i)
+      faceVertices[i] = cohesiveVertices[offset+i];
 
-    // Compute orientation at each vertex
     int iBasis = 0;
     const ALE::Obj<sieve_type::traits::coneSequence>& cone = 
       sieve->cone(*c_iter);
@@ -140,102 +177,42 @@
 	v_iter != vEnd;
 	++v_iter, ++iBasis) {
       // Compute Jacobian and determinant of Jacobian at vertex
-      double_array vertex(&verticesRef[iBasis*cellDim], cellDim);
-      cellGeometry.jacobian(&jacobian, &jacobianDet, cellVertices, vertex);
+      double_array vertex(&verticesRef[iBasis*cohesiveDim], cohesiveDim);
+      cellGeometry.jacobian(&jacobian, &jacobianDet, faceVertices, vertex);
 
       // Compute orientation
       orientFn(&vertexOrientation, jacobian, jacobianDet, upDir);
       
       // Update orientation
-      orientation->updatePoint(*v_iter, &vertexOrientation[0]);
+      _orientation->updatePoint(*v_iter, &vertexOrientation[0]);
     } // for
   } // for
 
   // Assemble orientation information
-  //orientation->complete();
+  // FIX THIS
+  //const ALE::Obj<Mesh>& bundle = orientation.b;
+  //ALE::Distribution<Mesh>::completeSection(bundle, orientation);
 
   // Loop over vertices, make orientation information unit magnitude
-  const ALE::Obj<Mesh::label_sequence>& vertices = 
-    (*_faultMesh)->depthStratum(0);
-  const Mesh::label_sequence::iterator vertFaultBegin = vertices->begin();
-  const Mesh::label_sequence::iterator vertFaultEnd = vertices->end();
   double_array vertexDir(orientationSize);
-  for (Mesh::label_sequence::iterator v_iter=vertFaultBegin;
-       v_iter != vertFaultEnd;
+  for (std::set<Mesh::point_type>::const_iterator v_iter=vertCohesiveBegin;
+       v_iter != vertCohesiveEnd;
        ++v_iter) {
     const real_section_type::value_type* vertexOrient = 
-      orientation->restrictPoint(*v_iter);
+      _orientation->restrictPoint(*v_iter);
     
-    assert(cellDim*spaceDim == orientationSize);
-    for (int iDim=0, index=0; iDim < cellDim; ++iDim, index+=cellDim) {
+    assert(cohesiveDim*spaceDim == orientationSize);
+    for (int iDim=0, index=0; iDim < cohesiveDim; ++iDim, index+=cohesiveDim) {
       double mag = 0;
       for (int jDim=0; jDim < spaceDim; ++jDim)
-	mag *= vertexOrient[index*cellDim+jDim];
-      for (int jDim=0; jDim < cellDim; ++jDim)
-	vertexDir[index*cellDim+jDim] = vertexOrient[index*cellDim+jDim] / mag;
+	mag *= vertexOrient[index*cohesiveDim+jDim];
+      for (int jDim=0; jDim < cohesiveDim; ++jDim)
+	vertexDir[index*cohesiveDim+jDim] = 
+	  vertexOrient[index*cohesiveDim+jDim] / mag;
     } // for
-    orientation->updatePoint(*v_iter, &vertexDir[0]);
+    _orientation->updatePoint(*v_iter, &vertexDir[0]);
   } // for
 
-  _constraintVert.clear();
-  // Create set of vertices associated with Lagrange multiplier constraints
-  const ALE::Obj<ALE::Mesh::label_sequence>& cellsCohesive = 
-    mesh->getLabelStratum("material-id", id());
-  assert(!cellsCohesive.isNull());
-  const Mesh::label_sequence::iterator cellsCohesiveBegin =
-    cellsCohesive->begin();
-  const Mesh::label_sequence::iterator cellsCohesiveEnd =
-    cellsCohesive->end();
-  for (Mesh::label_sequence::iterator c_iter=cellsCohesiveBegin;
-       c_iter != cellsCohesiveEnd;
-       ++c_iter) {
-    // Vertices for each cohesive cell are in groups of N.
-    // 0 to N-1: vertices on negative side of the fault
-    // N-1 to 2N-1: vertices on positive side of the fault
-    // 2N to 3N-1: vertices associated with constraint forces
-    const ALE::Obj<sieve_type::traits::coneSequence>& cone = 
-      sieve->cone(*c_iter);
-    assert(!cone.isNull());
-    const sieve_type::traits::coneSequence::iterator vBegin = cone->begin();
-    const sieve_type::traits::coneSequence::iterator vEnd = cone->end();
-    const int coneSize = cone->size();
-    assert(coneSize % 3 == 0);
-    sieve_type::traits::coneSequence::iterator v_iter = vBegin;
-    // Skip over non-constraint vertices
-    for (int i=0, numSkip=2*coneSize/3; i < numSkip; ++i)
-      ++v_iter;
-    // Add constraint vertices to set
-    for(int i=0, numConstraintVert=coneSize/3; 
-	i < numConstraintVert; 
-	++i, ++v_iter)
-      _constraintVert.insert(*v_iter);
-  } // for
-
-  // Create orientation section over vertices associated with Lagrange
-  // multiplier constraints
-  _orientation = new real_section_type(mesh->comm(), mesh->debug());
-  assert(!_orientation.isNull());
-  const std::set<Mesh::point_type>::const_iterator vertCohesiveBegin = 
-    _constraintVert.begin();
-  const std::set<Mesh::point_type>::const_iterator vertCohesiveEnd = 
-    _constraintVert.end();
-  for (std::set<Mesh::point_type>::const_iterator v_iter=vertCohesiveBegin;
-       v_iter != vertCohesiveEnd;
-       ++v_iter)
-    _orientation->setFiberDimension(*v_iter, orientationSize);
-  mesh->allocate(_orientation);
-  
-  // Transfer orientation information from fault vertices to vertices
-  // associated with Lagrange multiplier constraints
-  Mesh::label_sequence::iterator vFault_iter=vertFaultBegin;
-  for (std::set<Mesh::point_type>::const_iterator vConstraint_iter=vertCohesiveBegin;
-       vConstraint_iter != vertCohesiveEnd;
-       ++vConstraint_iter, ++vFault_iter) {
-    const real_section_type::value_type* vertexOrient = 
-      orientation->restrictPoint(*vFault_iter);
-    _orientation->updatePoint(*vConstraint_iter, vertexOrient);
-  } // for
-  
   _eqsrc->initialize(mesh, *_faultMesh, _constraintVert, cs);
 } // initialize
 
@@ -250,8 +227,6 @@
   assert(!residual.isNull());
   assert(0 != fields);
   assert(!mesh.isNull());
-  assert(0 != _faultMesh);
-  assert(!_faultMesh->isNull());
 
   // Subtract constraint forces (which are disp at the constraint
   // DOF) to residual; contributions are at DOF of normal vertices (i and j)
@@ -310,8 +285,6 @@
   assert(0 != mat);
   assert(0 != fields);
   assert(!mesh.isNull());
-  assert(0 != _faultMesh);
-  assert(!_faultMesh->isNull());
 
   // Add constraint information to Jacobian matrix; these are the
   // direction cosines. Entries are associated with vertices ik, jk,
@@ -330,9 +303,9 @@
   const ALE::Obj<real_section_type>& disp = fields->getHistoryItem(1);
   assert(!disp.isNull());  
 
-  const int cellDim = _quadrature->cellDim();
+  const int cohesiveDim = _quadrature->cellDim();
   const int spaceDim = _quadrature->spaceDim();
-  const int orientationSize = cellDim*spaceDim;
+  const int orientationSize = cohesiveDim*spaceDim;
 
   // Allocate matrix for cell values (if necessary)
   _initCellMatrix();
@@ -395,6 +368,7 @@
     if (err)
       throw std::runtime_error("Update to PETSc Mat failed.");
   } // for
+  _needNewJacobian = false;
 } // integrateJacobian
   
 // ----------------------------------------------------------------------

Modified: short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am	2007-06-08 23:11:50 UTC (rev 7116)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am	2007-06-09 01:30:56 UTC (rev 7117)
@@ -35,6 +35,8 @@
 	GeometryHex3D.hh \
 	Quadrature.hh \
 	Quadrature.icc \
+	Quadrature0D.hh \
+	Quadrature0D.icc \
 	Quadrature1D.hh \
 	Quadrature1D.icc \
 	Quadrature1Din2D.hh \

Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc	2007-06-08 23:11:50 UTC (rev 7116)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc	2007-06-09 01:30:56 UTC (rev 7117)
@@ -82,7 +82,7 @@
       0 == basisDeriv ||
       0 == quadPtsRef ||
       0 == quadWts ||
-      cellDim < 1 || cellDim > 3 ||
+      cellDim < 0 || cellDim > 3 ||
       numBasis < 1 ||
       numQuadPts < 1 ||
       spaceDim < 1 || spaceDim > 3) {
@@ -103,48 +103,107 @@
     throw std::runtime_error(msg.str());
   } // if
 
-  int size = numBasis * cellDim;
-  _vertices.resize(size);
-  for (int i=0; i < size; ++i)
-    _vertices[i] = vertices[i];
+  if (cellDim > 0) {
+    int size = numBasis * cellDim;
+    _vertices.resize(size);
+    for (int i=0; i < size; ++i)
+      _vertices[i] = vertices[i];
 
-  size = numBasis * numQuadPts;
-  _basis.resize(size);
-  for (int i=0; i < size; ++i)
-    _basis[i] = basis[i];
+    size = numBasis * numQuadPts;
+    _basis.resize(size);
+    for (int i=0; i < size; ++i)
+      _basis[i] = basis[i];
 
-  size = numBasis * numQuadPts * cellDim;
-  _basisDeriv.resize(size);
-  for (int i=0; i < size; ++i)
-    _basisDeriv[i] = basisDeriv[i];
+    size = numBasis * numQuadPts * cellDim;
+    _basisDeriv.resize(size);
+    for (int i=0; i < size; ++i)
+      _basisDeriv[i] = basisDeriv[i];
 
-  size = numQuadPts * cellDim;
-  _quadPtsRef.resize(size);
-  for (int i=0; i < size; ++i)
-    _quadPtsRef[i] = quadPtsRef[i];
+    size = numQuadPts * cellDim;
+    _quadPtsRef.resize(size);
+    for (int i=0; i < size; ++i)
+      _quadPtsRef[i] = quadPtsRef[i];
 
-  size = numQuadPts;
-  _quadWts.resize(size);
-  for (int i=0; i < size; ++i)
-    _quadWts[i] = quadWts[i];
+    size = numQuadPts;
+    _quadWts.resize(size);
+    for (int i=0; i < size; ++i)
+      _quadWts[i] = quadWts[i];
 
-  _cellDim = cellDim;
-  _numBasis = numBasis;
-  _numQuadPts = numQuadPts;
-  _spaceDim = spaceDim;
+    _cellDim = cellDim;
+    _numBasis = numBasis;
+    _numQuadPts = numQuadPts;
+    _spaceDim = spaceDim;
 
-  // Allocate for Jacobian and its inverse
-  size = numQuadPts * cellDim * spaceDim;
-  _jacobian.resize(size);
-  _jacobianInv.resize(size);
+    // Allocate for Jacobian and its inverse
+    size = numQuadPts * cellDim * spaceDim;
+    _jacobian.resize(size);
+    _jacobianInv.resize(size);
 
-  // Allocate for Jacobian determinant
-  size = numQuadPts;
-  _jacobianDet.resize(size);
+    // Allocate for Jacobian determinant
+    size = numQuadPts;
+    _jacobianDet.resize(size);
 
-  // Allocate for quad pts
-  size = numQuadPts*spaceDim;
-  _quadPts.resize(size);
+    // Allocate for quad pts
+    size = numQuadPts*spaceDim;
+    _quadPts.resize(size);
+  } else {
+    if (1 != numBasis ||
+	1 != numQuadPts ||
+	1 != spaceDim) {
+      std::ostringstream msg;
+      msg << "0-D quadrature only works in 1-D and is limited to 1 basis "
+	  << "function and 1 quadrature point.\n"
+	  << "Values:\n"
+	  << "  cell dimension: " << cellDim << "\n"
+	  << "  spatial dimension: " << spaceDim << "\n"
+	  << "  # vertices per cell: " << numBasis << "\n"
+	  << "  # quadrature points: " << numQuadPts << "\n";
+      throw std::runtime_error(msg.str());
+    } // if
+
+    int size = 1;
+    _vertices.resize(size);
+    for (int i=0; i < size; ++i)
+      _vertices[i] = vertices[i];
+
+    size = 1;
+    _basis.resize(size);
+    for (int i=0; i < size; ++i)
+      _basis[i] = basis[i];
+
+    size = 1;
+    _basisDeriv.resize(size);
+    for (int i=0; i < size; ++i)
+      _basisDeriv[i] = basisDeriv[i];
+
+    size = 1;
+    _quadPtsRef.resize(size);
+    for (int i=0; i < size; ++i)
+      _quadPtsRef[i] = quadPtsRef[i];
+
+    size = 1;
+    _quadWts.resize(size);
+    for (int i=0; i < size; ++i)
+      _quadWts[i] = quadWts[i];
+
+    _cellDim = cellDim;
+    _numBasis = numBasis;
+    _numQuadPts = numQuadPts;
+    _spaceDim = spaceDim;
+
+    // Allocate for Jacobian and its inverse
+    size = 1;
+    _jacobian.resize(size);
+    _jacobianInv.resize(size);
+
+    // Allocate for Jacobian determinant
+    size = 1;
+    _jacobianDet.resize(size);
+
+    // Allocate for quad pts
+    size = spaceDim;
+    _quadPts.resize(size);
+  } // else
 } // initialize
 
 // ----------------------------------------------------------------------

Added: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature0D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature0D.cc	2007-06-08 23:11:50 UTC (rev 7116)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature0D.cc	2007-06-09 01:30:56 UTC (rev 7117)
@@ -0,0 +1,68 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+//                           Brad T. Aagaard
+//                        U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#include <portinfo>
+
+#include "Quadrature0D.hh" // implementation of class methods
+
+#include "pylith/utils/array.hh" // USES double_array
+
+#include <assert.h> // USES assert()
+
+// ----------------------------------------------------------------------
+// Constructor
+pylith::feassemble::Quadrature0D::Quadrature0D(void)
+{ // constructor
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor
+pylith::feassemble::Quadrature0D::~Quadrature0D(void)
+{ // destructor
+} // destructor
+  
+// ----------------------------------------------------------------------
+// Copy constructor.
+pylith::feassemble::Quadrature0D::Quadrature0D(const Quadrature0D& q) :
+  Quadrature(q)
+{ // copy constructor
+} // copy constructor
+
+// ----------------------------------------------------------------------
+// Compute geometric quantities for a cell at quadrature points.
+void
+pylith::feassemble::Quadrature0D::computeGeometry(
+			      const ALE::Obj<Mesh>& mesh,
+			      const ALE::Obj<real_section_type>& coordinates,
+			      const Mesh::point_type& cell)
+{ // computeGeometry
+  assert(0 == _cellDim);
+  assert(1 == _numQuadPts);
+  assert(1 == _numBasis);
+
+  _resetGeometry();
+  
+  // Get coordinates of cell's vertices
+  const real_section_type::value_type* vertCoords = 
+    mesh->restrict(coordinates, cell);
+  assert(1 == coordinates->getFiberDimension(*mesh->depthStratum(0)->begin()));
+
+  for (int i=0; i < _spaceDim; ++i)
+    _quadPts[i] = vertCoords[i];
+
+  _jacobian[0] = 1.0;
+  _jacobianDet[0] = 1.0;
+  _jacobianInv[0] = 1.0;
+} // computeGeometry
+
+
+// End of file 

Added: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature0D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature0D.hh	2007-06-08 23:11:50 UTC (rev 7116)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature0D.hh	2007-06-09 01:30:56 UTC (rev 7117)
@@ -0,0 +1,80 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+//                           Brad T. Aagaard
+//                        U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+/**
+ * @file pylith/feassemble/Quadrature0D.hh
+ *
+ * @brief Quadrature for 0-D finite-elements.
+ *
+ * Need Quadrature in 0-D for integration of boundary condition for
+ * 1-D meshes.
+ */
+
+#if !defined(pylith_feassemble_quadrature0d_hh)
+#define pylith_feassemble_quadrature0d_hh
+
+#include "Quadrature.hh"
+
+namespace pylith {
+  namespace feassemble {
+    class Quadrature0D;
+    class TestQuadrature0D;
+  } // feassemble
+} // pylith
+
+class pylith::feassemble::Quadrature0D : public Quadrature
+{ // Quadrature0D
+  friend class TestQuadrature0D; // unit testing
+
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public :
+
+  /// Constructor
+  Quadrature0D(void);
+
+  /// Destructor
+  ~Quadrature0D(void);
+
+  /// Create a copy of this object.
+  Quadrature* clone(void) const;
+
+  /** Compute geometric quantities for a cell at quadrature points.
+   *
+   * @param mesh Finite-element mesh
+   * @param coordinates Section containing vertex coordinates
+   * @param cell Finite-element cell
+   */
+  void computeGeometry(const ALE::Obj<Mesh>& mesh,
+		       const ALE::Obj<real_section_type>& coordinates,
+		       const Mesh::point_type& cell);
+
+// PROTECTED METHODS ////////////////////////////////////////////////////
+protected :
+
+  /** Copy constructor.
+   *
+   * @param q Quadrature to copy
+   */
+  Quadrature0D(const Quadrature0D& q);
+
+// PRIVATE METHODS //////////////////////////////////////////////////////
+private :
+
+  const Quadrature0D& operator=(const Quadrature0D&); ///< Not implemented
+
+}; // Quadrature0D
+
+#include "Quadrature0D.icc" // inline methods
+
+#endif // pylith_feassemble_quadrature0d_hh
+
+// End of file 

Added: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature0D.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature0D.icc	2007-06-08 23:11:50 UTC (rev 7116)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature0D.icc	2007-06-09 01:30:56 UTC (rev 7117)
@@ -0,0 +1,26 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+//                           Brad T. Aagaard
+//                        U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#if !defined(pylith_feassemble_quadrature0d_hh)
+#error "Quadrature0D.icc must be included only from Quadrature0D.hh"
+#else
+
+// Create a copy of this object.
+inline
+pylith::feassemble::Quadrature*
+pylith::feassemble::Quadrature0D::clone(void) const {
+  return new Quadrature0D(*this);
+} // clone
+
+#endif
+
+// End of file

Modified: short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe.src
===================================================================
--- short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe.src	2007-06-08 23:11:50 UTC (rev 7116)
+++ short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe.src	2007-06-09 01:30:56 UTC (rev 7117)
@@ -26,6 +26,7 @@
 #include "pylith/feassemble/GeometryHex3D.hh"
 
 #include "pylith/feassemble/Quadrature.hh"
+#include "pylith/feassemble/Quadrature0D.hh"
 #include "pylith/feassemble/Quadrature1D.hh"
 #include "pylith/feassemble/Quadrature1Din2D.hh"
 #include "pylith/feassemble/Quadrature1Din3D.hh"
@@ -668,6 +669,33 @@
 
 
 # ----------------------------------------------------------------------
+cdef class Quadrature0D(Quadrature):
+
+  def __init__(self):
+    """Constructor."""
+    # create shim for constructor
+    #embed{ void* Quadrature0D_constructor()
+    void* result = 0;
+    try {
+      result = (void*)(new pylith::feassemble::Quadrature0D);
+      assert(0 != result);
+    } catch (const std::exception& err) {
+      PyErr_SetString(PyExc_RuntimeError,
+                      const_cast<char*>(err.what()));
+    } catch (...) {
+      PyErr_SetString(PyExc_RuntimeError,
+                      "Caught unknown C++ exception.");
+    } // try/catch
+    return result;
+    #}embed
+
+    Quadrature.__init__(self)
+    self.thisptr = Quadrature0D_constructor()
+    self.handle = self._createHandle()
+    return
+
+
+# ----------------------------------------------------------------------
 cdef class Quadrature1D(Quadrature):
 
   def __init__(self):

Modified: short/3D/PyLith/trunk/pylith/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/pylith/Makefile.am	2007-06-08 23:11:50 UTC (rev 7116)
+++ short/3D/PyLith/trunk/pylith/Makefile.am	2007-06-09 01:30:56 UTC (rev 7117)
@@ -36,6 +36,7 @@
 	feassemble/ReferenceCell.py \
 	feassemble/quadrature/__init__.py \
 	feassemble/quadrature/Quadrature.py \
+	feassemble/quadrature/Quadrature0D.py \
 	feassemble/quadrature/Quadrature1D.py \
 	feassemble/quadrature/Quadrature1Din2D.py \
 	feassemble/quadrature/Quadrature1Din3D.py \

Added: short/3D/PyLith/trunk/pylith/feassemble/quadrature/Quadrature0D.py
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/quadrature/Quadrature0D.py	2007-06-08 23:11:50 UTC (rev 7116)
+++ short/3D/PyLith/trunk/pylith/feassemble/quadrature/Quadrature0D.py	2007-06-09 01:30:56 UTC (rev 7117)
@@ -0,0 +1,52 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+#                           Brad T. Aagaard
+#                        U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file pylith/feassemble/quadrature/Qudrature0D.py
+##
+## @brief Python object implementing 0-D integration in 1-D space
+## using numerical quadrature.
+##
+## Factory: quadrature
+
+from Quadrature import Quadrature
+
+# Quadrature0D class
+class Quadrature0D(Quadrature):
+  """
+  Python object for integrating over 0-D finite-elements in a 1-D
+  domain using quadrature.
+
+  Factory: quadrature.
+  """
+
+  def __init__(self, name="quadrature0d"):
+    """
+    Constructor.
+    """
+    Quadrature.__init__(self, name)
+    import pylith.feassemble.feassemble as bindings
+    self.cppHandle = bindings.Quadrature0D()
+    self.spaceDim = 1
+    self.cellDim = 0
+    return
+
+
+# FACTORIES ////////////////////////////////////////////////////////////
+
+def quadrature():
+  """
+  Factory associated with Quadrature0D.
+  """
+  return Quadrature0D()
+
+
+# End of file 

Modified: short/3D/PyLith/trunk/pylith/feassemble/quadrature/__init__.py
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/quadrature/__init__.py	2007-06-08 23:11:50 UTC (rev 7116)
+++ short/3D/PyLith/trunk/pylith/feassemble/quadrature/__init__.py	2007-06-09 01:30:56 UTC (rev 7117)
@@ -16,6 +16,7 @@
 ## initialization
 
 __all__ = ['Quadrature',
+           'Quadrature0D',
            'Quadrature1D',
            'Quadrature1Din2D',
            'Quadrature1Din3D',

Modified: short/3D/PyLith/trunk/unittests/libtests/faults/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/Makefile.am	2007-06-08 23:11:50 UTC (rev 7116)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/Makefile.am	2007-06-09 01:30:56 UTC (rev 7117)
@@ -26,16 +26,20 @@
 	TestFault.cc \
 	TestFaultCohesive.cc \
 	TestFaultCohesiveKin.cc \
-	TestFaultCohesiveKinLine2.cc \
+	TestFaultCohesiveKinTri3.cc \
 	test_faults.cc
 
+#	TestFaultCohesiveKinLine2.cc
+
+
 noinst_HEADERS = \
 	TestBruneSlipFn.hh \
 	TestEqKinSrc.hh \
 	TestFault.hh \
 	TestFaultCohesive.hh \
 	TestFaultCohesiveKin.hh \
-	TestFaultCohesiveKinLine2.hh
+	TestFaultCohesiveKinLine2.hh \
+	TestFaultCohesiveKinTri3.hh
 
 # Source files associated with testing data
 testfaults_SOURCES += \
@@ -49,7 +53,10 @@
 	data/CohesiveDataHex8.cc \
 	data/CohesiveDataHex8Lagrange.cc \
 	data/CohesiveDataTet4.cc \
-	data/CohesiveDataTet4Lagrange.cc
+	data/CohesiveDataTet4Lagrange.cc \
+	data/CohesiveKinData.cc \
+	data/CohesiveKinDataLine2.cc \
+	data/CohesiveKinDataTri3.cc
 
 noinst_HEADERS += \
 	data/CohesiveData.hh \
@@ -62,7 +69,10 @@
 	data/CohesiveDataHex8.hh \
 	data/CohesiveDataHex8Lagrange.hh \
 	data/CohesiveDataTet4.hh \
-	data/CohesiveDataTet4Lagrange.hh
+	data/CohesiveDataTet4Lagrange.hh \
+	data/CohesiveKinData.hh \
+	data/CohesiveKinDataLine2.hh \
+	data/CohesiveKinDataTri3.hh
 
 testfaults_LDFLAGS = $(PETSC_LIB) $(PYTHON_BLDLIBRARY)
 

Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc	2007-06-08 23:11:50 UTC (rev 7116)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc	2007-06-09 01:30:56 UTC (rev 7117)
@@ -16,13 +16,42 @@
 
 #include "pylith/faults/FaultCohesiveKin.hh" // USES FaultCohesiveKin
 
+#include "data/CohesiveKinData.hh" // USES CohesiveKinData
+
 #include "pylith/faults/EqKinSrc.hh" // USES EqKinSrc
+#include "pylith/faults/BruneSlipFn.hh" // USES BruneSlipFn
+#include "pylith/feassemble/Quadrature.hh" // USES Quadrature
+#include "pylith/topology/FieldsManager.hh" // USES FieldsManager
+#include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
 
-#include <stdexcept> // TEMPORARY
+#include "spatialdata/geocoords/CSCart.hh" // USES CSCart
+#include "spatialdata/spatialdb/SimpleDB.hh" // USES SimpleDB
+#include "spatialdata/spatialdb/SimpleIOAscii.hh" // USES SimpleIOAscii
+
+#include <stdexcept> // USES runtime_error
+
 // ----------------------------------------------------------------------
 CPPUNIT_TEST_SUITE_REGISTRATION( pylith::faults::TestFaultCohesiveKin );
 
 // ----------------------------------------------------------------------
+// Setup testing data.
+void
+pylith::faults::TestFaultCohesiveKin::setUp(void)
+{ // setUp
+  _data = 0;
+  _quadrature = 0;
+} // setUp
+
+// ----------------------------------------------------------------------
+// Tear down testing data.
+void
+pylith::faults::TestFaultCohesiveKin::tearDown(void)
+{ // tearDown
+  delete _data; _data = 0;
+  delete _quadrature; _quadrature = 0;
+} // tearDown
+
+// ----------------------------------------------------------------------
 // Test constructor.
 void
 pylith::faults::TestFaultCohesiveKin::testConstructor(void)
@@ -56,7 +85,44 @@
 void
 pylith::faults::TestFaultCohesiveKin::testInitialize(void)
 { // testInitialize
-  throw std::logic_error("Unit test not implemented.");
+  ALE::Obj<Mesh> mesh;
+  FaultCohesiveKin fault;
+  _initialize(&mesh, &fault);
+  
+  // Check set of constraint vertices
+  CPPUNIT_ASSERT_EQUAL(_data->numConstraintVert,
+		       int(fault._constraintVert.size()));
+  const std::set<Mesh::point_type>::const_iterator vertConstraintBegin = 
+    fault._constraintVert.begin();
+  const std::set<Mesh::point_type>::const_iterator vertConstraintEnd = 
+    fault._constraintVert.end();
+  int iVertex = 0;
+  for (std::set<Mesh::point_type>::const_iterator v_iter=vertConstraintBegin;
+       v_iter != vertConstraintEnd;
+       ++v_iter)
+    CPPUNIT_ASSERT_EQUAL(_data->constraintVertices[iVertex],
+			 *v_iter);
+
+  // Check orientation
+  iVertex = 0;
+  const int cellDim = _data->cellDim;
+  const int spaceDim = _data->spaceDim;
+  const int orientationSize = (cellDim > 0) ? cellDim*spaceDim : 1;
+  for (std::set<Mesh::point_type>::const_iterator v_iter=vertConstraintBegin;
+       v_iter != vertConstraintEnd;
+       ++v_iter) {
+    const int fiberDim = fault._orientation->getFiberDimension(*v_iter);
+    CPPUNIT_ASSERT_EQUAL(orientationSize, fiberDim);
+    const real_section_type::value_type* vertexOrient = 
+      fault._orientation->restrictPoint(*v_iter);
+    CPPUNIT_ASSERT(0 != vertexOrient);
+    const double tolerance = 1.0e-06;
+    for (int i=0; i < orientationSize; ++i) {
+      const int index = iVertex*orientationSize+i;
+      CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->orientation[index],
+				   vertexOrient[i], tolerance);
+    } // for
+  } // for
 } // testInitialize
 
 // ----------------------------------------------------------------------
@@ -64,7 +130,67 @@
 void
 pylith::faults::TestFaultCohesiveKin::testIntegrateResidual(void)
 { // testIntegrateResidual
-  throw std::logic_error("Unit test not implemented.");
+  ALE::Obj<Mesh> mesh;
+  FaultCohesiveKin fault;
+  _initialize(&mesh, &fault);
+
+  // Setup fields
+  topology::FieldsManager fields(mesh);
+  fields.addReal("residual");
+  fields.addReal("dispTpdt");
+  fields.addReal("dispT");
+  const char* history[] = { "dispTpdt", "dispT" };
+  const int historySize = 2;
+  fields.createHistory(history, historySize);
+  
+  const ALE::Obj<real_section_type>& residual = fields.getReal("residual");
+  CPPUNIT_ASSERT(!residual.isNull());
+  const int spaceDim = _data->spaceDim;
+  residual->setFiberDimension(mesh->depthStratum(0), spaceDim);
+  mesh->allocate(residual);
+  residual->zero();
+  fields.copyLayout("residual");
+
+  const ALE::Obj<real_section_type>& dispTpdt = fields.getReal("dispTpdt");
+  const ALE::Obj<real_section_type>& dispT = fields.getReal("dispT");
+  CPPUNIT_ASSERT(!dispTpdt.isNull());
+  CPPUNIT_ASSERT(!dispT.isNull());
+
+  const ALE::Obj<Mesh::label_sequence>& vertices = mesh->depthStratum(0);
+  CPPUNIT_ASSERT(!vertices.isNull());
+  const Mesh::label_sequence::iterator vBegin = vertices->begin();
+  const Mesh::label_sequence::iterator vEnd = vertices->end();
+  int iVertex = 0;
+  for (Mesh::label_sequence::iterator v_iter=vBegin;
+       v_iter != vEnd;
+       ++v_iter) {
+    dispTpdt->updatePoint(*v_iter, &_data->fieldTpdt[iVertex*spaceDim]);
+    dispT->updatePoint(*v_iter, &_data->fieldT[iVertex*spaceDim]);
+  } // for
+  
+  // Call integrateResidual()
+  fault.integrateResidual(residual, &fields, mesh);
+
+  // Check values
+  const double* valsE = _data->valsResidual;
+  iVertex = 0;
+  const int fiberDimE = spaceDim;
+  const double tolerance = 1.0e-06;
+  for (Mesh::label_sequence::iterator v_iter=vBegin;
+       v_iter != vEnd;
+       ++v_iter) {
+    const int fiberDim = residual->getFiberDimension(*v_iter);
+    CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
+    const real_section_type::value_type* vals = 
+      residual->restrictPoint(*v_iter);
+    for (int i=0; i < fiberDimE; ++i) {
+      const int index = iVertex*spaceDim+i;
+      if (valsE[index] > tolerance)
+	CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valsE[index], tolerance);
+      else
+	CPPUNIT_ASSERT_DOUBLES_EQUAL(valsE[index], vals[i], tolerance);
+    } // for
+  } // for
 } // testIntegrateResidual
 
 // ----------------------------------------------------------------------
@@ -72,7 +198,91 @@
 void
 pylith::faults::TestFaultCohesiveKin::testIntegrateJacobian(void)
 { // testIntegrateJacobian
-  throw std::logic_error("Unit test not implemented.");
+
+  ALE::Obj<Mesh> mesh;
+  FaultCohesiveKin fault;
+  _initialize(&mesh, &fault);
+
+  // Setup fields
+  topology::FieldsManager fields(mesh);
+  fields.addReal("residual");
+  fields.addReal("dispTpdt");
+  fields.addReal("dispT");
+  const char* history[] = { "dispTpdt", "dispT" };
+  const int historySize = 2;
+  fields.createHistory(history, historySize);
+  
+  const ALE::Obj<real_section_type>& residual = fields.getReal("residual");
+  CPPUNIT_ASSERT(!residual.isNull());
+  const int spaceDim = _data->spaceDim;
+  residual->setFiberDimension(mesh->depthStratum(0), spaceDim);
+  mesh->allocate(residual);
+  residual->zero();
+  fields.copyLayout("residual");
+
+  const ALE::Obj<real_section_type>& dispTpdt = fields.getReal("dispTpdt");
+  const ALE::Obj<real_section_type>& dispT = fields.getReal("dispT");
+  CPPUNIT_ASSERT(!dispTpdt.isNull());
+  CPPUNIT_ASSERT(!dispT.isNull());
+
+  const ALE::Obj<Mesh::label_sequence>& vertices = mesh->depthStratum(0);
+  CPPUNIT_ASSERT(!vertices.isNull());
+  const Mesh::label_sequence::iterator vBegin = vertices->begin();
+  const Mesh::label_sequence::iterator vEnd = vertices->end();
+  int iVertex = 0;
+  for (Mesh::label_sequence::iterator v_iter=vBegin;
+       v_iter != vEnd;
+       ++v_iter) {
+    dispTpdt->updatePoint(*v_iter, &_data->fieldTpdt[iVertex*spaceDim]);
+    dispT->updatePoint(*v_iter, &_data->fieldT[iVertex*spaceDim]);
+  } // for
+  
+  PetscMat jacobian;
+  PetscErrorCode err = MeshCreateMatrix(mesh, dispTpdt, MATMPIBAIJ, &jacobian);
+  CPPUNIT_ASSERT(0 == err);
+
+  fault.integrateJacobian(&jacobian, &fields, mesh);
+  CPPUNIT_ASSERT_EQUAL(false, fault.needNewJacobian());
+
+  err = MatAssemblyBegin(jacobian, MAT_FINAL_ASSEMBLY);
+  CPPUNIT_ASSERT(0 == err);
+  err = MatAssemblyEnd(jacobian, MAT_FINAL_ASSEMBLY);
+  CPPUNIT_ASSERT(0 == err);
+
+  const double* valsE = _data->valsJacobian;
+  const int nrowsE = dispT->sizeWithBC();
+  const int ncolsE = nrowsE;
+
+  int nrows = 0;
+  int ncols = 0;
+  MatGetSize(jacobian, &nrows, &ncols);
+  CPPUNIT_ASSERT_EQUAL(nrowsE, nrows);
+  CPPUNIT_ASSERT_EQUAL(ncolsE, ncols);
+
+  PetscMat jDense;
+  PetscMat jSparseAIJ;
+  MatConvert(jacobian, MATSEQAIJ, MAT_INITIAL_MATRIX, &jSparseAIJ);
+  MatConvert(jSparseAIJ, MATSEQDENSE, MAT_INITIAL_MATRIX, &jDense);
+
+  double_array vals(nrows*ncols);
+  int_array rows(nrows);
+  int_array cols(ncols);
+  for (int iRow=0; iRow < nrows; ++iRow)
+    rows[iRow] = iRow;
+  for (int iCol=0; iCol < ncols; ++iCol)
+    cols[iCol] = iCol;
+  MatGetValues(jDense, nrows, &rows[0], ncols, &cols[0], &vals[0]);
+  const double tolerance = 1.0e-06;
+  for (int iRow=0; iRow < nrows; ++iRow)
+    for (int iCol=0; iCol < ncols; ++iCol) {
+      const int index = ncols*iRow+iCol;
+      if (fabs(valsE[index]) > 1.0)
+	CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[index]/valsE[index], tolerance);
+      else
+	CPPUNIT_ASSERT_DOUBLES_EQUAL(valsE[index], vals[index], tolerance);
+    } // for
+  MatDestroy(jDense);
+  MatDestroy(jSparseAIJ);
 } // testIntegrateJacobian
 
 // ----------------------------------------------------------------------
@@ -85,7 +295,29 @@
   // formation does not eliminate any DOF from the system of
   // equations.
 
-  throw std::logic_error("Unit test not implemented.");
+  ALE::Obj<Mesh> mesh;
+  FaultCohesiveKin fault;
+  _initialize(&mesh, &fault);
+  
+  const int fiberDim = 3;
+  const ALE::Obj<real_section_type>& field = 
+    new real_section_type(mesh->comm(), mesh->debug());
+  CPPUNIT_ASSERT(!field.isNull());
+
+  field->setFiberDimension(mesh->depthStratum(0), fiberDim);
+  fault.setConstraintSizes(field, mesh);
+  mesh->allocate(field);
+
+  const ALE::Obj<Mesh::label_sequence>& vertices = mesh->depthStratum(0);
+  CPPUNIT_ASSERT(!vertices.isNull());
+  const Mesh::label_sequence::iterator vBegin = vertices->begin();
+  const Mesh::label_sequence::iterator vEnd = vertices->end();
+  for (Mesh::label_sequence::iterator v_iter=vBegin;
+       v_iter != vEnd;
+       ++v_iter) {
+    CPPUNIT_ASSERT_EQUAL(fiberDim, field->getFiberDimension(*v_iter));
+    CPPUNIT_ASSERT_EQUAL(0, field->getConstraintDimension(*v_iter));
+  } // for
 } // testSetConstraintSizes
 
 // ----------------------------------------------------------------------
@@ -93,7 +325,49 @@
 void
 pylith::faults::TestFaultCohesiveKin::testSetField(void)
 { // testSetField
-  throw std::logic_error("Unit test not implemented.");
+  ALE::Obj<Mesh> mesh;
+  FaultCohesiveKin fault;
+  _initialize(&mesh, &fault);
+
+  // Setup fields
+  topology::FieldsManager fields(mesh);
+  fields.addReal("disp");
+  
+  const ALE::Obj<real_section_type>& disp = fields.getReal("disp");
+  CPPUNIT_ASSERT(!disp.isNull());
+  const int spaceDim = _data->spaceDim;
+  disp->setFiberDimension(mesh->depthStratum(0), spaceDim);
+  mesh->allocate(disp);
+  disp->zero();
+
+  const double t = 2.134;
+  fault.setField(t, disp, mesh);
+
+  // Check values
+  const double* valsE = _data->valsSlip;
+  int iVertex = 0;
+  const int fiberDimE = spaceDim;
+  const double tolerance = 1.0e-06;
+
+  const ALE::Obj<Mesh::label_sequence>& vertices = mesh->depthStratum(0);
+  CPPUNIT_ASSERT(!vertices.isNull());
+  const Mesh::label_sequence::iterator vBegin = vertices->begin();
+  const Mesh::label_sequence::iterator vEnd = vertices->end();
+  for (Mesh::label_sequence::iterator v_iter=vBegin;
+       v_iter != vEnd;
+       ++v_iter) {
+    const int fiberDim = disp->getFiberDimension(*v_iter);
+    CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
+    const real_section_type::value_type* vals = 
+      disp->restrictPoint(*v_iter);
+    for (int i=0; i < fiberDimE; ++i) {
+      const int index = iVertex*spaceDim+i;
+      if (valsE[index] > tolerance)
+	CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valsE[index], tolerance);
+      else
+	CPPUNIT_ASSERT_DOUBLES_EQUAL(valsE[index], vals[i], tolerance);
+    } // for
+  } // for
 } // testSetField
 
 // ----------------------------------------------------------------------
@@ -102,6 +376,64 @@
 pylith::faults::TestFaultCohesiveKin::_initialize(ALE::Obj<ALE::Mesh>* mesh,
 					FaultCohesiveKin* const fault) const
 { // _initialize
+  CPPUNIT_ASSERT(0 != _data);
+  CPPUNIT_ASSERT(0 != mesh);
+  CPPUNIT_ASSERT(0 != fault);
+  CPPUNIT_ASSERT(0 != _quadrature);
+
+  try {
+    meshio::MeshIOAscii iohandler;
+    iohandler.filename(_data->meshFilename);
+    iohandler.read(mesh);
+    CPPUNIT_ASSERT(!mesh->isNull());
+    (*mesh)->getFactory()->clear();
+
+    spatialdata::geocoords::CSCart cs;
+    cs.setSpaceDim((*mesh)->getDimension());
+    cs.initialize();
+
+    _quadrature->initialize(_data->verticesRef, 
+			    _data->basis, _data->basisDeriv, _data->quadPts,
+			    _data->quadWts, _data->cellDim, _data->numBasis,
+			    _data->numQuadPts, _data->spaceDim);
+
+    // Setup earthquake source
+    spatialdata::spatialdb::SimpleDB dbFinalSlip("final slip");
+    spatialdata::spatialdb::SimpleIOAscii ioFinalSlip;
+    ioFinalSlip.filename(_data->finalSlipFilename);
+    dbFinalSlip.ioHandler(&ioFinalSlip);
+  
+    spatialdata::spatialdb::SimpleDB dbSlipTime("slip time");
+    spatialdata::spatialdb::SimpleIOAscii ioSlipTime;
+    ioSlipTime.filename(_data->slipTimeFilename);
+    dbSlipTime.ioHandler(&ioSlipTime);
+  
+    spatialdata::spatialdb::SimpleDB dbPeakRate("peak rate");
+    spatialdata::spatialdb::SimpleIOAscii ioPeakRate;
+    ioPeakRate.filename(_data->peakRateFilename);
+    dbPeakRate.ioHandler(&ioPeakRate);
+
+    BruneSlipFn slipFn;
+    slipFn.dbFinalSlip(&dbFinalSlip);
+    slipFn.dbSlipTime(&dbSlipTime);
+    slipFn.dbPeakRate(&dbPeakRate);
+  
+    EqKinSrc eqsrc;
+    eqsrc.slipfn(&slipFn);
+  
+    fault->id(_data->id);
+    fault->label(_data->label);
+    fault->quadrature(_quadrature);
+    fault->eqsrc(&eqsrc);
+    fault->adjustTopology(*mesh);
+
+    const double upDirVals[] = { 0.0, 0.0, 1.0 };
+    double_array upDir(upDirVals, 3);
+
+    fault->initialize(*mesh, &cs, upDir); 
+  } catch (const ALE::Exception& err) {
+    throw std::runtime_error(err.msg());
+  } // catch
 } // _initialize
 
 

Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.hh	2007-06-08 23:11:50 UTC (rev 7116)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.hh	2007-06-09 01:30:56 UTC (rev 7117)
@@ -31,8 +31,12 @@
     class TestFaultCohesiveKin;
 
     class FaultCohesiveKin; // USES FaultCohesiveKin
-    class FaultCohesiveKinData; // HOLDSA FaultCohesiveKinData
+    class CohesiveKinData; // HOLDSA CohesiveKinData
   } // faults
+
+  namespace feassemble {
+    class Quadrature; // HOLDSA Quadrature
+  } // feassemble
 } // pylith
 
 /// C++ unit testing for FaultCohesiveKin
@@ -51,11 +55,18 @@
   // PROTECTED MEMBERS //////////////////////////////////////////////////
 protected :
 
-  FaultCohesiveKinData* _data; ///< Data for testing
+  CohesiveKinData* _data; ///< Data for testing
+  feassemble::Quadrature* _quadrature; ///< Data used in testing
 
   // PUBLIC METHODS /////////////////////////////////////////////////////
 public :
 
+  /// Setup testing data.
+  void setUp(void);
+
+  /// Tear down testing data.
+  void tearDown(void);
+
   /// Test constructor.
   void testConstructor(void);
 

Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKinLine2.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKinLine2.cc	2007-06-08 23:11:50 UTC (rev 7116)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKinLine2.cc	2007-06-09 01:30:56 UTC (rev 7117)
@@ -14,8 +14,11 @@
 
 #include "TestFaultCohesiveKinLine2.hh" // Implementation of class methods
 
-//#include "data/FaultCohesiveKinDataLine2.hh" // USES DirichletDataLine2
+#include "data/CohesiveKinDataLine2.hh" // USES CohesiveKinDataLine2
 
+#include "pylith/feassemble/Quadrature0D.hh" // USES Quadrature1D
+#include "pylith/feassemble/GeometryLine1D.hh" // USES GeometryLine1D
+
 // ----------------------------------------------------------------------
 CPPUNIT_TEST_SUITE_REGISTRATION( pylith::faults::TestFaultCohesiveKinLine2 );
 
@@ -24,16 +27,12 @@
 void
 pylith::faults::TestFaultCohesiveKinLine2::setUp(void)
 { // setUp
-  //_data = new FaultCohesiveKinDataLine2();
+  _data = new CohesiveKinDataLine2();
+  _quadrature = new feassemble::Quadrature0D();
+  CPPUNIT_ASSERT(0 != _quadrature);
+  feassemble::GeometryLine1D geometry;
+  _quadrature->refGeometry(&geometry);
 } // setUp
 
-// ----------------------------------------------------------------------
-// Tear down testing data.
-void
-pylith::faults::TestFaultCohesiveKinLine2::tearDown(void)
-{ // tearDown
-  //delete _data;
-} // tearDown
 
-
 // End of file 

Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKinLine2.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKinLine2.hh	2007-06-08 23:11:50 UTC (rev 7116)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKinLine2.hh	2007-06-09 01:30:56 UTC (rev 7117)
@@ -51,9 +51,6 @@
   /// Setup testing data.
   void setUp(void);
 
-  /// Tear down testing data.
-  void tearDown(void);
-
 }; // class TestFaultCohesiveKinLine2
 
 #endif // pylith_faults_testfaultcohesiveline2_hh

Added: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKinTri3.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKinTri3.cc	2007-06-08 23:11:50 UTC (rev 7116)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKinTri3.cc	2007-06-09 01:30:56 UTC (rev 7117)
@@ -0,0 +1,38 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+//                           Brad T. Aagaard
+//                        U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "TestFaultCohesiveKinTri3.hh" // Implementation of class methods
+
+#include "data/CohesiveKinDataTri3.hh" // USES CohesiveKinDataTri3
+
+#include "pylith/feassemble/Quadrature1D.hh" // USES Quadrature1D
+#include "pylith/feassemble/GeometryLine1D.hh" // USES GeometryLine1D
+
+// ----------------------------------------------------------------------
+CPPUNIT_TEST_SUITE_REGISTRATION( pylith::faults::TestFaultCohesiveKinTri3 );
+
+// ----------------------------------------------------------------------
+// Setup testing data.
+void
+pylith::faults::TestFaultCohesiveKinTri3::setUp(void)
+{ // setUp
+  _data = new CohesiveKinDataTri3();
+  _quadrature = new feassemble::Quadrature1D();
+  CPPUNIT_ASSERT(0 != _quadrature);
+  feassemble::GeometryLine1D geometry;
+  _quadrature->refGeometry(&geometry);
+} // setUp
+
+
+// End of file 

Added: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKinTri3.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKinTri3.hh	2007-06-08 23:11:50 UTC (rev 7116)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKinTri3.hh	2007-06-09 01:30:56 UTC (rev 7117)
@@ -0,0 +1,59 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+//                           Brad T. Aagaard
+//                        U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/**
+ * @file unittests/libtests/faults/TestFaultCohesiveKinTri3.hh
+ *
+ * @brief C++ TestFaultCohesiveKinTri3 object.
+ *
+ * C++ unit testing for FaultCohesiveKin for mesh with 2-D triangular cells.
+ */
+
+#if !defined(pylith_faults_testfaultcohesivekintri3_hh)
+#define pylith_faults_testfaultcohesivekintri3_hh
+
+#include "TestFaultCohesiveKin.hh" // ISA TestFaultCohesiveKin
+
+/// Namespace for pylith package
+namespace pylith {
+  namespace faults {
+    class TestFaultCohesiveKinTri3;
+  } // bc
+} // pylith
+
+/// C++ unit testing for FaultCohesiveKin for mesh with 1-D line cells.
+class pylith::faults::TestFaultCohesiveKinTri3 : public TestFaultCohesiveKin
+{ // class TestFaultCohesiveKinTri3
+
+  // CPPUNIT TEST SUITE /////////////////////////////////////////////////
+  CPPUNIT_TEST_SUITE( TestFaultCohesiveKinTri3 );
+
+  CPPUNIT_TEST( testInitialize );
+  CPPUNIT_TEST( testIntegrateResidual );
+  CPPUNIT_TEST( testIntegrateJacobian );
+  CPPUNIT_TEST( testSetConstraintSizes );
+  CPPUNIT_TEST( testSetField );
+
+  CPPUNIT_TEST_SUITE_END();
+
+  // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+  /// Setup testing data.
+  void setUp(void);
+
+}; // class TestFaultCohesiveKinTri3
+
+#endif // pylith_faults_testfaultcohesivetri3_hh
+
+
+// End of file 

Added: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinData.cc	2007-06-08 23:11:50 UTC (rev 7116)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinData.cc	2007-06-09 01:30:56 UTC (rev 7117)
@@ -0,0 +1,51 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+//                           Brad T. Aagaard
+//                        U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#include "CohesiveKinData.hh"
+
+// ----------------------------------------------------------------------
+// Constructor
+pylith::faults::CohesiveKinData::CohesiveKinData(void) :
+  meshFilename(0),
+  spaceDim(0),
+  cellDim(0),
+  numBasis(0),
+  numQuadPts(0),
+  quadPts(0),
+  quadWts(0),
+  basis(0),
+  basisDeriv(0),
+  verticesRef(0),
+  id(0),
+  label(0),
+  finalSlipFilename(0),
+  slipTimeFilename(0),
+  peakRateFilename(0),
+  fieldTpdt(0),
+  fieldT(0),
+  orientation(0),
+  constraintVertices(0),
+  valsResidual(0),
+  valsSlip(0),
+  valsJacobian(0),
+  numConstraintVert(0)
+{ // constructor
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor
+pylith::faults::CohesiveKinData::~CohesiveKinData(void)
+{ // destructor
+} // destructor
+
+
+// End of file

Added: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinData.hh	2007-06-08 23:11:50 UTC (rev 7116)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinData.hh	2007-06-09 01:30:56 UTC (rev 7117)
@@ -0,0 +1,81 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+//                           Brad T. Aagaard
+//                        U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#if !defined(pylith_faults_cohesivekindata_hh)
+#define pylith_faults_cohesivekindata_hh
+
+namespace pylith {
+  namespace faults {
+     class CohesiveKinData;
+  } // pylith
+} // faults
+
+class pylith::faults::CohesiveKinData
+{
+
+// PUBLIC METHODS ///////////////////////////////////////////////////////
+public :
+  
+  /// Constructor
+  CohesiveKinData(void);
+
+  /// Destructor
+  ~CohesiveKinData(void);
+
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public:
+
+  char* meshFilename; ///< Filename for input mesh
+
+  /// @name Quadrature information
+  //@{
+  int spaceDim; ///< Number of dimensions in vertex coordinates
+  int cellDim; ///< Number of dimensions associated with cell
+  int numBasis; ///< Number of vertices in cell
+  int numQuadPts; ///< Number of quadrature points
+  double* quadPts; ///< Coordinates of quad pts in ref cell
+  double* quadWts; ///< Weights of quadrature points
+  double* basis; ///< Basis fns at quadrature points
+  double* basisDeriv; ///< Derivatives of basis fns at quad pts
+  double* verticesRef; ///< Coordinates of vertices in ref cell (dual basis)
+  //@}
+
+  /// @name Fault information
+  //@{
+  int id; ///< Fault material identifier
+  char* label; ///< Label for fault
+  char* finalSlipFilename; ///< Name of db for final slip
+  char* slipTimeFilename; ///< Name of db for slip time
+  char* peakRateFilename; ///< Name of db for peak rate
+  //@}
+
+  /// @name Input fields
+  //@{
+  double* fieldTpdt; ///< Input field at time t+dt.
+  double* fieldT; ///< Input field at time t.
+  //@}
+
+  /// @name Calculated values.
+  //@{
+  double* orientation; ///< Expected values for fault orientation.
+  int* constraintVertices; ///< Expected points for constraint vertices
+  double* valsResidual; ///< Expected values from residual calculation.
+  double* valsSlip; ///< Expected values from settting field.
+  double* valsJacobian; ///< Expected values from Jacobian calculation.
+  int numConstraintVert; ///< Number of constraint vertices
+  //@}
+
+};
+
+#endif // pylith_faults_cohesivekindata_hh
+
+// End of file

Added: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataLine2.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataLine2.cc	2007-06-08 23:11:50 UTC (rev 7116)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataLine2.cc	2007-06-09 01:30:56 UTC (rev 7117)
@@ -0,0 +1,153 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+//                           Brad T. Aagaard
+//                        U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+/* Original mesh
+ *
+ *   Cells are 0-1, vertices are 2-4.
+ *
+ *   2 -------- 3 -------- 4
+ *
+ * After adding cohesive elements
+ *   2 -------- 3-6-5 -------- 4
+ */
+
+#include "CohesiveKinDataLine2.hh"
+
+const char* pylith::faults::CohesiveKinDataLine2::_meshFilename =
+  "data/line2.mesh";
+
+const int pylith::faults::CohesiveKinDataLine2::_spaceDim = 1;
+
+const int pylith::faults::CohesiveKinDataLine2::_cellDim = 0;
+
+const int pylith::faults::CohesiveKinDataLine2::_numBasis = 1;
+
+const int pylith::faults::CohesiveKinDataLine2::_numQuadPts = 1;
+
+const double pylith::faults::CohesiveKinDataLine2::_quadPts[] = {
+  0.0,
+};
+
+const double pylith::faults::CohesiveKinDataLine2::_quadWts[] = {
+  1.0,
+};
+
+const double pylith::faults::CohesiveKinDataLine2::_basis[] = {
+  1.0,
+};
+
+const double pylith::faults::CohesiveKinDataLine2::_basisDeriv[] = {
+  1.0
+};
+
+const double pylith::faults::CohesiveKinDataLine2::_verticesRef[] = {
+  0.0
+};
+
+const int pylith::faults::CohesiveKinDataLine2::_id = 10;
+
+const char* pylith::faults::CohesiveKinDataLine2::_label = "fault";
+
+const char* pylith::faults::CohesiveKinDataLine2::_finalSlipFilename = 
+  "data/line2_finalslip.spatialdb";
+
+const char* pylith::faults::CohesiveKinDataLine2::_slipTimeFilename = 
+  "data/line2_sliptime.spatialdb";
+
+const char* pylith::faults::CohesiveKinDataLine2::_peakRateFilename = 
+  "data/line2_peakrate.spatialdb";
+
+const double pylith::faults::CohesiveKinDataLine2::_fieldTpdt[] = {
+  0.0,
+  0.0,
+  0.0,
+  0.0,
+  0.0
+};
+
+const double pylith::faults::CohesiveKinDataLine2::_fieldT[] = {
+  0.0,
+  0.0,
+  0.0,
+  0.0,
+  0.0
+};
+
+const int pylith::faults::CohesiveKinDataLine2::_numConstraintVert = 1;
+
+const double pylith::faults::CohesiveKinDataLine2::_orientation[] = {
+  1.0
+};
+
+const int pylith::faults::CohesiveKinDataLine2::_constraintVertices[] = {
+  6
+};
+
+const double pylith::faults::CohesiveKinDataLine2::_valsResidual[] = {
+  0.0,
+  0.0,
+  0.0,
+  0.0,
+  0.0,
+  0.0
+};
+
+const double pylith::faults::CohesiveKinDataLine2::_valsSlip[] = {
+  0.0,
+  0.0,
+  0.0,
+  0.0,
+  0.0,
+  1.2
+};
+
+const double pylith::faults::CohesiveKinDataLine2::_valsJacobian[] = {
+  0.0,
+  0.0,
+  0.0,
+  0.0,
+  0.0,
+  0.0
+};
+
+pylith::faults::CohesiveKinDataLine2::CohesiveKinDataLine2(void)
+{ // constructor
+  meshFilename = const_cast<char*>(_meshFilename);
+  spaceDim = _spaceDim;
+  cellDim = _cellDim;
+  numBasis = _numBasis;
+  numQuadPts = _numQuadPts;
+  quadPts = const_cast<double*>(_quadPts);
+  quadWts = const_cast<double*>(_quadWts);
+  basis = const_cast<double*>(_basis);
+  basisDeriv = const_cast<double*>(_basisDeriv);
+  verticesRef = const_cast<double*>(_verticesRef);
+  id = _id;
+  label = const_cast<char*>(_label);
+  finalSlipFilename = const_cast<char*>(_finalSlipFilename);
+  slipTimeFilename = const_cast<char*>(_slipTimeFilename);
+  peakRateFilename = const_cast<char*>(_peakRateFilename);
+  fieldTpdt = const_cast<double*>(_fieldTpdt);
+  fieldT = const_cast<double*>(_fieldT);
+  orientation = const_cast<double*>(_orientation);
+  constraintVertices = const_cast<int*>(_constraintVertices);
+  valsResidual = const_cast<double*>(_valsResidual);
+  valsSlip = const_cast<double*>(_valsSlip);
+  valsJacobian = const_cast<double*>(_valsJacobian);
+  numConstraintVert = _numConstraintVert;  
+} // constructor
+
+pylith::faults::CohesiveKinDataLine2::~CohesiveKinDataLine2(void)
+{}
+
+
+// End of file

Added: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataLine2.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataLine2.hh	2007-06-08 23:11:50 UTC (rev 7116)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataLine2.hh	2007-06-09 01:30:56 UTC (rev 7117)
@@ -0,0 +1,74 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+//                           Brad T. Aagaard
+//                        U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#if !defined(pylith_faults_cohesivekindataline2_hh)
+#define pylith_faults_cohesivekindataline2_hh
+
+#include "CohesiveKinData.hh"
+
+namespace pylith {
+  namespace faults {
+     class CohesiveKinDataLine2;
+  } // pylith
+} // faults
+
+class pylith::faults::CohesiveKinDataLine2 : public CohesiveKinData
+{
+
+// PUBLIC METHODS ///////////////////////////////////////////////////////
+public: 
+
+  /// Constructor
+  CohesiveKinDataLine2(void);
+
+  /// Destructor
+  ~CohesiveKinDataLine2(void);
+
+// PRIVATE MEMBERS //////////////////////////////////////////////////////
+private:
+
+  static const char* _meshFilename; ///< Filename of input mesh
+
+  static const int _spaceDim; ///< Number of dimensions in vertex coordinates
+  static const int _cellDim; ///< Number of dimensions associated with cell
+
+  static const int _numBasis; ///< Number of vertices in cell
+  static const int _numQuadPts; ///< Number of quadrature points
+  static const double _quadPts[]; ///< Coordinates of quad pts in ref cell
+  static const double _quadWts[]; ///< Weights of quadrature points
+  static const double _basis[]; ///< Basis fns at quadrature points
+  static const double _basisDeriv[]; ///< Derivatives of basis fns at quad pts
+  static const double _verticesRef[]; ///< Coordinates of vertices in ref cell (dual basis)
+
+  static const int _id; ///< Fault material identifier
+  static const char* _label; ///< Label for fault
+  static const char* _finalSlipFilename; ///< Name of db for final slip
+  static const char* _slipTimeFilename; ///< Name of db for slip time
+  static const char* _peakRateFilename; ///< Name of db for peak rate
+  //@}
+
+  static const double _fieldTpdt[]; ///< Input field at time t+dt.
+  static const double _fieldT[]; ///< Input field at time t.
+
+  static const double _orientation[]; ///< Expected values for fault orientation.
+  static const int _constraintVertices[]; ///< Expected points for constraint vertices
+  static const double _valsResidual[]; ///< Expected values from residual calculation.
+  static const double _valsSlip[]; ///< Expected values from settting field.
+  static const double _valsJacobian[]; ///< Expected values from Jacobian calculation.
+  static const int _numConstraintVert; ///< Number of constraint vertices
+
+};
+
+#endif // pylith_faults_cohesivekindataline2_hh
+
+
+// End of file

Added: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTri3.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTri3.cc	2007-06-08 23:11:50 UTC (rev 7116)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTri3.cc	2007-06-09 01:30:56 UTC (rev 7117)
@@ -0,0 +1,192 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+//                           Brad T. Aagaard
+//                        U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+/* Original mesh
+ *
+ * Cells are 0-1, vertices are 2-5.
+ *
+ *              3
+ *             /|\
+ *            / | \
+ *           /  |  \
+ *          /   |   \
+ *         2    |    5
+ *          \   |   /
+ *           \  |  /
+ *            \ | /
+ *             \|/
+ *              4
+ *
+ *
+ * After adding cohesive elements
+ *
+ * Cells are 0-1, 8, vertices are 2-7.
+ *
+ *              3 -7- 6
+ *             /|     |\
+ *            / |     | \
+ *           /  |     |  \
+ *          /   |     |   \
+ *         2    |     |    5
+ *          \   |     |   /
+ *           \  |     |  /
+ *            \ |     | /
+ *             \|     |/
+ *              4 -9- 8
+ */
+
+#include "CohesiveKinDataTri3.hh"
+
+const char* pylith::faults::CohesiveKinDataTri3::_meshFilename =
+  "data/tri3.mesh";
+
+const int pylith::faults::CohesiveKinDataTri3::_spaceDim = 2;
+
+const int pylith::faults::CohesiveKinDataTri3::_cellDim = 1;
+
+const int pylith::faults::CohesiveKinDataTri3::_numBasis = 2;
+
+const int pylith::faults::CohesiveKinDataTri3::_numQuadPts = 1;
+
+const double pylith::faults::CohesiveKinDataTri3::_quadPts[] = {
+  0.0,
+};
+
+const double pylith::faults::CohesiveKinDataTri3::_quadWts[] = {
+  2.0,
+};
+
+const double pylith::faults::CohesiveKinDataTri3::_basis[] = {
+  0.5,
+  0.5
+};
+
+const double pylith::faults::CohesiveKinDataTri3::_basisDeriv[] = {
+  -0.5,
+   0.5
+};
+
+const double pylith::faults::CohesiveKinDataTri3::_verticesRef[] = {
+  -1.0, 1.0
+};
+
+const int pylith::faults::CohesiveKinDataTri3::_id = 10;
+
+const char* pylith::faults::CohesiveKinDataTri3::_label = "fault";
+
+const char* pylith::faults::CohesiveKinDataTri3::_finalSlipFilename = 
+  "data/tri3_finalslip.spatialdb";
+
+const char* pylith::faults::CohesiveKinDataTri3::_slipTimeFilename = 
+  "data/tri3_sliptime.spatialdb";
+
+const char* pylith::faults::CohesiveKinDataTri3::_peakRateFilename = 
+  "data/tri3_peakrate.spatialdb";
+
+const double pylith::faults::CohesiveKinDataTri3::_fieldTpdt[] = {
+  0.0, 0.0,
+  0.0, 0.0,
+  0.0, 0.0,
+  0.0, 0.0,
+  0.0, 0.0,
+  0.0, 0.0,
+  0.0, 0.0,
+  0.0, 0.0,
+};
+
+const double pylith::faults::CohesiveKinDataTri3::_fieldT[] = {
+  0.0, 0.0,
+  0.0, 0.0,
+  0.0, 0.0,
+  0.0, 0.0,
+  0.0, 0.0,
+  0.0, 0.0,
+  0.0, 0.0,
+  0.0, 0.0,
+};
+
+const int pylith::faults::CohesiveKinDataTri3::_numConstraintVert = 2;
+
+const double pylith::faults::CohesiveKinDataTri3::_orientation[] = {
+  0.0, 1.0,
+  0.0, 1.0
+};
+
+const int pylith::faults::CohesiveKinDataTri3::_constraintVertices[] = {
+  7, 9
+};
+
+const double pylith::faults::CohesiveKinDataTri3::_valsResidual[] = {
+  0.0, 0.0,
+  0.0, 0.0,
+  0.0, 0.0,
+  0.0, 0.0,
+  0.0, 0.0,
+  0.0, 0.0,
+  0.0, 0.0,
+  0.0, 0.0,
+};
+
+const double pylith::faults::CohesiveKinDataTri3::_valsSlip[] = {
+  0.0, 0.0,
+  0.0, 0.0,
+  0.0, 0.0,
+  0.0, 0.0,
+  0.0, 0.0,
+  0.0, 0.0,
+  0.0, 0.0,
+  0.0, 0.0,
+};
+
+const double pylith::faults::CohesiveKinDataTri3::_valsJacobian[] = {
+  0.0, 0.0,
+  0.0, 0.0,
+  0.0, 0.0,
+  0.0, 0.0,
+  0.0, 0.0,
+  0.0, 0.0,
+  0.0, 0.0,
+  0.0, 0.0,
+};
+
+pylith::faults::CohesiveKinDataTri3::CohesiveKinDataTri3(void)
+{ // constructor
+  meshFilename = const_cast<char*>(_meshFilename);
+  spaceDim = _spaceDim;
+  cellDim = _cellDim;
+  numBasis = _numBasis;
+  numQuadPts = _numQuadPts;
+  quadPts = const_cast<double*>(_quadPts);
+  quadWts = const_cast<double*>(_quadWts);
+  basis = const_cast<double*>(_basis);
+  basisDeriv = const_cast<double*>(_basisDeriv);
+  verticesRef = const_cast<double*>(_verticesRef);
+  id = _id;
+  label = const_cast<char*>(_label);
+  finalSlipFilename = const_cast<char*>(_finalSlipFilename);
+  slipTimeFilename = const_cast<char*>(_slipTimeFilename);
+  peakRateFilename = const_cast<char*>(_peakRateFilename);
+  fieldTpdt = const_cast<double*>(_fieldTpdt);
+  fieldT = const_cast<double*>(_fieldT);
+  orientation = const_cast<double*>(_orientation);
+  constraintVertices = const_cast<int*>(_constraintVertices);
+  valsResidual = const_cast<double*>(_valsResidual);
+  valsSlip = const_cast<double*>(_valsSlip);
+  valsJacobian = const_cast<double*>(_valsJacobian);
+  numConstraintVert = _numConstraintVert;  
+} // constructor
+
+pylith::faults::CohesiveKinDataTri3::~CohesiveKinDataTri3(void)
+{}
+
+
+// End of file

Added: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTri3.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTri3.hh	2007-06-08 23:11:50 UTC (rev 7116)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTri3.hh	2007-06-09 01:30:56 UTC (rev 7117)
@@ -0,0 +1,74 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+//                           Brad T. Aagaard
+//                        U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#if !defined(pylith_faults_cohesivekindatatri3_hh)
+#define pylith_faults_cohesivekindatatri3_hh
+
+#include "CohesiveKinData.hh"
+
+namespace pylith {
+  namespace faults {
+     class CohesiveKinDataTri3;
+  } // pylith
+} // faults
+
+class pylith::faults::CohesiveKinDataTri3 : public CohesiveKinData
+{
+
+// PUBLIC METHODS ///////////////////////////////////////////////////////
+public: 
+
+  /// Constructor
+  CohesiveKinDataTri3(void);
+
+  /// Destructor
+  ~CohesiveKinDataTri3(void);
+
+// PRIVATE MEMBERS //////////////////////////////////////////////////////
+private:
+
+  static const char* _meshFilename; ///< Filename of input mesh
+
+  static const int _spaceDim; ///< Number of dimensions in vertex coordinates
+  static const int _cellDim; ///< Number of dimensions associated with cell
+
+  static const int _numBasis; ///< Number of vertices in cell
+  static const int _numQuadPts; ///< Number of quadrature points
+  static const double _quadPts[]; ///< Coordinates of quad pts in ref cell
+  static const double _quadWts[]; ///< Weights of quadrature points
+  static const double _basis[]; ///< Basis fns at quadrature points
+  static const double _basisDeriv[]; ///< Derivatives of basis fns at quad pts
+  static const double _verticesRef[]; ///< Coordinates of vertices in ref cell (dual basis)
+
+  static const int _id; ///< Fault material identifier
+  static const char* _label; ///< Label for fault
+  static const char* _finalSlipFilename; ///< Name of db for final slip
+  static const char* _slipTimeFilename; ///< Name of db for slip time
+  static const char* _peakRateFilename; ///< Name of db for peak rate
+  //@}
+
+  static const double _fieldTpdt[]; ///< Input field at time t+dt.
+  static const double _fieldT[]; ///< Input field at time t.
+
+  static const double _orientation[]; ///< Expected values for fault orientation.
+  static const int _constraintVertices[]; ///< Expected points for constraint vertices
+  static const double _valsResidual[]; ///< Expected values from residual calculation.
+  static const double _valsSlip[]; ///< Expected values from settting field.
+  static const double _valsJacobian[]; ///< Expected values from Jacobian calculation.
+  static const int _numConstraintVert; ///< Number of constraint vertices
+
+};
+
+#endif // pylith_faults_cohesivekindatatri3_hh
+
+
+// End of file

Modified: short/3D/PyLith/trunk/unittests/libtests/faults/test_faults.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/test_faults.cc	2007-06-08 23:11:50 UTC (rev 7116)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/test_faults.cc	2007-06-09 01:30:56 UTC (rev 7117)
@@ -11,6 +11,7 @@
 //
 
 #include "petsc.h"
+#include "petscmesh.h"
 
 #include <cppunit/extensions/TestFactoryRegistry.h>
 

Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/Makefile.am	2007-06-08 23:11:50 UTC (rev 7116)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/Makefile.am	2007-06-09 01:30:56 UTC (rev 7117)
@@ -51,6 +51,7 @@
 	TestGeometryHex3D.cc \
 	TestIntegrator.cc \
 	TestQuadrature.cc \
+	TestQuadrature0D.cc \
 	TestQuadrature1D.cc \
 	TestQuadrature1Din2D.cc \
 	TestQuadrature2D.cc \
@@ -89,6 +90,7 @@
 	TestGeometryHex3D.hh \
 	TestIntegrator.hh \
 	TestQuadrature.hh \
+	TestQuadrature0D.hh \
 	TestQuadrature1D.hh \
 	TestQuadrature1Din2D.hh \
 	TestQuadrature1Din3D.hh \

Added: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature0D.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature0D.cc	2007-06-08 23:11:50 UTC (rev 7116)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature0D.cc	2007-06-09 01:30:56 UTC (rev 7117)
@@ -0,0 +1,98 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+//                           Brad T. Aagaard
+//                        U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "TestQuadrature0D.hh" // Implementation of class methods
+
+#include "pylith/feassemble/Quadrature0D.hh"
+
+// ----------------------------------------------------------------------
+CPPUNIT_TEST_SUITE_REGISTRATION( pylith::feassemble::TestQuadrature0D );
+
+// ----------------------------------------------------------------------
+// Test constructor
+void
+pylith::feassemble::TestQuadrature0D::testConstructor(void)
+{ // testConstructor
+  Quadrature0D quadrature;
+} // testConstructor
+
+// ----------------------------------------------------------------------
+// Test computeGeometry().
+void
+pylith::feassemble::TestQuadrature0D::testPoint(void)
+{ // testPoint
+  const int cellDim = 0;
+  const int numBasis = 1;
+  const int numQuadPts = 1;
+  const int spaceDim = 1;
+  const double verticesRef[] = { 0.0 };
+  const double basis[] = { 1.0 };
+  const double basisDeriv[] = { 1.0 };
+  const double quadPtsRef[] = { 0.0 };
+  const double quadWts[] = { 1.0 };
+
+  const int numVertices = 1;
+  const int numCells = 1;
+  const double vertCoords[] = { 1.1 };
+  const int cells[] = { 0 };
+  const double quadPts[] = { 1.1 };
+  const double jacobian[] = { 1.0 };
+  const double jacobianInv[] = { 1.0 };
+  const double jacobianDet[] = { 1.0 };
+
+  const double minJacobian = 1.0e-06;
+  
+  Quadrature0D quadrature;
+
+  quadrature.minJacobian(minJacobian);
+  quadrature.initialize(verticesRef, basis, basisDeriv, quadPtsRef, quadWts,
+		    cellDim, numBasis, numQuadPts, spaceDim);
+
+  // Create mesh with test cell
+  ALE::Obj<Mesh> mesh = new Mesh(PETSC_COMM_WORLD, cellDim);
+  CPPUNIT_ASSERT(!mesh.isNull());
+  ALE::Obj<sieve_type> sieve = new sieve_type(mesh->comm());
+  CPPUNIT_ASSERT(!sieve.isNull());
+
+  const bool interpolate = false;
+  ALE::SieveBuilder<Mesh>::buildTopology(sieve, cellDim, numCells,
+		     (int*) cells, numVertices, interpolate, numBasis);
+  mesh->setSieve(sieve);
+  mesh->stratify();
+  ALE::SieveBuilder<Mesh>::buildCoordinates(mesh, spaceDim, vertCoords);
+  
+  // Check values from computeGeometry()
+  const ALE::Obj<Mesh::label_sequence>& cellsMesh = mesh->heightStratum(0);
+  CPPUNIT_ASSERT(!cellsMesh.isNull());
+  const Mesh::label_sequence::iterator e_iter = cellsMesh->begin(); 
+  const ALE::Obj<Mesh::real_section_type>& coordinates = 
+    mesh->getRealSection("coordinates");
+  CPPUNIT_ASSERT(!coordinates.isNull());
+  quadrature.computeGeometry(mesh, coordinates, *e_iter);
+
+  CPPUNIT_ASSERT(1 == numCells);
+
+  const double tolerance = 1.0e-06;
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(quadPts[0], quadrature._quadPts[0], 
+			       tolerance);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(jacobian[0], quadrature._jacobian[0], 
+			       tolerance);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(jacobianInv[0], quadrature._jacobianInv[0], 
+			       tolerance);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(jacobianDet[0], quadrature._jacobianDet[0], 
+			       tolerance);
+} // testPoint
+
+
+// End of file 

Added: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature0D.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature0D.hh	2007-06-08 23:11:50 UTC (rev 7116)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature0D.hh	2007-06-09 01:30:56 UTC (rev 7117)
@@ -0,0 +1,56 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+//                           Brad T. Aagaard
+//                        U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/**
+ * @file unittests/libtests/feassemble/TestQuadrature0D.hh
+ *
+ * @brief C++ TestQuadrature0D object
+ *
+ * C++ unit testing for Quadrature0D.
+ */
+
+#if !defined(pylith_feassemble_testquadrature0d_hh)
+#define pylith_feassemble_testquadrature0d_hh
+
+#include "TestQuadrature.hh"
+
+/// Namespace for pylith package
+namespace pylith {
+  namespace feassemble {
+    class TestQuadrature0D;
+  } // feassemble
+} // pylith
+
+/// C++ unit testing for Quadrature0D
+class pylith::feassemble::TestQuadrature0D : public TestQuadrature
+{ // class TestQuadrature0D
+
+  // CPPUNIT TEST SUITE /////////////////////////////////////////////////
+  CPPUNIT_TEST_SUITE( TestQuadrature0D );
+  CPPUNIT_TEST( testConstructor );
+  CPPUNIT_TEST( testPoint );
+  CPPUNIT_TEST_SUITE_END();
+
+  // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+  /// Test constructor
+  void testConstructor(void);
+
+  /// Test initialize() & computeGeometry().
+  void testPoint(void);
+
+}; // class TestQuadrature0D
+
+#endif // pylith_feassemble_testquadrature0d_hh
+
+// End of file 

Modified: short/3D/PyLith/trunk/unittests/pytests/feassemble/TestQuadrature.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/feassemble/TestQuadrature.py	2007-06-08 23:11:50 UTC (rev 7116)
+++ short/3D/PyLith/trunk/unittests/pytests/feassemble/TestQuadrature.py	2007-06-09 01:30:56 UTC (rev 7117)
@@ -17,6 +17,7 @@
 import unittest
 import numpy
 
+from pylith.feassemble.quadrature.Quadrature0D import Quadrature0D
 from pylith.feassemble.quadrature.Quadrature1D import Quadrature1D
 from pylith.feassemble.quadrature.Quadrature1Din2D import Quadrature1Din2D
 from pylith.feassemble.quadrature.Quadrature1Din3D import Quadrature1Din3D
@@ -122,6 +123,11 @@
     """
     Test constructors for quadrature objects.
     """
+    q = Quadrature0D()
+    self.assertEqual(1, q.spaceDim)
+    self.assertEqual(0, q.cellDim)
+    self.failIfEqual(None, q.cppHandle)
+    
     q = Quadrature1D()
     self.assertEqual(1, q.spaceDim)
     self.assertEqual(1, q.cellDim)