[cig-commits] r4567 - in short/3D/PyLith/trunk: . libsrc/feassemble libsrc/meshio unittests unittests/libtests unittests/libtests/feassemble

baagaard at geodynamics.org baagaard at geodynamics.org
Sun Sep 17 20:35:27 PDT 2006


Author: baagaard
Date: 2006-09-17 20:35:26 -0700 (Sun, 17 Sep 2006)
New Revision: 4567

Added:
   short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.icc
   short/3D/PyLith/trunk/unittests/
   short/3D/PyLith/trunk/unittests/Makefile.am
   short/3D/PyLith/trunk/unittests/libtests/
   short/3D/PyLith/trunk/unittests/libtests/Makefile.am
   short/3D/PyLith/trunk/unittests/libtests/feassemble/
   short/3D/PyLith/trunk/unittests/libtests/feassemble/Makefile.am
   short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.cc
   short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.hh
   short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature1D.cc
   short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature1D.hh
   short/3D/PyLith/trunk/unittests/libtests/feassemble/test_feassemble.cc
Modified:
   short/3D/PyLith/trunk/Makefile.am
   short/3D/PyLith/trunk/configure.ac
   short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am
   short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc
   short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh
   short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.icc
   short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.cc
   short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.hh
   short/3D/PyLith/trunk/libsrc/meshio/MeshIO.cc
   short/3D/PyLith/trunk/libsrc/meshio/MeshIOAscii.cc
Log:
Worked on quadrature implementation. Added Quadrature1D (1-D quadrature in 1-D domain). Started unit testing of quadrature.

Modified: short/3D/PyLith/trunk/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/Makefile.am	2006-09-18 03:30:20 UTC (rev 4566)
+++ short/3D/PyLith/trunk/Makefile.am	2006-09-18 03:35:26 UTC (rev 4567)
@@ -13,12 +13,14 @@
 ACLOCAL_AMFLAGS = -I ./m4
 
 SUBDIRS = \
-	pylith \
 	libsrc \
 	modulesrc \
+	pylith \
 	applications
 
-#modulesrc
+if ENABLE_TESTING
+  SUBDIRS += unittests
+endif
 
 if ENABLE_DOCUMENTATION
   SUBDIRS += doc

Modified: short/3D/PyLith/trunk/configure.ac
===================================================================
--- short/3D/PyLith/trunk/configure.ac	2006-09-18 03:30:20 UTC (rev 4566)
+++ short/3D/PyLith/trunk/configure.ac	2006-09-18 03:35:26 UTC (rev 4567)
@@ -131,6 +131,9 @@
                 modulesrc/Makefile
                 modulesrc/meshio/Makefile
 		applications/Makefile
+		unittests/Makefile
+		unittests/libtests/Makefile
+		unittests/libtests/feassemble/Makefile
                 doc/Makefile])
 
 #                 libsrc/materials/Makefile

Modified: short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am	2006-09-18 03:30:20 UTC (rev 4566)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am	2006-09-18 03:35:26 UTC (rev 4567)
@@ -16,10 +16,14 @@
 lib_LTLIBRARIES = libpylithfeassemble.la
 
 libpylithfeassemble_la_SOURCES = \
-	Quadrature.cc
+	Quadrature.cc \
+	Quadrature1D.cc
 
 subpkginclude_HEADERS = \
-	Quadrature.hh
+	Quadrature.hh \
+	Quadrature.icc \
+	Quadrature1D.hh \
+	Quadrature1D.icc
 
 noinst_HEADERS =
 

Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc	2006-09-18 03:30:20 UTC (rev 4566)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc	2006-09-18 03:35:26 UTC (rev 4567)
@@ -10,19 +10,31 @@
 // ======================================================================
 //
 
+#include <portinfo>
+
 #include "Quadrature.hh" // implementation of class methods
 
+#include <string.h> // USES memcpy()
+#include <assert.h> // USES assert()
 #include <stdexcept> // USES std::runtime_error
+#include <sstream> // USES std::ostringstream
 
 // ----------------------------------------------------------------------
 // Constructor
 pylith::feassemble::Quadrature::Quadrature(void) :
-  _pBasisFns(0),
-  _pBasisFnsDeriv(0),
-  _pQuadPts(0),
-  _pQuadWts(0),
-  _numDims(0),
-  _numCorners(0)
+  _jacobianTol(0),
+  _basis(0),
+  _basisDeriv(0),
+  _quadPtsRef(0),
+  _quadPts(0),
+  _quadWts(0),
+  _jacobian(0),
+  _jacobianInv(0),
+  _jacobianDet(0),
+  _cellDim(0),
+  _numCorners(0),
+  _numQuadPts(0),
+  _spaceDim(0)
 { // constructor
 } // constructor
 
@@ -30,59 +42,156 @@
 // Destructor
 pylith::feassemble::Quadrature::~Quadrature(void)
 { // destructor
-  delete[] _pBasisFns; _pBasisFns = 0;
-  delete[] _pBasisFnsDeriv; _pBasisFnsDeriv = 0;
-  delete[] _pQuadPts; _pQuadPts = 0;
-  delete[] _pQuadWts; _pQuadWts = 0;
+  delete[] _basis; _basis = 0;
+  delete[] _basisDeriv; _basisDeriv = 0;
+  delete[] _quadPtsRef; _quadPtsRef = 0;
+  delete[] _quadPts; _quadPts = 0;
+  delete[] _quadWts; _quadWts = 0;
+  delete[] _jacobian; _jacobian = 0;
+  delete[] _jacobianInv; _jacobianInv = 0;
+  delete[] _jacobianDet; _jacobianDet = 0;
 } // destructor
   
 // ----------------------------------------------------------------------
