[cig-commits] r14027 - in short/3D/PyLith/branches/pylith-swig/libsrc: . feassemble

brad at geodynamics.org brad at geodynamics.org
Fri Feb 6 11:16:54 PST 2009


Author: brad
Date: 2009-02-06 11:16:53 -0800 (Fri, 06 Feb 2009)
New Revision: 14027

Added:
   short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/QuadratureBase.cc
   short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/QuadratureBase.hh
   short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/QuadratureBase.icc
Modified:
   short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am
   short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Makefile.am
   short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature.hh
   short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/feassemblefwd.hh
Log:
Worked on quadrature.

Modified: short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am	2009-02-06 16:26:28 UTC (rev 14026)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am	2009-02-06 19:16:53 UTC (rev 14027)
@@ -39,6 +39,7 @@
 	feassemble/GeometryQuad2D.cc \
 	feassemble/GeometryQuad3D.cc \
 	feassemble/GeometryHex3D.cc \
+	feassemble/QuadratureBase.cc \
 	meshio/BinaryIO.cc \
 	meshio/GMVFile.cc \
 	meshio/GMVFileAscii.cc \

Modified: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Makefile.am	2009-02-06 16:26:28 UTC (rev 14026)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Makefile.am	2009-02-06 19:16:53 UTC (rev 14027)
@@ -35,6 +35,8 @@
 	GeometryQuad2D.hh \
 	GeometryQuad3D.hh \
 	GeometryHex3D.hh \
+	QuadratureBase.hh \
+	QuadratureBase.icc \
 	Quadrature.hh \
 	Quadrature.icc \
 	Quadrature0D.hh \

Modified: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature.hh	2009-02-06 16:26:28 UTC (rev 14026)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Quadrature.hh	2009-02-06 19:16:53 UTC (rev 14027)
@@ -16,32 +16,28 @@
  * @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.
+ * This object contains the informatio needed to perform numerical
+ * quadrature over a finite-element cell. It inherits quadrature
+ * information over the reference cell from the QuadratureBase object.
+
+ * 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.
  */
 
 #if !defined(pylith_feassemble_quadrature_hh)
 #define pylith_feassemble_quadrature_hh
 
-#include "pylith/utils/array.hh" // HASA double_array
+// Include directives ---------------------------------------------------
+#include "QuadratureBase.hh" // ISA QuadratureBase
 
-namespace pylith {
-  namespace feassemble {
-    class Quadrature;
+#include "pylith/topology/topologyfwd.hh" // forward declarations
 
-    class CellGeometry; // HOLDSA CellGeometry
-    class TestQuadrature;
-  } // feassemble
-} // pylith
+#include "pylith/utils/array.hh" // HASA double_array
 