+// Copy constructor
+pylith::feassemble::Quadrature::Quadrature(const Quadrature& q) :
+  _jacobianTol(q._jacobianTol),
+  _basis(0),
+  _basisDeriv(0),
+  _quadPtsRef(0),
+  _quadPts(0),
+  _quadWts(0),
+  _jacobian(0),
+  _jacobianInv(0),
+  _jacobianDet(0),
+  _cellDim(q._cellDim),
+  _numCorners(q._numCorners),
+  _numQuadPts(q._numQuadPts),
+  _spaceDim(q._spaceDim)
+{ // copy constructor
+  if (0 != q._basis) {
+    const int size = _numCorners * _numQuadPts;
+    _basis = (size > 0) ? new double[size] : 0;
+    memcpy(_basis, q._basis, size*sizeof(double));
+  } // if
+
+  if (0 != q._basisDeriv) {
+    const int size = _numCorners * _numQuadPts * _spaceDim;
+    _basisDeriv = (size > 0) ? new double[size] : 0;
+    memcpy(_basisDeriv, q._basisDeriv, size*sizeof(double));
+  } // if
+
+  if (0 != q._quadPtsRef) {
+    const int size = _numQuadPts * _cellDim;
+    _quadPtsRef = (size > 0) ? new double[size] : 0;
+    memcpy(_quadPtsRef, q._quadPtsRef, size*sizeof(double));
+  } // if
+
+  if (0 != q._quadPts) {
+    const int size = _numQuadPts*_spaceDim;
+    _quadPts = (size > 0) ? new double[size] : 0;
+    memcpy(_quadPts, q._quadPts, size*sizeof(double));
+  } // if
+
+  if (0 != q._quadWts) {
+    const int size = _numQuadPts;
+    _quadWts = (size > 0) ? new double[size] : 0;
+    memcpy(_quadWts, q._quadWts, size*sizeof(double));
+  } // if
+  
+  if (0 != q._jacobian) {
+    const int size = _numQuadPts*_cellDim*_spaceDim;
+    _jacobian = (size > 0) ? new double[size] : 0;
+    memcpy(_jacobian, q._jacobian, size*sizeof(double));
+  } // if
+
+  if (0 != q._jacobianInv) {
+    const int size = _numQuadPts*_cellDim*_spaceDim;
+    _jacobianInv = (size > 0) ? new double[size] : 0;
+    memcpy(_jacobianInv, q._jacobianInv, size*sizeof(double));
+  } // if
+
+  if (0 != q._jacobianDet) {
+    const int size = _numQuadPts;
+    _jacobianDet = (size > 0) ? new double[size] : 0;
+    memcpy(_jacobianDet, q._jacobianDet, size*sizeof(double));
+  } // if
+} // copy constructor
+
+// ----------------------------------------------------------------------
 // Set basis functions and their derivatives and coordinates and
 //   weights of the quadrature points.
 void
-pylith::feassemble::Quadrature::initialize(const double* pBasisFns,
-					   const double* pBasisFnsDeriv,
-					   const double* pQuadPts,
-					   const double* pQuadWts,
-					   const int numDims,
+pylith::feassemble::Quadrature::initialize(const double* basis,
+					   const double* basisDeriv,
+					   const double* quadPtsRef,
+					   const double* quadWts,
+					   const int cellDim,
 					   const int numCorners,
-					   const int numQuadPts)
+					   const int numQuadPts,
+					   const int spaceDim)
 { // initialize
-  if (0 == pBasisFns ||
-      0 == pBasisFnsDeriv ||
-      0 == pQuadPts ||
-      0 == pQuadWts ||
-      numDims < 1 || numDims > 3)
-    throw std::runtime_error("Basis functions and their derivatives "
-			     "and quadrature points must all be specified.");
+  if (0 == basis ||
+      0 == basisDeriv ||
+      0 == quadPtsRef ||
+      0 == quadWts ||
+      cellDim < 1 || cellDim > 3 ||
+      numCorners < 1 ||
+      numQuadPts < 1 ||
+      spaceDim < 1 || spaceDim > 3) {
+    std::ostringstream msg;
+    msg << "Incompatible values for quadrature information. Basis functions,\n"
+	<< "their derivatives, and coordinates and weights of quadrature\n"
+	<< "points must all be specified.\n"
+	<< "Values:\n"
+	<< "  basis pointer: " << basis << "\n"
+	<< "  basis derivatites pointer: " << basisDeriv << "\n"
+	<< "  quadrature points pointer: " << quadPtsRef << "\n"
+	<< "  quadrature weights pointer: " << quadWts << "\n"
+	<< "  space dimension: " << spaceDim << "\n"
+	<< "  # vertices per cell: " << numCorners << "\n"
+	<< "  # quadrature points: " << numQuadPts << "\n"
+	<< "  dimension of coordinate space: " << spaceDim << "\n";
+    throw std::runtime_error(msg.str());
+  } // if
 
   int size = numCorners * numQuadPts;
-  delete[] _pBasisFns;
-  _pBasisFns = (size > 0) ? new double[size] : 0;
+  delete[] _basis; _basis = (size > 0) ? new double[size] : 0;
   for (int i=0; i < size; ++i)
-    _pBasisFns[i] = pBasisFns[i];
+    _basis[i] = basis[i];
 
-  size = numCorners * numQuadPts * numDims;
-  delete[] _pBasisFnsDeriv;
-  _pBasisFnsDeriv = (size > 0) ? new double[size] : 0;
+  size = numCorners * numQuadPts * spaceDim;
+  delete[] _basisDeriv; _basisDeriv = (size > 0) ? new double[size] : 0;
   for (int i=0; i < size; ++i)
-    _pBasisFnsDeriv[i] = pBasisFnsDeriv[i];
+    _basisDeriv[i] = basisDeriv[i];
 
-  size = numQuadPts * numDims;
-  delete[] _pQuadPts;
-  _pQuadPts = (size > 0) ? new double[size] : 0;
+  size = numQuadPts * cellDim;
+  delete[] _quadPtsRef; _quadPtsRef = (size > 0) ? new double[size] : 0;
   for (int i=0; i < size; ++i)
-    _pQuadPts[i] = pQuadPts[i];
+    _quadPtsRef[i] = quadPtsRef[i];
 
   size = numQuadPts;
-  delete[] _pQuadWts;
-  _pQuadWts = (size > 0) ? new double[size] : 0;
+  delete[] _quadWts; _quadWts = (size > 0) ? new double[size] : 0;
   for (int i=0; i < size; ++i)
-    _pQuadWts[i] = pQuadWts[i];
+    _quadWts[i] = quadWts[i];
 
-  _numDims = numDims;
+  _cellDim = cellDim;
   _numCorners = numCorners;
   _numQuadPts = numQuadPts;
+  _spaceDim = spaceDim;
+
+  // Allocate for Jacobian and its inverse
+  size = numQuadPts*cellDim*spaceDim;
+  delete[] _jacobian; _jacobian = (size > 0) ? new double[size] : 0;
+  delete[] _jacobianInv; _jacobianInv = (size > 0) ? new double[size] : 0;
+
+  // Allocate for Jacobian determinant
+  size = numQuadPts;
+  delete[] _jacobianDet; _jacobianDet = (size > 0) ? new double[size] : 0;
+
+  // Allocate for quad pts
+  size = numQuadPts*spaceDim;
+  delete[] _quadPts; _quadPts = (size > 0) ? new double[size] : 0;
 } // initialize
 
 // End of file 

Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh	2006-09-18 03:30:20 UTC (rev 4566)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh	2006-09-18 03:35:26 UTC (rev 4567)
@@ -10,15 +10,24 @@
 // ======================================================================
 //
 
-// @file pylith/feassemble/Quadrature.hh
+/**
+ * @file pylith/feassemble/Quadrature.hh
+ *
+ * @brief Abstract base class for integrating over finite-elements
+ * using quadrature.
+ *
+ * This object contains the basis functions and their derivatives
+ * evaluated at the quadrature points of the reference element, and
+ * the coordinates and weights of the quadrature points. Given a cell
+ * this object will compute the cell's Jacobian, the determinant of
+ * the Jacobian, the inverse of the Jacobian, and the coordinates in
+ * the domain of the cell's quadrature points. The Jacobian and its
+ * inverse are computed at the quadrature points.
+ *
+ * The memory for the Jacobian and its associated information are
+ * managed locally.
+ */
 
-// @brief Abstract base class for integrating over finite-elements
-// using quadratures.
-//
-// This object holds the basis functions and their derivatives
-// evaluated at the quadrature points, and the coordinates and weights
-// of the quadrature points.
-
 #if !defined(pylith_feassemble_quadrature_hh)
 #define pylith_feassemble_quadrature_hh
 
@@ -32,65 +41,64 @@
 
 class pylith::feassemble::Quadrature
 { // Quadrature
-  
-// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+  friend class TestQuadrature; // unit testing
+
+// PUBLIC METHODS ///////////////////////////////////////////////////////
 public :
 
   /// Constructor
   Quadrature(void);
 
   /// Destructor
-  virtual ~Quadrature(void);
+  virtual
+  ~Quadrature(void);
 
-  /** Compute geometric quantities for a cell.
-   *
-   * @param coordinates Section containing vertex coordinates
-   * @param cell Finite-element cell
-   * @param pV Array of ???
-   * @param pJacobian Jacobian evaluated at quadrature points
-   *   size = numDims*numDims*numQuadPts
-   *   index = iQuadPt*numDims*numDims + iJacobian
-   * @param pJacobianInv Inverse Jacobian evaluated at quadrature points
-   *   quadrature pts
-   *   size = numDims*numDims*numQuadPts
-   *   index = iQuadPt*numDims*numDims + iJacobian
-   * @param jacobianDet Determinant of Jacobian
-   */
-  virtual void compute(const ALE::Obj<ALE::Mesh::section_type>& coordinates,
-		       const ALE::Mesh::point_type& cell,
-		       double* pV,
-		       double* pJacobian,
-		       double* pJacobianInv,
-		       double& jacobianDet) = 0;
+  /// Create a copy of this object.
+  virtual
+  Quadrature* clone(void) const = 0;
 
-  /** Set basis functions and their derivatives and coordinates and
+  /** Set basis functions and their derivatives, and coordinates and
    *  weights of the quadrature points.
    *
-   * @param pBasisFns Array of basis functions evaluated at quadrature pts
-   *   index = iVertex*numDims + iDimension
+   * @param basis Array of basis functions evaluated at quadrature pts
+   *   N0Qp0, N0Qp1, ...
+   *   N1Qp0, N1Qp1, ...
+   *   ...
+   *   size = numCorners * numQuadPts
+   *   index = iBasis*numQuadPts + iQuadPt
    *
-   * @param pBasisFnsDeriv Array of basis function derivaties evaluated 
+   * @param basisDeriv Array of basis function derivaties evaluated 
    *   at quadrature pts
-   *   index = iVertex*numDims + iDimension??
+   *   N0xQp0, N0yQp0, N0zQp0, N0xQp1, N0yQp1, N0zQp1, ... 
+   *   N1xQp0, N1yQp0, N1zQp0, N1xQp1, N1yQp1, N1zQp1, ...
+   *   ...
+   * size = numCorners * numQuadPts * spaceDim
+   * index = iVertex*numQuadPts*spaceDim + iQuadPt*spaceDim + iDim
    *
-   * @param pQuadPts Array of coordinates of quadrature points in 
-   *   reference element
-   *   index = iQuadPt*numDims + iDimension
+   * @param quadPts Array of coordinates of quadrature points in 
+   *   reference cell
+   *   Qp0x, Qp0y, Qp0z
+   *   Qp1x, Qp1y, Qp1z
+   *   size = numQuadPts * numDims
+   *   index = iQuadPts*numDims + iDim
    *
-   * @param pQuadWts Array of weights of quadrature points
+   * @param quadWts Array of weights of quadrature points
+   *   WtQp0, WtQp1, ...
    *   index = iQuadPt
    *
-   * @param numDims Number of dimensions
+   * @param cellDim Number of dimensions in reference cell
    * @param numCorners Number of vertices in a cell
    * @param numQuadPts Number of quadrature points
+   * @param spaceDim Number of dimensions in coordinates of cell vertices
    */
-  void initialize(const double* pBasisFns,
-		  const double* pBasisFnsDeriv,
-		  const double* pQuadPts,
-		  const double* pQuadWts,
-		  const int numDims,
+  void initialize(const double* basis,
+		  const double* basisDeriv,
+		  const double* quadPtsRef,
+		  const double* quadWts,
+		  const int cellDim,
 		  const int numCorners,
-		  const int numQuadPts);
+		  const int numQuadPts,
+		  const int spaceDim);
 
   /** Set tolerance for minimum allowable Jacobian.
    *
@@ -98,98 +106,125 @@
    */
   void jacobianTolerance(const double tolerance);
 
-  /** Set tolerance for minimum allowable Jacobian.
+  /** Get tolerance for minimum allowable Jacobian.
    *
-   * @param tolerance Minimum allowable value for Jacobian
+   * @returns Minimum allowable value for Jacobian
    */
   double jacobianTolerance(void);
 
-  /** Get number of dimensions.
-   *
-   * @returns Number of dimensions
-   */
-  int numDims(void) const;
+// PROTECTED METHODS ////////////////////////////////////////////////////
+protected :
 
-  /** Get number of vertices in a cell.
+  /** Copy constructor.
    *
-   * @returns Number of vertices in a cell
+   * @param q Quadrature to copy
    */
-  int numCorners(void) const;
+  Quadrature(const Quadrature& q);
 
-  /** Get number of quadrature points.
+  /** Compute geometric quantities for a cell.
    *
-   * @param returns Number of quadrature points
+   * @param coordinates Section containing vertex coordinates
+   * @param cell Finite-element cell
    */
-  int numQuadPts(void) const;
+  virtual 
+  void _computeGeometry(const ALE::Obj<ALE::Mesh::section_type>& coordinates,
+			const ALE::Mesh::point_type& cell) = 0;
 
+// PRIVATE METHODS //////////////////////////////////////////////////////
+private :
+
+  const Quadrature& operator=(const Quadrature&); ///< Not implemented
+
 // PROTECTED MEMBERS ////////////////////////////////////////////////////
 protected :
 
-  /** Get basis functions evaluated at quadrature points.
+  double _jacobianTol; ///< Tolerance for minium allowable Jacobian determinant
+  
+  /** Array of basis functions evaluated at the quadrature points.
    *
-   * @returns Array of basis functions evaluated at quadrature points.
+   * N0Qp0, N0Qp1, ...
+   * N1Qp0, N1Qp1, ...
+   *
+   * size = numCorners * numQuadPts
+   * index = iBasis*numQuadPts + iQuadPt
    */
-  const double* basisFns(void) const;
+  double* _basis; ///< Array of basis fns evaluated at quad pts
 
-  /** Get derivatives of basis functions evaluated at quadrature points.
+  /** Array of basis functions evaluated at the quadrature points.
    *
-   * @returns Array of basis functions evaluated at quadrature points.
+   * N0xQp0, N0yQp0, N0zQp0, N0xQp1, N0yQp1, N0zQp1, ... 
+   * N1xQp0, N1yQp0, N1zQp0, N1xQp1, N1yQp1, N1zQp1, ...
+   *
+   * size = numCorners * numQuadPts * spaceDim
+   * index = iVertex*numQuadPts*spaceDim + iQuadPt*spaceDim + iDim
    */
-  const double* basisFnsDeriv(void) const;
+  double* _basisDeriv;
 
-  /** Get coordinates of quadrature points.
+  /** Array of coordinates of quadrature points in reference cell.
    *
-   * @returns Array of coordinates
+   * Reference coordinates: (p,q,r)
+   *
+   * Qp0p, Qp0q, Qp0r
+   * Qp1p, Qp1q, Qp1r
+   *
+   * size = numQuadPts * cellDim
+   * index = iQuadPts*cellDim + iDim
    */
-  const double* quadPts(void) const;
+  double* _quadPtsRef;
 
-  /** Get weights of quadrature points.
+  /** Array of coordinates of quadrature points in cell (NOT reference cell).
    *
-   * @returns Array of weights
+   * Qp0x, Qp0y, Qp0z
+   * Qp1x, Qp1y, Qp1z
+   *
+   * size = numQuadPts * spaceDim
+   * index = iQuadPts*spaceDim + iDim
    */
-  const double* quadWts(void) const;
+  double* _quadPts;
 
-// PRIVATE MEMBERS //////////////////////////////////////////////////////
-private :
-
-  double _jacobianTol; ///< Tolernace for small Jacobian determinant
-  
-  /** Array of basis functions evaluated at the quadrature points.
+  /** Array of weights of quadrature points.
    *
-   * N1Qp1, N1Qp2, ...
-   * N2Qp1, N2Qp2, ...
-   *
-   * size = numCorners * numQuadPts
-   * index = iBasis*numQuadPts + iQuadPt
+   * WtQp0, WtQp1, ...
+   * size = numQuadPts
+   * index = iQuadPt
    */
-  double* _pBasisFns; ///< Array of basis fns evaluated at quad pts
+  double* _quadWts;
 
-  /** Array of basis functions evaluated at the quadrature points.
+  /** Array of Jacobian evaluated at quadrature points.
    *
-   * N1xQp1, N1yQp1, N1zQp1, N1xQp2, N1yQp2, N1zQp2, ... 
-   * N2xQp1, N2yQp1, N2zQp1, N2xQp2, N2yQp2, N2zQp2, ...
+   * Qp0J00, Qp0J01, Qp0J02, ...
+   * Qp1J00, Qp1J01, Qp1J02, ...
+   * ...
    *
-   * size = numCorners * numQuadPts * numDeriv
-   * index = iBasis*numQuadPts*numDeriv + iQuadPt*numDeriv + iDeriv
+   * size = numQuadPts*cellDim*spaceDim
+   * index = iQuadPt*cellDim*spaceDim + iRow*spaceDim + iCol
    */
-  double* _pBasisFnsDeriv;
+  double* _jacobian;
 
-  /** Array of coordinates of quadrature points.
+  /** Array of Jacobian inverses evaluated at quadrature points.
    *
-   * size = numQuadPts * numDims
-   * index = iQuadPts*numDims + iDim
+   * Qp0Jinv00, Qp0Jinv01, Qp0Jinv02, ...
+   * Qp1Jinv00, Qp1Jinv01, Qp1Jinv02, ...
+   * ...
+   *
+   * size = numQuadPts*spaceDim*cellDim
+   * index = iQuadPt*spaceDim*cellDim + iRow*cellDim + iCol
    */
-  double* _pQuadPts; ///< Array of coordinates of quadrature points
+  double* _jacobianInv;
 
-  /** Array of weights of quadrature points.
+  /** Array of determinant of Jacobian evaluated at quadrature points.
    *
-   * WtQp1, WtQp2, ...
+   * JdetQp0, JdetQp1, ...
+   *
+   * size = numQuadPts
+   * index = iQuadPt
    */
-  double* _pQuadWts; ///< Array of weights of quadrature points
+  double* _jacobianDet;
 
-  int _numDims; ///< Number of dimensions
+  int _cellDim; ///< Number of dimensions in reference cell
   int _numCorners; ///< Number of vertices in cell
   int _numQuadPts; ///< Number of quadrature points
+  int _spaceDim; ///< Number of dimensions in coordinates of cell vertices
 
 }; // Quadrature
 

Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.icc	2006-09-18 03:30:20 UTC (rev 4566)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.icc	2006-09-18 03:35:26 UTC (rev 4567)
@@ -14,30 +14,20 @@
 #error "Quadrature.icc must be included only from Quadrature.hh"
 #else
 