+// Quadrature -----------------------------------------------------------
+template<typename mesh_type>
 class pylith::feassemble::Quadrature
 { // Quadrature
   friend class TestQuadrature; // unit testing
@@ -60,74 +56,6 @@
   virtual
   Quadrature* clone(void) const = 0;
 
-  /** Set basis functions and their derivatives, and coordinates and
-   *  weights of the quadrature points.
-   *
-   * @param basis Array of basis functions evaluated at quadrature pts
-   *   N0Qp0, N1Qp0, ...
-   *   N0Qp1, N1Qp1, ...
-   *   ...
-   *   size = numQuadPts * numBasis
-   *   index = iQuadPt*numBasis + iBasis
-   *
-   * @param basisDerivRef Array of basis function derivaties evaluated at
-   * quadrature pts, where derivatives are with respect to cell's
-   * local coordinates.
-   *   N0pQp0, N0qQp0, N0rQp0, N1pQp0, N1qQp0, N1rQp0, ... 
-   *   N0pQp1, N0qQp1, N0rQp1, N1pQp1, N1qQp1, N1rQp1, ...
-   *   ...
-   *   size = numQuadPts * numBasis * cellDim
-   *   index = iQuadPt*numBasis*cellDim + iBasis*cellDim + iDim
-   *
-   * @param quadPts Array of coordinates of quadrature points in 
-   *   reference cell
-   *   Qp0p, Qp0q, Qp0r
-   *   Qp1p, Qp1q, Qp1r
-   *   size = numQuadPts * numDims
-   *   index = iQuadPt*numDims + iDim
-   *
-   * @param quadWts Array of weights of quadrature points
-   *   WtQp0, WtQp1, ...
-   *   index = iQuadPt
-   *
-   * @param cellDim Number of dimensions in reference cell
-   * @param numBasis Number of basis functions for a cell
-   * @param numQuadPts Number of quadrature points
-   * @param spaceDim Number of dimensions in coordinates of cell vertices
-   */
-  void initialize(const double* basis,
-		  const double* basisDerivRef,
-		  const double* quadPtsRef,
-		  const double* quadWts,
-		  const int cellDim,
-		  const int numBasis,
-		  const int numQuadPts,
-		  const int spaceDim);
-
-  /** Set geometry associated with reference cell.
-   *
-   * @param geometry Geometry of reference cell.
-   */
-  void refGeometry(CellGeometry* const geometry);
-
-  /** Get geometry associated with reference cell.
-   *
-   * @returns Geometry of reference cell.
-   */
-  const CellGeometry& refGeometry(void) const;
-
-  /** Set minimum allowable determinant of Jacobian.
-   *
-   * @param tolerance Minimum allowable value for Jacobian
-   */
-  void minJacobian(const double min);
-
-  /** Get minimum allowable determinant of Jacobian.
-   *
-   * @returns Minimum allowable value for Jacobian
-   */
-  double minJacobian(void) const;
-
   /** Set flag for checking ill-conditioning.
    *
    * @param flag True to check for ill-conditioning, false otherwise.
@@ -140,30 +68,12 @@
    */
   bool checkConditioning(void) const;
 
-  /** Get coordinates of quadrature points in reference cell.
-   *
-   * @returns Array of coordinates of quadrature points in reference cell.
-   */
-  const double_array& quadPtsRef(void) const;
-
   /** Get coordinates of quadrature points in cell (NOT reference cell).
    *
    * @returns Array of coordinates of quadrature points in cell
    */
   const double_array& quadPts(void) const;
 
-  /** Get weights of quadrature points.
-   *
-   * @returns Weights of quadrature points
-   */
-  const double_array& quadWts(void) const;
-
-  /** Get basis fns evaluated at quadrature points.
-   *
-   * @returns Array of basis fns evaluated at quadrature points
-   */
-  const double_array& basis(void) const;
-
   /** Get derivatives of basis fns evaluated at quadrature points.
    *
    * @returns Array of derivatives of basis fns evaluated at
@@ -183,78 +93,24 @@
    */
   const double_array& jacobianDet(void) const;
 
-  /** Get Jacobian inverses evaluated at quadrature points.
-   *
-   * @returns Array of Jacobian inverses evaluated at quadrature points.
-   */
-  const double_array& jacobianInv(void) const;
-
-  /** Get number of dimensions in reference cell.
-   *
-   * @returns Number of dimensions in reference cell
-   */
-  int cellDim(void) const;
-
-  /** Get number of basis functions for cell.
-   *
-   * @returns Number of basis functions for cell
-   */
-  int numBasis(void) const;
-
-  /** Get number of quadrature points.
-   *
-   * @returns Number of quadrature points
-   */
-  int numQuadPts(void) const;
-
-  /** Get number of dimensions in coordinates of cell vertices.
-   *
-   * @returns Number of dimensions in coordinates of cell vertices
-   */
-  int spaceDim(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
-   */
-  virtual 
-  void computeGeometry(const double* vertCoords,
-                       const int coordDim,
-                       const SieveMesh::point_type& cell) = 0;
-
-  template<typename mesh_type>
-  void computeGeometry(const ALE::Obj<mesh_type>& mesh,
-                       const ALE::Obj<mesh_type::real_section_type>& coordinates,
-                       const mesh_type::point_type& cell) {
-    computeGeometry(mesh->restrictClosure(coordinates, cell),
-                    coordinates->getFiberDimension(*mesh->depthStratum(0)->begin()),
-                    cell);
-  };
-
-  /** Reset the precomputation structures. */
+  /// Reset the precomputation structures.
   void resetPrecomputation(void);
 
   /** Precompute geometric quantities for each cell.
    *
    * @param mesh Finite-element mesh
-   * @param coordinates Section containing vertex coordinates
+   * @param cells Finite-element cells for geometry.
    */
-  void precomputeGeometry(const topology::Mesh& mesh,
-                          const ALE::Obj<MeshRealSection>& coordinates,
-                          const ALE::Obj<SieveMesh::label_sequence>& cells);
+  void precomputeGeometry(const mesh_type& mesh,
+             const ALE::Obj<mesh_type::SieveMesh::label_sequence>& cells);
 
   /** Retrieve precomputed geometric quantities for a cell.
    *
    * @param mesh Finite-element mesh
-   * @param coordinates Section containing vertex coordinates
    * @param cell Finite-element cell
    */
-  void retrieveGeometry(const ALE::Obj<SieveMesh>& mesh,
-                        const ALE::Obj<MeshRealSection>& coordinates,
-                        const SieveMesh::point_type& cell,
-                        const int c);
+  void retrieveGeometry(const ALE::Obj<mesh_type::SieveMesh>& mesh,
+                        const mesh_type::SieveMesh::point_type& cell);
 
 // PROTECTED METHODS ////////////////////////////////////////////////////
 protected :
@@ -271,141 +127,57 @@
    * @param cell Label of finite-element cell
    */
   void _checkJacobianDet(const double det,
-			 const int cell) const;
+			 const mesh_type::SieveMesh::point_type& cell) const;
 
   /// Set entries in geometry arrays to zero.
   void _resetGeometry(void);
 
-// PRIVATE METHODS //////////////////////////////////////////////////////
-private :
+  /** Compute geometric quantities for a cell at quadrature points.
+   *
+   * @param coordsVertices Coordinates of vertices of cell.
+   * @param spaceDim Number of coordinates per vertex.
+   * @param cell Finite-element cell
+   */
+  virtual 
+  void _computeGeometry(const double* coordsVertices,
+			const int spaceDim,
+			const mesh_type::SieveMesh::point_type& cell) = 0;
 
-  const Quadrature& operator=(const Quadrature&); ///< Not implemented
-
 // PROTECTED MEMBERS ////////////////////////////////////////////////////
 protected :
 
-  double _minJacobian; ///< Minium allowable Jacobian determinant
-  
-  /** Array of coordinates of quadrature points in reference cell.
-   *
-   * Reference coordinates: (p,q,r)
-   *
-   * Qp0p, Qp0q, Qp0r
-   * Qp1p, Qp1q, Qp1r
-   *
-   * size = numQuadPts * cellDim
-   * index = iQuadPts*cellDim + iDim
+  /** Fields and visitors for precomputing geometry information for
+   * cells associated with this quadrature.
    */
-  double_array _quadPtsRef;
-
-  /** Array of coordinates of quadrature points in cell (global coordinates).
-   *
-   * Qp0x, Qp0y, Qp0z
-   * Qp1x, Qp1y, Qp1z
-   *
-   * size = numQuadPts * spaceDim
-   * index = iQuadPts*spaceDim + iDim
+  /** :OPTIMIZATION:
+   * Need to check speed of retrieving geometry using Fields and Visitors.
+   * Can switch to Sieve sections and visitors if too slow.
    */
-  double_array _quadPts;
+  Field<mesh_type>* _quadPtsField; ///< Coordinates of quad pts.
+  Visitor<Field<mesh_type> >* _quadPtsVisitor; ///< Visitor for quad pts.
 
-  /** Array of weights of quadrature points.
-   *
-   * WtQp0, WtQp1, ...
-   * size = numQuadPts
-   * index = iQuadPt
-   */
-  double_array _quadWts;
+  Field<mesh_type>* _jacobianField; ///< Jacobian at quad pts.
+  Visitor<Field<mesh_type> >* _jacobianVisitor; ///< Visitor for Jacobian.
 
-  /** Array of basis functions evaluated at the quadrature points.
-   *
-   * N0Qp0, N1Qp0, ...
-   * N0Qp1, N1Qp1, ...
-   *
-   * size = numQuadPts * numBasis
-   * index = iQuadPt*numBasis + iBasis
-   */
-  double_array _basis;
+  Field<mesh_type>* _jacobianDetField; ///< |J| at quad pts.
+  Visitor<Field<mesh_type> >* _jacobianDetVisitor; ///< Visitor for |J|.
 
-  /** Array of basis function derivatives evaluated at the quadrature
-   * points, where derivatives are with respect to cell's local
-   * coordinates.
-   *
-   * N0pQp0, N0qQp0, N0rQp0, N1pQp0, N1qQp0, N1rQp0, ... 
-   * N0pQp1, N0qQp1, N0rQp1, N1pQp1, N1qQp1, N1rQp1, ...
-   *
-   * size = numQuadPts * numBasis * cellDim
-   * index = iQuadPt*numBasis*cellDim + iBasis*cellDim + iDim
-   */
-  double_array _basisDerivRef;
+  Field<mesh_type>* _basisDerivField; ///< Deriv. of basis fns at quad pts.
+  Visitor<Field<mesh_type> >* _basisDerivVisitor; ///< Visitor for derivatives.
 
-  /** Array of basis function derivatives evaluated at the quadrature
-   * points, where derivatives are with respect to global coordinates.
-   *
-   * N0xQp0, N0yQp0, N0zQp0, N1xQp0, N1yQp0, N1zQp0, ... 
-   * N0xQp1, N0yQp1, N0zQp1, N1xQp1, N1yQp1, N1zQp1, ...
-   *
-   * size = numQuadPts * numBasis * spaceDim
-   * index = iQuadPt*numBasis*spaceDim + iBasis*spaceDim + iDim
-   */
-  double_array _basisDeriv;
+  bool _checkConditioning; ///< True if checking for ill-conditioning.
 
-  /** Array of Jacobians evaluated at quadrature points.
-   *
-   * Qp0J00, Qp0J01, Qp0J02, ...
-   * Qp1J00, Qp1J01, Qp1J02, ...
-   * ...
-   *
-   * size = numQuadPts*cellDim*spaceDim
-   * index = iQuadPt*cellDim*spaceDim + iRow*spaceDim + iCol
-   */
-  double_array _jacobian;
+// NOT IMPLEMENTED //////////////////////////////////////////////////////
+private :
 
-  /** Array of determinant of Jacobian evaluated at quadrature points.
-   *
-   * JdetQp0, JdetQp1, ...
-   *
-   * size = numQuadPts
-   * index = iQuadPt
-   */
-  double_array _jacobianDet;
+  const Quadrature& operator=(const Quadrature&); ///< Not implemented
 
-  /** Array of Jacobian inverses evaluated at quadrature points.
-   *
-   * Qp0Jinv00, Qp0Jinv01, Qp0Jinv02, ...
-   * Qp1Jinv00, Qp1Jinv01, Qp1Jinv02, ...
-   * ...
-   *
-   * size = numQuadPts*spaceDim*cellDim
-   * index = iQuadPt*spaceDim*cellDim + iRow*cellDim + iCol
-   */
-  double_array _jacobianInv;
-
-  int _cellDim; ///< Number of dimensions in reference cell
-  int _numBasis; ///< Number of basis functions (and vertices) for cell
-  int _numQuadPts; ///< Number of quadrature points
-  int _spaceDim; ///< Number of dimensions in coordinates of cell vertices
-
-  CellGeometry* _geometry; ///< Geometry of reference cell
-
-  /* Precomputation sections */
-  int _qTag, _jTag, _jDTag, _jITag, _bTag;
-  ALE::Obj<MeshRealSection> _quadPtsPre;
-  ALE::Obj<ALE::ISieveVisitor::RestrictVisitor<MeshRealSection> > _quadPtsPreV;
-  ALE::Obj<MeshRealSection> _jacobianPre;
-  ALE::Obj<ALE::ISieveVisitor::RestrictVisitor<MeshRealSection> > _jacobianPreV;
-  ALE::Obj<MeshRealSection> _jacobianDetPre;
-  ALE::Obj<ALE::ISieveVisitor::RestrictVisitor<MeshRealSection> > _jacobianDetPreV;
-  ALE::Obj<MeshRealSection> _jacobianInvPre;
-  ALE::Obj<ALE::ISieveVisitor::RestrictVisitor<MeshRealSection> > _jacobianInvPreV;
-  ALE::Obj<MeshRealSection> _basisDerivPre;
-  ALE::Obj<ALE::ISieveVisitor::RestrictVisitor<MeshRealSection> > _basisDerivPreV;
-
-  bool _precomputed; ///< True if we have computed geometry info
-  bool _checkConditioning; ///< True if checking for ill-conditioning
 }; // Quadrature
 
 #include "Quadrature.icc" // inline methods