-// Get basis functions evaluated at quadrature points.
-const double*
-pylith::feassemble::Quadrature::basisFns(void) const {
-  return _pBasisFns;
+// Get tolerance for minimum allowable Jacobian.
+inline
+double
+pylith::feassemble::Quadrature::jacobianTolerance(void) {
+  return _jacobianTol;
 }
 
-// Get derivatives of basis functions evaluated at quadrature points.
-const double*
-pylith::feassemble::Quadrature::basisFnsDeriv(void) const {
-  return _pBasisFnsDeriv;
+// Set tolerance for minimum allowable Jacobian.
+inline
+void
+pylith::feassemble::Quadrature::jacobianTolerance(const double tolerance) {
+  _jacobianTol = tolerance;
 }
 
-// Get coordinates of quadrature points.
-const double*
-pylith::feassemble::Quadrature::quadPts(void) const {
-  return _pQuadPts;
-}
-
-// Get weights of quadrature points.
-const double*
-pylith::feassemble::Quadrature::quadWts(void) const {
-  return _pQuadWts;
-}
-
 #endif
 
 // End of file

Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.cc	2006-09-18 03:30:20 UTC (rev 4566)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.cc	2006-09-18 03:35:26 UTC (rev 4567)
@@ -10,88 +10,87 @@
 // ======================================================================
 //
 
-#include "ElemGeometry1D.hh" // implementation of class methods
+#include <portinfo>
 
+#include "Quadrature1D.hh" // implementation of class methods
+
+#include <assert.h> // USES assert()
+#include <stdexcept> // USES std::runtime_error
+#include <sstream> // USES std::ostringstream
+
 // ----------------------------------------------------------------------
 // Constructor
-pylith::feassemble::ElemGeometry1D::ElemGeometry1D(void)
+pylith::feassemble::Quadrature1D::Quadrature1D(void)
 { // constructor
 } // constructor
 
 // ----------------------------------------------------------------------
 // Destructor
-pylith::feassemble::ElemGeometry1D::~ElemGeometry1D(void)
+pylith::feassemble::Quadrature1D::~Quadrature1D(void)
 { // destructor
 } // destructor
   
 // ----------------------------------------------------------------------
-/** Compute geometry of element object.
- *
- */
-void compute(const Obj<section_type>& coordinates,
-	     const point_type& cell,
-	     value_type* pV,
-	     value_type* pJacobian,
-	     valye_type* pJacobianInv,
-	     value_type& jacobianDet)
-{ // compute
-  assert(0 != pV);
-  assert(0 != pJacobian);
-  assert(0 != pJacobianInv);
+// Copy constructor.
+pylith::feassemble::Quadrature1D::Quadrature1D(const Quadrature1D& q) :
+  Quadrature(q)
+{ // copy constructor
+} // copy constructor
 
-  const int ndims = 1;
-  const topology_type::patch_type patch  = 0;
-  const value_type* coords = coordinates->restrict(patch, cell);
+// ----------------------------------------------------------------------
+// Compute geometric quantities for a cell.
+void
+pylith::feassemble::Quadrature1D::_computeGeometry(
+		       const ALE::Obj<ALE::Mesh::section_type>& coordinates,
+		       const ALE::Mesh::point_type& cell)
+{ // _computeGeometry
+  assert(1 == _cellDim);
+  assert(1 == _spaceDim);
+  assert(0 != cell);
+  assert(0 != _basisDeriv);
+  assert(0 != _quadPtsRef);
+  assert(0 != _quadPts);
+  assert(0 != _quadWts);
+  assert(0 != _jacobian);
+  assert(0 != _jacobianInv);
 
-  // Matt: What is pV???
-  // Are they the coordinates of the vertices of the cell?
-  pV[0] = coords[0];
+  // Get coordinates of cell's vertices
+  const ALE::Mesh::topology_type::patch_type patch  = 0;
+  const ALE::Mesh::section_type::value_type* vertCoords = 
+    coordinates->restrict(patch, cell);
+  //assert(1 == coordinates.GetFiberDimension(patch, *vertices->begin()));
 
-  
+  // Loop over quadrature points
+  for (int iQuadPt=0; iQuadPt < _numQuadPts; ++iQuadPt) {
 
+    // Compute coordinates of quadrature point in cell
+    // x = sum[i=1,n] (Ni * xi)
+    _quadPts[iQuadPt] = 0.0;
+    for (int iVertex=0; iVertex < _numCorners; ++iVertex)
+      _quadPts[iQuadPt] += 
+	_basis[iVertex*_numQuadPts+iQuadPt]*vertCoords[iVertex];
 
-  value_type* invDet;
+    // Compute Jacobian at quadrature point
+    // J = dx/dp = sum[i=1,n] (dNi/dp * xi)
+    _jacobian[iQuadPt] = 0.0;
+    for (int iVertex=0; iVertex < _numCorners; ++iVertex)
+      _jacobian[iQuadPt] += 
+	_basisDeriv[iVertex*_numQuadPts+iQuadPt] * vertCoords[iVertex];
 
+    // Compute determinant of Jacobian at quadrature point
+    // |J| = j00
+    const double det = _jacobian[iQuadPt];
+    if (det < _jacobianTol) {
+      std::ostringstream msg;
+      msg << "Determinant of Jacobian is below minimum permissible value!\n";
+      throw std::runtime_error(msg.str());
+    } // for
+    _jacobianDet[iQuadPt] = _jacobian[iQuadPt];
 
+    // Compute inverse of Jacobian at quadrature point
+    // Jinv = 1/j00
+    _jacobianInv[iQuadPt] = 1.0/_jacobianDet[iQuadPt];
+  } // for
+} // _computeGeometry
 
-
-
-  
-  for(int d = 0; d < _dim; d++) {
-    for(int f = 0; f < _dim; f++) {
-      J[d*_dim+f] = 0.5*(coords[(f+1)*_dim+d] - coords[0*_dim+d]);
-    }
-  }
-  if (_dim == 1) {
-    detJ = J[0];
-  } else if (_dim == 2) {
-    detJ = J[0]*J[3] - J[1]*J[2];
-  } else if (_dim == 3) {
-    detJ = J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
-      J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
-      J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
-  }
-  if (invJ) {
-    invDet = 1.0/detJ;
-    if (_dim == 2) {
-      invJ[0] =  invDet*J[3];
-      invJ[1] = -invDet*J[1];
-      invJ[2] = -invDet*J[2];
-      invJ[3] =  invDet*J[0];
-    } else if (_dim == 3) {
-      // FIX: This may be wrong
-      invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]);
-      invJ[0*3+1] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]);
-      invJ[0*3+2] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
-      invJ[1*3+0] = invDet*(J[0*3+1]*J[2*3+2] - J[0*3+2]*J[2*3+1]);
-      invJ[1*3+1] = invDet*(J[0*3+2]*J[2*3+0] - J[0*3+0]*J[2*3+2]);
-      invJ[1*3+2] = invDet*(J[0*3+0]*J[2*3+1] - J[0*3+1]*J[2*3+0]);
-      invJ[2*3+0] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]);
-      invJ[2*3+1] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]);
-      invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]);
-    }
-  }
-  
-} // compute
-
 // End of file 

Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.hh	2006-09-18 03:30:20 UTC (rev 4566)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.hh	2006-09-18 03:35:26 UTC (rev 4567)
@@ -10,6 +10,12 @@
 // ======================================================================
 //
 
+/**
+ * @file pylith/feassemble/Quadrature1D.hh
+ *
+ * @brief Quadrature for 1-D finite-elements.
+ */
+
 #if !defined(pylith_feassemble_quadrature1d_hh)
 #define pylith_feassemble_quadrature1d_hh
 
@@ -18,12 +24,14 @@
 namespace pylith {
   namespace feassemble {
     class Quadrature1D;
+    class TestQuadrature1D;
   } // feassemble
 } // pylith
 
-class pylith::feassemble::Quadrature1D
+class pylith::feassemble::Quadrature1D : public Quadrature
 { // Quadrature1D
-  
+  friend class TestQuadrature1D; // unit testing
+
 // PUBLIC MEMBERS ///////////////////////////////////////////////////////
 public :
 
@@ -31,20 +39,38 @@
   Quadrature1D(void);
 
   /// Destructor
-  ~Quadrature1D(void);
+  virtual ~Quadrature1D(void);
 
-  /** Compute geometry of element object.
+  /// Create a copy of this object.
+  virtual
+  Quadrature* clone(void) const;
+
+// PROTECTED METHODS ////////////////////////////////////////////////////
+protected :
+
+  /** Copy constructor.
    *
+   * @param q Quadrature to copy
    */
-  void compute(const ALE::Obj<ALE::Mesh::section_type>& coordinates,
-	       const ALE::Mesh::point_type& cell,
-	       value_type* pV,
-	       value_type* pJacobian,
-	       valye_type* pJacobianInv,
-	       value_type* pJacobianDet);
+  Quadrature1D(const Quadrature1D& q);
 
+  /** Compute geometric quantities for a cell.
+   *
+   * @param coordinates Section containing vertex coordinates
+   * @param cell Finite-element cell
+   */
+  void _computeGeometry(const ALE::Obj<ALE::Mesh::section_type>& coordinates,
+			const ALE::Mesh::point_type& cell);
+
+// PRIVATE METHODS //////////////////////////////////////////////////////
+private :
+
+  const Quadrature& operator=(const Quadrature&); ///< Not implemented
+
 }; // Quadrature1D
 