+#include "Quadrature.cc" // template methods
 
 #endif // pylith_feassemble_quadrature_hh
 
+
 // End of file 

Added: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/QuadratureBase.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/QuadratureBase.cc	                        (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/QuadratureBase.cc	2009-02-06 19:16:53 UTC (rev 14027)
@@ -0,0 +1,185 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+//                           Brad T. Aagaard
+//                        U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#include <portinfo>
+
+#include "QuadratureBase.hh" // implementation of class methods
+
+#include "CellGeometry.hh" // USES CellGeometry
+
+#include <cstring> // USES memcpy()
+#include <cassert> // USES assert()
+#include <stdexcept> // USES std::runtime_error
+#include <sstream> // USES std::ostringstream
+
+// ----------------------------------------------------------------------
+// Constructor
+pylith::feassemble::QuadratureBase::QuadratureBase(void) :
+  _minJacobian(0),
+  _cellDim(0),
+  _numBasis(0),
+  _numQuadPts(0),
+  _spaceDim(0),
+  _geometry(0)
+{ // constructor
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor
+pylith::feassemble::QuadratureBase::~QuadratureBase(void)
+{ // destructor
+  delete _geometry; _geometry = 0;
+} // destructor
+  
+// ----------------------------------------------------------------------
+// Copy constructor
+pylith::feassemble::QuadratureBase::QuadratureBase(const QuadratureBase& q) :
+  _minJacobian(q._minJacobian),
+  _quadPtsRef(q._quadPtsRef),
+  _quadWts(q._quadWts),
+  _basis(q._basis),
+  _basisDerivRef(q._basisDerivRef),
+  _cellDim(q._cellDim),
+  _numBasis(q._numBasis),
+  _numQuadPts(q._numQuadPts),
+  _spaceDim(q._spaceDim),
+  _geometry(0)
+{ // copy constructor
+  if (0 != q._geometry)
+    _geometry = q._geometry->clone();
+} // copy constructor
+
+// ----------------------------------------------------------------------
+// Set basis functions and their derivatives and coordinates and
+//   weights of the quadrature points.
+void
+pylith::feassemble::QuadratureBase::initialize(const double* basis,
+					       const double* basisDerivRef,
+					       const double* quadPtsRef,
+					       const double* quadWts,
+					       const int cellDim,
+					       const int numBasis,
+					       const int numQuadPts,
+					       const int spaceDim)
+{ // initialize
+  if (0 == basis ||
+      0 == basisDerivRef ||
+      0 == quadPtsRef ||
+      0 == quadWts ||
+      cellDim < 0 || cellDim > 3 ||
+      numBasis < 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: " << basisDerivRef << "\n"
+	<< "  quadrature points pointer: " << quadPtsRef << "\n"
+	<< "  quadrature weights pointer: " << quadWts << "\n"
+	<< "  space dimension: " << spaceDim << "\n"
+	<< "  # basis functions: " << numBasis << "\n"
+	<< "  # quadrature points: " << numQuadPts << "\n"
+	<< "  dimension of coordinate space: " << spaceDim << "\n";
+    throw std::runtime_error(msg.str());
+  } // if
+
+  if (cellDim > 0) {
+    int size = numBasis * numQuadPts; assert(size > 0);
+    _basis.resize(size);
+    for (int i=0; i < size; ++i)
+      _basis[i] = basis[i];
+
+    size = numQuadPts * numBasis * cellDim; assert(size > 0);
+    _basisDerivRef.resize(size);
+    for (int i=0; i < size; ++i)
+      _basisDerivRef[i] = basisDerivRef[i];
+
+    size = numQuadPts * cellDim; assert(size > 0);
+    _quadPtsRef.resize(size);
+    for (int i=0; i < size; ++i)
+      _quadPtsRef[i] = quadPtsRef[i];
+
+    size = numQuadPts; assert(size > 0);
+    _quadWts.resize(size);
+    for (int i=0; i < size; ++i)
+      _quadWts[i] = quadWts[i];
+
+    _cellDim = cellDim;
+    _numBasis = numBasis;
+    _numQuadPts = numQuadPts;
+    _spaceDim = spaceDim;
+
+  } 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"
+	  << "  # basis functions: " << numBasis << "\n"
+	  << "  # quadrature points: " << numQuadPts << "\n";
+      throw std::runtime_error(msg.str());
+    } // if
+
+    int size = 1;
+    _basis.resize(size);
+    for (int i=0; i < size; ++i)
+      _basis[i] = basis[i];
+
+    size = 1;
+    _basisDerivRef.resize(size);
+    for (int i=0; i < size; ++i)
+      _basisDerivRef[i] = basisDerivRef[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;
+
+  } // else
+} // initialize
+
+// ----------------------------------------------------------------------
+// Set geometry associated with reference cell.
+void
+pylith::feassemble::QuadratureBase::refGeometry(CellGeometry* const geometry)
+{ // refGeometry
+  delete _geometry; _geometry = (0 != geometry) ? geometry->clone() : 0;
+} // refGeometry
+
+// ----------------------------------------------------------------------
+// Get geometry associated with reference cell.
+const pylith::feassemble::CellGeometry&
+pylith::feassemble::QuadratureBase::refGeometry(void) const
+{ // refGeometry
+  assert(0 != _geometry);
+  return *_geometry;
+} // refGeometry
+
+
+// End of file 

Added: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/QuadratureBase.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/QuadratureBase.hh	                        (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/QuadratureBase.hh	2009-02-06 19:16:53 UTC (rev 14027)
@@ -0,0 +1,233 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+//                           Brad T. Aagaard
+//                        U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+/**
+ * @file pylith/feassemble/QuadratureBase.hh
+ *
+ * @brief Object with basic quadrature information for the reference cell.
+ *
+ * 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. 
+ *
+ * The Quadrature object manages the general functionality of the
+ * numerical quadrature.
+ */
+
+#if !defined(pylith_feassemble_quadraturebase_hh)
+#define pylith_feassemble_quadraturebase_hh
+
+// Include directives ---------------------------------------------------
+#include "feassemblefwd.hh" // forward declarations
+
+#include "pylith/utils/array.hh" // HASA double_array
+
+// Quadrature -----------------------------------------------------------
+class pylith::feassemble::QuadratureBase
+{ // Quadrature
+  friend class TestQuadratureBase; // unit testing
+
+// PUBLIC METHODS ///////////////////////////////////////////////////////
+public :
+
+  /// Constructor
+  QuadratureBase(void);
+
+  /// Destructor
+  ~QuadratureBase(void);
+
+  /** Set basis functions and their derivatives, and coordinates and
+   *  weights of the quadrature points.
+   *
+   * @param basis Array of basis functions evaluated at quadrature pts
+   *   N0Qp0, N1Qp0, ...
+   *   N0Qp1, N1Qp1, ...
+   *   ...
+   *   size = numQuadPts * numBasis
+   *   index = iQuadPt*numBasis + iBasis
+   *
+   * @param basisDerivRef Array of basis function derivaties evaluated at
+   * quadrature pts, where derivatives are with respect to cell's
+   * local coordinates.
+   *   N0pQp0, N0qQp0, N0rQp0, N1pQp0, N1qQp0, N1rQp0, ... 
+   *   N0pQp1, N0qQp1, N0rQp1, N1pQp1, N1qQp1, N1rQp1, ...
+   *   ...
+   *   size = numQuadPts * numBasis * cellDim
+   *   index = iQuadPt*numBasis*cellDim + iBasis*cellDim + iDim
+   *
+   * @param quadPts Array of coordinates of quadrature points in 
+   *   reference cell
+   *   Qp0p, Qp0q, Qp0r
+   *   Qp1p, Qp1q, Qp1r
+   *   size = numQuadPts * numDims
+   *   index = iQuadPt*numDims + iDim
+   *
+   * @param quadWts Array of weights of quadrature points
+   *   WtQp0, WtQp1, ...
+   *   index = iQuadPt
+   *
+   * @param cellDim Number of dimensions in reference cell
+   * @param numBasis Number of basis functions for a cell
+   * @param numQuadPts Number of quadrature points
+   * @param spaceDim Number of dimensions in coordinates of cell vertices
+   */
+  void initialize(const double* basis,
+		  const double* basisDerivRef,
+		  const double* quadPtsRef,
+		  const double* quadWts,
+		  const int cellDim,
+		  const int numBasis,
+		  const int numQuadPts,
+		  const int spaceDim);
+
+  /** Set geometry associated with reference cell.
+   *
+   * @param geometry Geometry of reference cell.
+   */
+  void refGeometry(CellGeometry* const geometry);
+
+  /** Get geometry associated with reference cell.
+   *
+   * @returns Geometry of reference cell.
+   */
+  const CellGeometry& refGeometry(void) const;
+
+  /** Set minimum allowable determinant of Jacobian.
+   *
+   * @param tolerance Minimum allowable value for Jacobian
+   */
+  void minJacobian(const double min);
+
+  /** Get minimum allowable determinant of Jacobian.
+   *
+   * @returns Minimum allowable value for Jacobian
+   */
+  double minJacobian(void) const;
+
+  /** Get coordinates of quadrature points in reference cell.
+   *
+   * @returns Array of coordinates of quadrature points in reference cell.
+   */
+  const double_array& quadPtsRef(void) const;
+
+  /** Get weights of quadrature points.
+   *
+   * @returns Weights of quadrature points
+   */
+  const double_array& quadWts(void) const;
+
+  /** Get basis fns evaluated at quadrature points.
+   *
+   * @returns Array of basis fns evaluated at quadrature points
+   */
+  const double_array& basis(void) const;
+
+  /** Get number of dimensions in reference cell.
+   *
+   * @returns Number of dimensions in reference cell
+   */
+  int cellDim(void) const;
+
+  /** Get number of basis functions for cell.
+   *
+   * @returns Number of basis functions for cell
+   */
+  int numBasis(void) const;
+
+  /** Get number of quadrature points.
+   *
+   * @returns Number of quadrature points
+   */
+  int numQuadPts(void) const;
+
+  /** Get number of dimensions in coordinates of cell vertices.
+   *
+   * @returns Number of dimensions in coordinates of cell vertices
+   */
+  int spaceDim(void) const;
+
+// PROTECTED METHODS ////////////////////////////////////////////////////
+protected :
+
+  /** Copy constructor.
+   *
+   * @param q Quadrature to copy
+   */
+  QuadratureBase(const QuadratureBase& q);
+
+// PROTECTED MEMBERS ////////////////////////////////////////////////////
+protected :
+
+  double _minJacobian; ///< Minium allowable Jacobian determinant
+  
+  /** Array of coordinates of quadrature points in reference cell.
+   *
+   * Reference coordinates: (p,q,r)
+   *
+   * Qp0p, Qp0q, Qp0r
+   * Qp1p, Qp1q, Qp1r
+   *
+   * size = numQuadPts * cellDim
+   * index = iQuadPts*cellDim + iDim
+   */
+  double_array _quadPtsRef;
+
+  /** Array of weights of quadrature points.
+   *
+   * WtQp0, WtQp1, ...
+   * size = numQuadPts
+   * index = iQuadPt
+   */
+  double_array _quadWts;
+
+  /** Array of basis functions evaluated at the quadrature points.
+   *
+   * N0Qp0, N1Qp0, ...
+   * N0Qp1, N1Qp1, ...
+   *
+   * size = numQuadPts * numBasis
+   * index = iQuadPt*numBasis + iBasis
+   */
+  double_array _basis;
+
+  /** Array of basis function derivatives evaluated at the quadrature
+   * points, where derivatives are with respect to cell's local
+   * coordinates.
+   *
+   * N0pQp0, N0qQp0, N0rQp0, N1pQp0, N1qQp0, N1rQp0, ... 
+   * N0pQp1, N0qQp1, N0rQp1, N1pQp1, N1qQp1, N1rQp1, ...
+   *
+   * size = numQuadPts * numBasis * cellDim
+   * index = iQuadPt*numBasis*cellDim + iBasis*cellDim + iDim
+   */
+  double_array _basisDerivRef;
+
+  int _cellDim; ///< Number of dimensions in reference cell
+  int _numBasis; ///< Number of basis functions (and vertices) for cell
+  int _numQuadPts; ///< Number of quadrature points
+  int _spaceDim; ///< Number of dimensions in coordinates of cell vertices
+
+  CellGeometry* _geometry; ///< Geometry of reference cell
+
+// NOT IMPLEMENTED //////////////////////////////////////////////////////
+private :
+
+  const QuadratureBase& operator=(const QuadratureBase&); ///< Not implemented
+
+}; // QuadratureBase
+
+#include "QuadratureBase.icc" // inline methods
+
+#endif // pylith_feassemble_quadraturebase_hh
+
+
+// End of file 

Added: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/QuadratureBase.icc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/QuadratureBase.icc	                        (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/QuadratureBase.icc	2009-02-06 19:16:53 UTC (rev 14027)
@@ -0,0 +1,82 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+//                           Brad T. Aagaard
+//                        U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#if !defined(pylith_feassemble_quadraturebase_hh)
+#error "QuadratureBase.icc must be included only from QuadratureBase.hh"
+#else
+
+// Get minimum allowable Jacobian.
+inline
+double
+pylith::feassemble::QuadratureBase::minJacobian(void) const {
+  return _minJacobian;
+}
+
+// Set minimum allowable Jacobian.
+inline
+void
+pylith::feassemble::QuadratureBase::minJacobian(const double min) {
+  _minJacobian = min;
+}
+
+// Get coordinates of quadrature points in reference cell.
+inline
+const pylith::double_array&
+pylith::feassemble::QuadratureBase::quadPtsRef(void) const {
+  return _quadPtsRef;
+}
+
+// Get weights of quadrature points.
+inline
+const pylith::double_array&
+pylith::feassemble::QuadratureBase::quadWts(void) const {
+  return _quadWts;
+}
+
+// Get basis fns evaluated at quadrature points.
+inline
+const pylith::double_array&
+pylith::feassemble::QuadratureBase::basis(void) const {
+  return _basis;
+}
+
+// Get number of dimensions in reference cell.
+inline
+int
+pylith::feassemble::QuadratureBase::cellDim(void) const {
+  return _cellDim;
+}
+
+// Get number of basis functions for cell.
+inline
+int
+pylith::feassemble::QuadratureBase::numBasis(void) const {
+  return _numBasis;
+}
+
+// Get number of quadrature points.
+inline
+int
+pylith::feassemble::QuadratureBase::numQuadPts(void) const {
+  return _numQuadPts;
+}
+
+// Get number of dimensions in coordinates of cell vertices.
+inline
+int
+pylith::feassemble::QuadratureBase::spaceDim(void) const {
+  return _spaceDim;
+}
+
+#endif
+
+// End of file

Modified: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/feassemblefwd.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/feassemblefwd.hh	2009-02-06 16:26:28 UTC (rev 14026)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/feassemblefwd.hh	2009-02-06 19:16:53 UTC (rev 14027)
@@ -38,17 +38,18 @@
     class GeometryTet3D;
     class GeometryHex3D;
 
-    class Quadrature;
-    class Quadrature0D;
-    class Quadrature1D;
-    class Quadrature1Din2D;
-    class Quadrature1Din3D;
-    class Quadrature2D;
-    class Quadrature2Din3D;
-    class Quadrature3D;
+    class QuadratureBase;
+    template<typename mesh_type> class Quadrature;
+    template<typename mesh_type> class Quadrature0D;
+    template<typename mesh_type> class Quadrature1D;
+    template<typename mesh_type> class Quadrature1Din2D;
+    template<typename mesh_type> class Quadrature1Din3D;
+    template<typename mesh_type> class Quadrature2D;
+    template<typename mesh_type> class Quadrature2Din3D;
+    template<typename mesh_type> class Quadrature3D;
 
     class Constraint;
-    class Integrator;
+    template<typename quadrature_type> class Integrator;
 
     class IntegratorElasticity;
     class ElasticityImplicit;



More information about the CIG-COMMITS mailing list