+#include "Quadrature1D.icc" // inline methods
+
 #endif // pylith_feassemble_quadrature1d_hh
 
 // End of file 

Added: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.icc	2006-09-18 03:30:20 UTC (rev 4566)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.icc	2006-09-18 03:35:26 UTC (rev 4567)
@@ -0,0 +1,26 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+//                           Brad T. Aagaard
+//                        U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#if !defined(pylith_feassemble_quadrature1d_hh)
+#error "Quadrature1D.icc must be included only from Quadrature1D.hh"
+#else
+
+// Create a copy of this object.
+inline
+pylith::feassemble::Quadrature*
+pylith::feassemble::Quadrature1D::clone(void) const {
+  return new Quadrature1D(*this);
+} // clone
+
+#endif
+
+// End of file

Modified: short/3D/PyLith/trunk/libsrc/meshio/MeshIO.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/MeshIO.cc	2006-09-18 03:30:20 UTC (rev 4566)
+++ short/3D/PyLith/trunk/libsrc/meshio/MeshIO.cc	2006-09-18 03:35:26 UTC (rev 4567)
@@ -10,6 +10,8 @@
 // ======================================================================
 //
 
+#include <portinfo>
+
 #include "MeshIO.hh" // implementation of class methods
 
 // ----------------------------------------------------------------------

Modified: short/3D/PyLith/trunk/libsrc/meshio/MeshIOAscii.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/MeshIOAscii.cc	2006-09-18 03:30:20 UTC (rev 4566)
+++ short/3D/PyLith/trunk/libsrc/meshio/MeshIOAscii.cc	2006-09-18 03:35:26 UTC (rev 4567)
@@ -10,7 +10,8 @@
 // ======================================================================
 //
 
-#include "MeshIO.hh" // MeshIOAscii ISA MeshIO
+#include <portinfo>
+
 #include "MeshIOAscii.hh" // implementation of class methods
 
 #include <fstream> // USES std::ifstream, std::ofstream

Added: short/3D/PyLith/trunk/unittests/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/Makefile.am	2006-09-18 03:30:20 UTC (rev 4566)
+++ short/3D/PyLith/trunk/unittests/Makefile.am	2006-09-18 03:35:26 UTC (rev 4567)
@@ -0,0 +1,16 @@
+# -*- Makefile -*-
+#
+# ----------------------------------------------------------------------
+#
+#                           Brad T. Aagaard
+#                        U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+SUBDIRS = \
+	libtests
+
+# End of file 

Added: short/3D/PyLith/trunk/unittests/libtests/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/Makefile.am	2006-09-18 03:30:20 UTC (rev 4566)
+++ short/3D/PyLith/trunk/unittests/libtests/Makefile.am	2006-09-18 03:35:26 UTC (rev 4567)
@@ -0,0 +1,16 @@
+# -*- Makefile -*-
+#
+# ----------------------------------------------------------------------
+#
+#                           Brad T. Aagaard
+#                        U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+SUBDIRS = \
+	feassemble
+
+# End of file 

Added: short/3D/PyLith/trunk/unittests/libtests/feassemble/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/Makefile.am	2006-09-18 03:30:20 UTC (rev 4566)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/Makefile.am	2006-09-18 03:35:26 UTC (rev 4567)
@@ -0,0 +1,40 @@
+# -*- Makefile -*-
+#
+# ----------------------------------------------------------------------
+#
+#                           Brad T. Aagaard
+#                        U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+subpackage = feassemble
+include $(top_srcdir)/subpackage.am
+
+TESTS = testfeassemble
+
+check_PROGRAMS = testfeassemble
+
+testfeassemble_SOURCES = \
+	TestQuadrature.cc \
+	TestQuadrature1D.cc \
+	test_feassemble.cc
+
+noinst_HEADERS = \
+	TestQuadrature.hh \
+	TestQuadrature1D.hh
+
+testfeassemble_LDFLAGS = $(PETSC_LIB)
+
+INCLUDES += $(PETSC_INCLUDE)
+
+testfeassemble_LDADD = \
+	-lcppunit -ldl \
+	$(top_builddir)/libsrc/feassemble/libpylithfeassemble.la
+
+# version
+# $Id$
+
+# End of file 

Added: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.cc	2006-09-18 03:30:20 UTC (rev 4566)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.cc	2006-09-18 03:35:26 UTC (rev 4567)
@@ -0,0 +1,163 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+//                           Brad T. Aagaard
+//                        U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "TestQuadrature.hh" // Implementation of class methods
+
+#include "pylith/feassemble/Quadrature1D.hh" // USES Quadrature1D
+
+#include <sstream> // USES std::stringstream
+
+// ----------------------------------------------------------------------
+CPPUNIT_TEST_SUITE_REGISTRATION( pylith::feassemble::TestQuadrature );
+
+// ----------------------------------------------------------------------
+// Test clone
+void
+pylith::feassemble::TestQuadrature::testClone(void)
+{ // testClone
+  Quadrature1D qOrig;
+  
+  const double jacobianTol = 1.0;
+
+  const int cellDim = 1;
+  const int numCorners = 2;
+  const int numQuadPts = 1;
+  const int spaceDim = 1;
+  const double basis[] = { 0.5, 0.5 };
+  const double basisDeriv[] = { -1.0, 1.0 };
+  const double quadPtsRef[] = { 0.5 };
+  const double quadWts[] = { 1.0 };
+
+  // Semi-random values manually set to check cloning
+  const double quadPts[] = { 2.0 };
+  const double jacobian[] = { 4.0 };
+  const double jacobianInv[] = { 0.25 };
+  const double jacobianDet[] = { 5.0 };
+
+  qOrig.jacobianTolerance(jacobianTol);
+  qOrig.initialize(basis, basisDeriv, quadPtsRef, quadWts,
+		   cellDim, numCorners, numQuadPts, spaceDim);
+  int size = 1;
+  memcpy(qOrig._quadPts, quadPts, size*sizeof(double));
+  memcpy(qOrig._jacobian, jacobian, size*sizeof(double));
+  memcpy(qOrig._jacobianInv, jacobianInv, size*sizeof(double));
+  memcpy(qOrig._jacobianDet, jacobianDet, size*sizeof(double));
+  
+  const Quadrature* qCopy = qOrig.clone();
+  CPPUNIT_ASSERT(0 != qCopy);
+
+  // Make sure values match those passed
+  CPPUNIT_ASSERT_EQUAL(jacobianTol, qCopy->_jacobianTol);
+
+  CPPUNIT_ASSERT_EQUAL(cellDim, qCopy->_cellDim);
+  CPPUNIT_ASSERT_EQUAL(numCorners, qCopy->_numCorners);
+  CPPUNIT_ASSERT_EQUAL(numQuadPts, qCopy->_numQuadPts);
+  CPPUNIT_ASSERT_EQUAL(spaceDim, qCopy->_spaceDim);
+
+  CPPUNIT_ASSERT(0 != qCopy->_basis);
+  size = numCorners * numQuadPts;
+  for (int i=0; i < size; ++i)
+    CPPUNIT_ASSERT_EQUAL(basis[i], qCopy->_basis[i]);
+
+  CPPUNIT_ASSERT(0 != qCopy->_basisDeriv);
+  size = numCorners * numQuadPts * spaceDim;
+  for (int i=0; i < size; ++i)
+    CPPUNIT_ASSERT_EQUAL(basisDeriv[i], qCopy->_basisDeriv[i]);
+
+  CPPUNIT_ASSERT(0 != qCopy->_quadPtsRef);
+  size = numQuadPts * cellDim;
+  for (int i=0; i < size; ++i)
+    CPPUNIT_ASSERT_EQUAL(quadPtsRef[i], qCopy->_quadPtsRef[i]);
+
+  CPPUNIT_ASSERT(0 != qCopy->_quadWts);
+  size = numQuadPts;
+  for (int i=0; i < size; ++i)
+    CPPUNIT_ASSERT_EQUAL(quadWts[i], qCopy->_quadWts[i]);
+
+  size = 1;
+
+  CPPUNIT_ASSERT(0 != qCopy->_quadPts);
+  for (int i=0; i < size; ++i)
+    CPPUNIT_ASSERT_EQUAL(quadPts[i], qCopy->_quadPts[i]);
+  CPPUNIT_ASSERT(0 != qCopy->_jacobian);
+  for (int i=0; i < size; ++i)
+    CPPUNIT_ASSERT_EQUAL(jacobian[i], qCopy->_jacobian[i]);
+  CPPUNIT_ASSERT(0 != qCopy->_jacobianInv);
+  for (int i=0; i < size; ++i)
+    CPPUNIT_ASSERT_EQUAL(jacobianInv[i], qCopy->_jacobianInv[i]);
+  CPPUNIT_ASSERT(0 != qCopy->_jacobianDet);
+  for (int i=0; i < size; ++i)
+    CPPUNIT_ASSERT_EQUAL(jacobianDet[i], qCopy->_jacobianDet[i]);
+} // testCopy
+
+// ----------------------------------------------------------------------
+// Test jacobianTolerance()
+void
+pylith::feassemble::TestQuadrature::testJacobianTol(void)
+{ // testJacobianTol
+  Quadrature1D q;
+  const double tolerance = 1.0;
+  q.jacobianTolerance(tolerance);
+  CPPUNIT_ASSERT_EQUAL(tolerance, q._jacobianTol);
+} // testJacobianTol
+
+// ----------------------------------------------------------------------
+// Test initialize()
+void
+pylith::feassemble::TestQuadrature::testInitialize(void)
+{ // initialize
+  Quadrature1D q;
+  
+  const int cellDim = 1;
+  const int numCorners = 2;
+  const int numQuadPts = 1;
+  const int spaceDim = 1;
+  const double basis[] = { 0.5, 0.5 };
+  const double basisDeriv[] = { -1.0, 1.0 };
+  const double quadPtsRef[] = { 0.5 };
+  const double quadWts[] = { 1.0 };
+  const double jacobianTol = 1.0;
+
+  q.initialize(basis, basisDeriv, quadPtsRef, quadWts,
+	       cellDim, numCorners, numQuadPts, spaceDim);
+  
+  CPPUNIT_ASSERT_EQUAL(cellDim, q._cellDim);
+  CPPUNIT_ASSERT_EQUAL(numCorners, q._numCorners);
+  CPPUNIT_ASSERT_EQUAL(numQuadPts, q._numQuadPts);
+  CPPUNIT_ASSERT_EQUAL(spaceDim, q._spaceDim);
+
+  int size = numCorners * numQuadPts;
+  for (int i=0; i < size; ++i)
+    CPPUNIT_ASSERT_EQUAL(basis[i], q._basis[i]);
+
+  size = numCorners * numQuadPts * spaceDim;
+  for (int i=0; i < size; ++i)
+    CPPUNIT_ASSERT_EQUAL(basisDeriv[i], q._basisDeriv[i]);
+
+  size = numQuadPts * cellDim;
+  for (int i=0; i < size; ++i)
+    CPPUNIT_ASSERT_EQUAL(quadPtsRef[i], q._quadPtsRef[i]);
+
+  size = numQuadPts;
+  for (int i=0; i < size; ++i)
+    CPPUNIT_ASSERT_EQUAL(quadWts[i], q._quadWts[i]);
+
+  // Make sure Jacobian stuff has been allocated
+  CPPUNIT_ASSERT(0 != q._jacobian);
+  CPPUNIT_ASSERT(0 != q._jacobianInv);
+  CPPUNIT_ASSERT(0 != q._jacobianDet);
+  CPPUNIT_ASSERT(0 != q._quadPts);
+} // initialize
+
+// End of file 

Added: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.hh	2006-09-18 03:30:20 UTC (rev 4566)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.hh	2006-09-18 03:35:26 UTC (rev 4567)
@@ -0,0 +1,60 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+//                           Brad T. Aagaard
+//                        U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/**
+ * @file unittests/libtests/feassemble/TestQuadrature.hh
+ *
+ * @brief C++ TestQuadrature object
+ *
+ * C++ unit testing for Quadrature.
+ */
+
+#if !defined(pylith_feassemble_testquadrature_hh)
+#define pylith_feassemble_testquadrature_hh
+
+#include <cppunit/extensions/HelperMacros.h>
+
+/// Namespace for spatialdata package
+namespace pylith {
+  namespace feassemble {
+    class TestQuadrature;
+  } // feassemble
+} // pylith
+
+/// C++ unit testing for Quadrature
+class pylith::feassemble::TestQuadrature : public CppUnit::TestFixture
+{ // class TestQuadrature
+
+  // CPPUNIT TEST SUITE /////////////////////////////////////////////////
+  CPPUNIT_TEST_SUITE( TestQuadrature );
+  CPPUNIT_TEST( testClone );
+  CPPUNIT_TEST( testJacobianTol );
+  CPPUNIT_TEST( testInitialize );
+  CPPUNIT_TEST_SUITE_END();
+
+  // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+  /// Test clone()
+  void testClone(void);
+
+  /// Test jacobianTolerance()
+  void testJacobianTol(void);
+
+  /// Test initialize()
+  void testInitialize(void);
+
+}; // class TestQuadrature
+
+#endif // pylith_feassemble_testquadrature_hh
+
+// End of file 

Added: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature1D.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature1D.cc	2006-09-18 03:30:20 UTC (rev 4566)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature1D.cc	2006-09-18 03:35:26 UTC (rev 4567)
@@ -0,0 +1,48 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+//                           Brad T. Aagaard
+//                        U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "TestQuadrature1D.hh" // Implementation of class methods
+
+#include "pylith/feassemble/Quadrature1D.hh" // USES Quadrature1D
+
+#include <sstream> // USES std::stringstream
+
+// ----------------------------------------------------------------------
+CPPUNIT_TEST_SUITE_REGISTRATION( pylith::feassemble::TestQuadrature1D );
+
+// ----------------------------------------------------------------------
+// Test constructor
+void
+pylith::feassemble::TestQuadrature1D::testConstructor(void)
+{ // testConstructor
+  Quadrature1D quadrature;
+} // testConstructor
+
+// ----------------------------------------------------------------------
+// Test computeGeometry() w/linear basis fns
+void
+pylith::feassemble::TestQuadrature1D::testLinear(void)
+{ // testLinear
+  CPPUNIT_ASSERT(false);
+} // testLinear
+
+// ----------------------------------------------------------------------
+// Test computeGeometry() w/quadratic basis fns
+void
+pylith::feassemble::TestQuadrature1D::testQuadratic(void)
+{ // testQuadratic
+  CPPUNIT_ASSERT(false);
+} // testQuadratic
+
+// End of file 

Added: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature1D.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature1D.hh	2006-09-18 03:30:20 UTC (rev 4566)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature1D.hh	2006-09-18 03:35:26 UTC (rev 4567)
@@ -0,0 +1,60 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+//                           Brad T. Aagaard
+//                        U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/**
+ * @file unittests/libtests/feassemble/TestQuadrature1D.hh
+ *
+ * @brief C++ TestQuadrature1D object
+ *
+ * C++ unit testing for Quadrature1D.
+ */
+
+#if !defined(pylith_feassemble_testquadrature1d_hh)
+#define pylith_feassemble_testquadrature1d_hh
+
+#include "TestQuadrature.hh"
+
+/// Namespace for spatialdata package
+namespace pylith {
+  namespace feassemble {
+    class TestQuadrature1D;
+  } // feassemble
+} // pylith
+
+/// C++ unit testing for Quadrature1D
+class pylith::feassemble::TestQuadrature1D : public TestQuadrature
+{ // class TestQuadrature1D
+
+  // CPPUNIT TEST SUITE /////////////////////////////////////////////////
+  CPPUNIT_TEST_SUITE( TestQuadrature1D );
+  CPPUNIT_TEST( testConstructor );
+  CPPUNIT_TEST( testLinear );
+  CPPUNIT_TEST( testQuadratic );
+  CPPUNIT_TEST_SUITE_END();
+
+  // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+  /// Test constructor
+  void testConstructor(void);
+
+  /// Test initialize() & computeGeometry() w/linear basis fns
+  void testLinear(void);
+
+  /// Test initialize() & computeGeometry() w/quadratic basis fns
+  void testQuadratic(void);
+
+}; // class TestQuadrature
+
+#endif // pylith_feassemble_testquadrature1d_hh
+
+// End of file 

Added: short/3D/PyLith/trunk/unittests/libtests/feassemble/test_feassemble.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/test_feassemble.cc	2006-09-18 03:30:20 UTC (rev 4566)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/test_feassemble.cc	2006-09-18 03:35:26 UTC (rev 4567)
@@ -0,0 +1,49 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+//                           Brad T. Aagaard
+//                        U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <cppunit/extensions/TestFactoryRegistry.h>
+
+#include <cppunit/BriefTestProgressListener.h>
+#include <cppunit/extensions/TestFactoryRegistry.h>
+#include <cppunit/TestResult.h>
+#include <cppunit/TestResultCollector.h>
+#include <cppunit/TestRunner.h>
+#include <cppunit/TextOutputter.h>
+
+int
+main(int argc,
+     char* argv[])
+{ // main
+  // Create event manager and test controller
+  CppUnit::TestResult controller;
+
+  // Add listener to collect test results
+  CppUnit::TestResultCollector result;
+  controller.addListener(&result);
+
+  // Add listener to show progress as tests run
+  CppUnit::BriefTestProgressListener progress;
+  controller.addListener(&progress);
+
+  // Add top suite to test runner
+  CppUnit::TestRunner runner;
+  runner.addTest(CppUnit::TestFactoryRegistry::getRegistry().makeTest());
+  runner.run(controller);
+
+  // Print tests
+  CppUnit::TextOutputter outputter(&result, std::cerr);
+  outputter.write();
+
+  return (result.wasSuccessful() ? 0 : 1);
+} // main
+
+// End of file



More information about the cig-commits mailing list