[cig-commits] r4629 - in short/3D/PyLith/trunk: . libsrc/feassemble
modulesrc modulesrc/feassemble pylith/feassemble
baagaard at geodynamics.org
baagaard at geodynamics.org
Tue Sep 26 16:24:51 PDT 2006
Author: baagaard
Date: 2006-09-26 16:24:51 -0700 (Tue, 26 Sep 2006)
New Revision: 4629
Added:
short/3D/PyLith/trunk/modulesrc/feassemble/
short/3D/PyLith/trunk/modulesrc/feassemble/Makefile.am
short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe
Modified:
short/3D/PyLith/trunk/configure.ac
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/Quadrature1Din2D.cc
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din3D.cc
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2D.cc
short/3D/PyLith/trunk/modulesrc/Makefile.am
short/3D/PyLith/trunk/pylith/feassemble/FIATCell.py
short/3D/PyLith/trunk/pylith/feassemble/Quadrature.py
short/3D/PyLith/trunk/pylith/feassemble/ReferenceCell.py
Log:
Added feassemble module. Initial implementation of quadrature for simplex cells using FIAT (not tested).
Modified: short/3D/PyLith/trunk/configure.ac
===================================================================
--- short/3D/PyLith/trunk/configure.ac 2006-09-26 23:23:28 UTC (rev 4628)
+++ short/3D/PyLith/trunk/configure.ac 2006-09-26 23:24:51 UTC (rev 4629)
@@ -129,6 +129,7 @@
libsrc/feassemble/Makefile
libsrc/meshio/Makefile
modulesrc/Makefile
+ modulesrc/feassemble/Makefile
modulesrc/meshio/Makefile
applications/Makefile
unittests/Makefile
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc 2006-09-26 23:23:28 UTC (rev 4628)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc 2006-09-26 23:24:51 UTC (rev 4629)
@@ -217,4 +217,17 @@
_quadPts[i] = 0.0;
} // _resetGeometry
+// ----------------------------------------------------------------------
+// Check determinant of Jacobian against minimum allowable value
+void
+pylith::feassemble::Quadrature::_checkJacobianDet(const double det) const
+{ // _checkJacobianDet
+ if (det < _minJacobian) {
+ std::ostringstream msg;
+ msg << "Determinant of Jacobian (" << det << ") is below minimum\n"
+ << "permissible value (" << _minJacobian << ")!\n";
+ throw std::runtime_error(msg.str());
+ } // if
+} // _checkJacobianDet
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh 2006-09-26 23:23:28 UTC (rev 4628)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh 2006-09-26 23:24:51 UTC (rev 4629)
@@ -72,8 +72,8 @@
* N0xQp0, N0yQp0, N0zQp0, N1xQp0, N1yQp0, N1zQp0, ...
* N0xQp1, N0yQp1, N0zQp1, N1xQp1, N1yQp1, N1zQp1, ...
* ...
- * size = numCorners * numQuadPts * spaceDim
- * index = iQuadPt*numCorners*spaceDim + iBasis*spaceDim + iDim
+ * size = numCorners * numQuadPts * cellDim
+ * index = iQuadPt*numCorners*cellDim + iBasis*cellDim + iDim
*
* @param quadPts Array of coordinates of quadrature points in
* reference cell
@@ -100,17 +100,17 @@
const int numQuadPts,
const int spaceDim);
- /** Set tolerance for minimum allowable Jacobian.
+ /** Set minimum allowable determinant of Jacobian.
*
* @param tolerance Minimum allowable value for Jacobian
*/
- void jacobianTolerance(const double tolerance);
+ void minJacobian(const double min);
- /** Get tolerance for minimum allowable Jacobian.
+ /** Get minimum allowable determinant of Jacobian.
*
* @returns Minimum allowable value for Jacobian
*/
- double jacobianTolerance(void);
+ double minJacobian(void);
// PROTECTED METHODS ////////////////////////////////////////////////////
protected :
@@ -133,6 +133,9 @@
/// Set entries in geometry arrays to zero.
void _resetGeometry(void);
+ /// Check determinant of Jacobian against minimum allowable value
+ void _checkJacobianDet(const double det) const;
+
// PRIVATE METHODS //////////////////////////////////////////////////////
private :
@@ -141,7 +144,7 @@
// PROTECTED MEMBERS ////////////////////////////////////////////////////
protected :
- double _jacobianTol; ///< Tolerance for minium allowable Jacobian determinant
+ double _minJacobian; ///< Minium allowable Jacobian determinant
/** Array of basis functions evaluated at the quadrature points.
*
@@ -158,8 +161,8 @@
* N0xQp0, N0yQp0, N0zQp0, N1xQp0, N1yQp0, N1zQp0, ...
* N0xQp1, N0yQp1, N0zQp1, N1xQp1, N1yQp1, N1zQp1, ...
*
- * size = numQuadPts * numCorners * spaceDim
- * index = iQuadPt*numCorners*spaceDim + iBasis*spaceDim + iDim
+ * size = numQuadPts * numCorners * cellDim
+ * index = iQuadPt*numCorners*cellDim + iBasis*cellDim + iDim
*/
double* _basisDeriv;
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.icc 2006-09-26 23:23:28 UTC (rev 4628)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.icc 2006-09-26 23:24:51 UTC (rev 4629)
@@ -14,18 +14,18 @@
#error "Quadrature.icc must be included only from Quadrature.hh"
#else
-// Get tolerance for minimum allowable Jacobian.
+// Get minimum allowable Jacobian.
inline
double
-pylith::feassemble::Quadrature::jacobianTolerance(void) {
- return _jacobianTol;
+pylith::feassemble::Quadrature::minJacobian(void) {
+ return _minJacobian;
}
-// Set tolerance for minimum allowable Jacobian.
+// Set minimum allowable Jacobian.
inline
void
-pylith::feassemble::Quadrature::jacobianTolerance(const double tolerance) {
- _jacobianTol = tolerance;
+pylith::feassemble::Quadrature::minJacobian(const double min) {
+ _minJacobian = min;
}
#endif
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.cc 2006-09-26 23:23:28 UTC (rev 4628)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.cc 2006-09-26 23:24:51 UTC (rev 4629)
@@ -15,8 +15,6 @@
#include "Quadrature1D.hh" // implementation of class methods
#include <assert.h> // USES assert()
-#include <stdexcept> // USES std::runtime_error
-#include <sstream> // USES std::ostringstream
// ----------------------------------------------------------------------
// Constructor
@@ -59,7 +57,8 @@
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()));
+ assert(1 == coordinates.GetFiberDimensionByDepth(patch,
+ *vertices->begin(), 0));
// Loop over quadrature points
for (int iQuadPt=0; iQuadPt < _numQuadPts; ++iQuadPt) {
@@ -79,12 +78,7 @@
// Compute determinant of Jacobian at quadrature point
// |J| = j00
const double det = _jacobian[iQuadPt];
- if (det < _jacobianTol) {
- std::ostringstream msg;
- msg << "Determinant of Jacobian (" << det << ") is below minimum\n"
- << "permissible value (" << _jacobianTol << ")!\n";
- throw std::runtime_error(msg.str());
- } // for
+ _checkJacobianDet(det);
_jacobianDet[iQuadPt] = _jacobian[iQuadPt];
// Compute inverse of Jacobian at quadrature point
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din2D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din2D.cc 2006-09-26 23:23:28 UTC (rev 4628)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din2D.cc 2006-09-26 23:24:51 UTC (rev 4629)
@@ -15,8 +15,6 @@
#include "Quadrature1Din2D.hh" // implementation of class methods
#include <assert.h> // USES assert()
-#include <stdexcept> // USES std::runtime_error
-#include <sstream> // USES std::ostringstream
// ----------------------------------------------------------------------
// Constructor
@@ -59,7 +57,8 @@
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()));
+ assert(2 == coordinates.GetFiberDimensionByDepth(patch,
+ *vertices->begin(), 0));
// Loop over quadrature points
for (int iQuadPt=0; iQuadPt < _numQuadPts; ++iQuadPt) {
@@ -97,12 +96,7 @@
for (int iDim=0, iJ=iQuadPt*_spaceDim; iDim < _spaceDim; ++iDim)
det += _jacobian[iJ+iDim]*_jacobian[iJ+iDim];
det = sqrt(det);
- if (det < _jacobianTol) {
- std::ostringstream msg;
- msg << "Determinant of Jacobian (" << det << ") is below minimum\n"
- << "permissible value (" << _jacobianTol << ")!\n";
- throw std::runtime_error(msg.str());
- } // for
+ _checkJacobianDet(det);
_jacobianDet[iQuadPt] = det;
// Compute inverse of Jacobian at quadrature point
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din3D.cc 2006-09-26 23:23:28 UTC (rev 4628)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din3D.cc 2006-09-26 23:24:51 UTC (rev 4629)
@@ -15,8 +15,6 @@
#include "Quadrature1Din3D.hh" // implementation of class methods
#include <assert.h> // USES assert()
-#include <stdexcept> // USES std::runtime_error
-#include <sstream> // USES std::ostringstream
// ----------------------------------------------------------------------
// Constructor
@@ -59,7 +57,8 @@
const ALE::Mesh::topology_type::patch_type patch = 0;
const ALE::Mesh::section_type::value_type* vertCoords =
coordinates->restrict(patch, cell);
- //assert(3 == coordinates.GetFiberDimension(patch, *vertices->begin()));
+ assert(3 == coordinates.GetFiberDimensionByDepth(patch,
+ *vertices->begin(), 0));
// Loop over quadrature points
for (int iQuadPt=0; iQuadPt < _numQuadPts; ++iQuadPt) {
@@ -99,12 +98,7 @@
for (int iDim=0, iJ=iQuadPt*_spaceDim; iDim < _spaceDim; ++iDim)
det += _jacobian[iJ+iDim]*_jacobian[iJ+iDim];
det = sqrt(det);
- if (det < _jacobianTol) {
- std::ostringstream msg;
- msg << "Determinant of Jacobian (" << det << ") is below minimum\n"
- << "permissible value (" << _jacobianTol << ")!\n";
- throw std::runtime_error(msg.str());
- } // for
+ _checkJacobianDet(det);
_jacobianDet[iQuadPt] = det;
// Compute inverse of Jacobian at quadrature point
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2D.cc 2006-09-26 23:23:28 UTC (rev 4628)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2D.cc 2006-09-26 23:24:51 UTC (rev 4629)
@@ -15,8 +15,6 @@
#include "Quadrature2D.hh" // implementation of class methods
#include <assert.h> // USES assert()
-#include <stdexcept> // USES std::runtime_error
-#include <sstream> // USES std::ostringstream
// ----------------------------------------------------------------------
// Constructor
@@ -59,7 +57,8 @@
const ALE::Mesh::topology_type::patch_type patch = 0;
const ALE::Mesh::section_type::value_type* vertCoords =
coordinates->restrict(patch, cell);
- //assert(3 == coordinates.GetFiberDimension(patch, *vertices->begin()));
+ assert(2 == coordinates.GetFiberDimensionByDepth(patch,
+ *vertices->begin(), 0));
// Loop over quadrature points
for (int iQuadPt=0; iQuadPt < _numQuadPts; ++iQuadPt) {
@@ -88,7 +87,7 @@
iVertex < _numCorners;
++iVertex)
for (int iRow=0,
- iB=iQuadPt*_numCorners*_spaceDim+iVertex*_spaceDim;
+ iB=iQuadPt*_numCorners*_spaceDim+iVertex*_cellDim;
iRow < _cellDim;
++iRow) {
const double deriv = _basisDeriv[iB+iRow];
@@ -106,12 +105,7 @@
const int i10 = iJ + 1*_spaceDim + 0;
const double det =
_jacobian[i00]*_jacobian[i11] - _jacobian[i01]*_jacobian[i10];
- if (det < _jacobianTol) {
- std::ostringstream msg;
- msg << "Determinant of Jacobian (" << det << ") is below minimum\n"
- << "permissible value (" << _jacobianTol << ")!\n";
- throw std::runtime_error(msg.str());
- } // for
+ _checkJacobianDet(det);
_jacobianDet[iQuadPt] = det;
// Compute inverse of Jacobian at quadrature point
Modified: short/3D/PyLith/trunk/modulesrc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/modulesrc/Makefile.am 2006-09-26 23:23:28 UTC (rev 4628)
+++ short/3D/PyLith/trunk/modulesrc/Makefile.am 2006-09-26 23:24:51 UTC (rev 4629)
@@ -11,6 +11,7 @@
#
SUBDIRS = \
+ feassemble \
meshio
# version
Added: short/3D/PyLith/trunk/modulesrc/feassemble/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/modulesrc/feassemble/Makefile.am 2006-09-26 23:23:28 UTC (rev 4628)
+++ short/3D/PyLith/trunk/modulesrc/feassemble/Makefile.am 2006-09-26 23:24:51 UTC (rev 4629)
@@ -0,0 +1,40 @@
+# -*- Makefile -*-
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+subpackage = feassemble
+include $(top_srcdir)/subpackage.am
+
+subpkgpyexec_LTLIBRARIES = feassemblemodule.la
+
+feassemblemodule_la_LDFLAGS = -module
+
+feassemblemodule_la_SOURCES = feassemble.pyxe
+
+nodist_feassemblemodule_la_SOURCES = \
+ feassemble.c feassemble_embed.cpp feassemble_embed.h
+
+feassemblemodule_la_LIBADD = \
+ $(top_builddir)/libsrc/feassemble/libpylithfeassemble.la $(PETSC_LIB)
+
+INCLUDES += -I$(PYTHON_INCDIR) $(PETSC_INCLUDE)
+
+feassemble.pyx feassemble_embed.cpp feassemble_embed.h: feassemble.pyxe
+ if [ "${VPATH}" != "" ]; then cp $< .; fi && pyrexembed feassemble.pyxe && if [ "${VPATH}" != "" ]; then rm -f feassemble.pyxe; fi
+feassemble_embed.cpp: feassemble_embed.h
+feassemble_embed.h: feassemble.pyx
+
+.pyx.c:
+ pyrexc $<
+
+CLEANFILES = feassemble.pyx feassemble.c *_embed.*
+
+# End of file
Added: short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe
===================================================================
--- short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe 2006-09-26 23:23:28 UTC (rev 4628)
+++ short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe 2006-09-26 23:24:51 UTC (rev 4629)
@@ -0,0 +1,210 @@
+#!/usr/bin/env python
+#
+# ======================================================================
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# {LicenseText}
+#
+# ======================================================================
+#
+
+#header{
+#include "pylith/feassemble/Quadrature.hh"
+#include "pylith/feassemble/Quadrature1D.hh"
+#include "pylith/feassemble/Quadrature1Din2D.hh"
+#include "pylith/feassemble/Quadrature1Din3D.hh"
+#include "pylith/feassemble/Quadrature2D.hh"
+#}header
+
+# ----------------------------------------------------------------------
+cdef extern from "Python.h":
+ object PyCObject_FromVoidPtr(void*, void (*destruct)(void*))
+ void* PyCObject_AsVoidPtr(object)
+
+cdef void* ptrFromHandle(obj):
+ """Extract pointer from PyCObject."""
+ return PyCObject_AsVoidPtr(obj.handle)
+
+cdef extern from "stdlib.h":
+ ctypedef unsigned long size_t
+ void* malloc(size_t size)
+ void free(void* mem)
+
+# ----------------------------------------------------------------------
+cdef class Quadrature:
+
+ cdef void* thisptr # Pointer to C++ object
+ cdef readonly object handle # PyCObject holding pointer to C++ object
+ cdef readonly object name # Identifier for object base type
+
+ def __init__(self):
+ """
+ Constructor.
+ """
+ self.handle = None
+ self.thisptr = NULL
+ self.name = "pylith_feassemble_Quadrature"
+ return
+
+
+ def initialize(self, basis, basisDeriv, quadPts, quadWts,
+ cellDim, numCorners, numQuadPts, spaceDim):
+ """
+ Set basis functions and their derivatives, and coordinates and
+ weights of quadrature points.
+
+ @param basis Basis functions evaluated at the quadrature points
+ @param basisDeriv Basis function derivatives evaluated at quad pts
+ @param quadPts Coordinates of quadrature points in reference cell
+ @param quadWts Weights of quadrature points
+ @param cellDim Dimension of reference cell
+ @param numCorners Number of vertices in reference cell
+ @param numQuadPts Number of quadrature points
+ @param spaceDim Number of dimensions associated with cell vertices
+ """
+ # create shim for method 'initialize'
+ #embed{ void Quadrature_initialize(void* pObj, double* basis, double* basisDeriv, double* quadPts, double* quadWts, int cellDim, int numCorners, int numQuadPts, int spaceDim)
+ ((pylith::feassemble::Quadrature*) pObj)->initialize(basis, basisDeriv,
+ quadPts, quadWts,
+ cellDim, numCorners,
+ numQuadPts, spaceDim);
+ #}embed
+
+ import spatialdata.utils.simplearray
+ basis = spatialdata.utils.simplearray.objAsSimpleArray(basis)
+ if not basis.isCompatible(nd=2, simpletype="double",
+ contiguous=True, notswapped=True):
+ raise TypeError, \
+ "Argument 'basis' must be a contiguous, 2-D array " \
+ "of type double."
+ basis = spatialdata.utils.simplearray.objAsSimpleArray(basisDeriv)
+ if not basisDeriv.isCompatible(nd=3, simpletype="double",
+ contiguous=True, notswapped=True):
+ raise TypeError, \
+ "Argument 'basisDeriv' must be a contiguous, 3-D array " \
+ "of type double."
+ quadPts = spatialdata.utils.simplearray.objAsSimpleArray(quadPts)
+ if not quadPts.isCompatible(nd=2, simpletype="double",
+ contiguous=True, notswapped=True):
+ raise TypeError, \
+ "Argument 'quadPts' must be a contiguous, 2-D array " \
+ "of type double."
+ quadWts = spatialdata.utils.simplearray.objAsSimpleArray(quadWts)
+ if not quadWts.isCompatible(nd=1, simpletype="double",
+ contiguous=True, notswapped=True):
+ raise TypeError, \
+ "Argument 'quadWts' must be a contiguous, 2-D array " \
+ "of type double."
+
+ cdef double* basisCpp
+ cdef double* basisDerivCpp
+ cdef double* quadPtsCpp
+ cdef double* quadWtsCpp
+ basisCpp = <double*> PyCObject_AsVoidPtr(basis.data)
+ basisDerivCpp = <double*> PyCObject_AsVoidPtr(basisDeriv.data)
+ quadPtsCpp = <double*> PyCObject_AsVoidPtr(quadPts.data)
+ quadWtsCpp = <double*> PyCObject_AsVoidPtr(quadWts.data)
+
+ Quadrature_initialize(self.thisptr, basisCpp, basisDerivCpp,
+ quadPtsCpp, quadWtsCpp,
+ cellDim, numCorners,
+ numQuadPts, spaceDim)
+ return
+
+
+ def _createHandle(self):
+ """Wrap pointer to C++ object in PyCObject."""
+ # create shim for destructor
+ #embed{ void Quadrature_destructor(void* pObj)
+ pylith::feassemble::Quadrature* pQ =
+ (pylith::feassemble::Quadrature*) pObj;
+ delete pQ;
+ #}embed
+ return PyCObject_FromVoidPtr(self.thisptr, Quadrature_destructor)
+
+
+ property minJacobian:
+ def __set__(self, min):
+ """Set minimum allowable Jacobian."""
+ # create shim for method 'minJacobian'
+ #embed{ void Quadrature_minJacobian_set(void* pObj, double min)
+ ((pylith::feassemble::Quadrature*) pObj)->minJacobian(min);
+ #}embed
+ Quadrature_minJacobian_set(self.thisptr, min)
+
+ def __get__(self):
+ """Get minimum allowable Jacobian."""
+ # create shim for method 'minJacobian'
+ #embed{ double Quadrature_minJacobian_get(void* pObj)
+ return ((pylith::feassemble::Quadrature*) pObj)->minJacobian();
+ #}embed
+ return Quadrature_minJacobian_get(self.thisptr)
+
+
+# ----------------------------------------------------------------------
+cdef class Quadrature1D(Quadrature):
+
+ def __init__(self):
+ """Constructor."""
+ # create shim for constructor
+ #embed{ void* Quadrature1D_constructor()
+ return (void*)(new pylith::feassemble::Quadrature1D);
+ #}embed
+
+ Quadrature.__init__(self)
+ self.thisptr = Quadrature1D_constructor()
+ self.handle = self._createHandle()
+ return
+
+
+# ----------------------------------------------------------------------
+cdef class Quadrature1Din2D(Quadrature):
+
+ def __init__(self):
+ """Constructor."""
+ # create shim for constructor
+ #embed{ void* Quadrature1Din2D_constructor()
+ return (void*)(new pylith::feassemble::Quadrature1Din2D);
+ #}embed
+
+ Quadrature.__init__(self)
+ self.thisptr = Quadrature1Din2D_constructor()
+ self.handle = self._createHandle()
+ return
+
+
+# ----------------------------------------------------------------------
+cdef class Quadrature1Din3D(Quadrature):
+
+ def __init__(self):
+ """Constructor."""
+ # create shim for constructor
+ #embed{ void* Quadrature1Din3D_constructor()
+ return (void*)(new pylith::feassemble::Quadrature1Din3D);
+ #}embed
+
+ Quadrature.__init__(self)
+ self.thisptr = Quadrature1Din3D_constructor()
+ self.handle = self._createHandle()
+ return
+
+
+# ----------------------------------------------------------------------
+cdef class Quadrature2D(Quadrature):
+
+ def __init__(self):
+ """Constructor."""
+ # create shim for constructor
+ #embed{ void* Quadrature2D_constructor()
+ return (void*)(new pylith::feassemble::Quadrature2D);
+ #}embed
+
+ Quadrature.__init__(self)
+ self.thisptr = Quadrature2D_constructor()
+ self.handle = self._createHandle()
+ return
+
+
+# End of file
Modified: short/3D/PyLith/trunk/pylith/feassemble/FIATCell.py
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/FIATCell.py 2006-09-26 23:23:28 UTC (rev 4628)
+++ short/3D/PyLith/trunk/pylith/feassemble/FIATCell.py 2006-09-26 23:24:51 UTC (rev 4629)
@@ -47,12 +47,19 @@
self.basis = numpy.array(basisFns.tabulate(quadpts)).transpose()
# Evaluate derivatives of basis functions at quadrature points
- # ADD STUFF HERE
+ import FIAT.shapes
+ dim = FIAT.shapes.dimension(basisFns.base.shape)
+ self.basisDeriv = numpy.array([basisFns.deriv_all(d).tabulate(quadpts) \
+ for d in range(dim)]).squeeze().transpose()
+ self.cellDim = dim
+ self.numCorners
+ self.numQuadPts = len(quadrature.get_weights())
+
self._info.line("Basis:")
self._info.line(self.basis)
- #print "Basis derivatives:"
- #print self.basisDeriv
+ self._info.line("Basis derivatives:")
+ self._info.line(self.basisDeriv)
self._info.line("Quad pts:")
self._info.line(quadrature.get_points())
self._info.line("Quad wts:")
Modified: short/3D/PyLith/trunk/pylith/feassemble/Quadrature.py
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/Quadrature.py 2006-09-26 23:23:28 UTC (rev 4628)
+++ short/3D/PyLith/trunk/pylith/feassemble/Quadrature.py 2006-09-26 23:24:51 UTC (rev 4629)
@@ -15,19 +15,9 @@
## @brief Python abstract base class for integrating over
## finite-elements using quadrature.
-# DESIGN QUESTION
-#
-# The dimension of the space associated with the coordinates of
-# the vertices is specific to the quadrature object. However, the
-# quadrature points/weights and basis functions are associated with
-# the reference cell not the actual cell.
-#
-# Where should the dimension of the space associated with the
-# coordinates of the vertices be specified? I don't think the user
-# will specify it, so that means it won't be in the inventory.
-
from pyre.components.Component import Component
+# ----------------------------------------------------------------------
# Quadrature class
class Quadrature(Component):
"""
@@ -44,18 +34,18 @@
## Python object for managing Quadrature facilities and properties.
##
## \b Properties
- ## @li \b jacobian_tolerance Minimum allowable determinant of Jacobian.
+ ## @li \b min_jacobian Minimum allowable determinant of Jacobian.
##
## \b Facilities
## @li \b cell Reference cell with basis functions and quadrature rules
import pyre.inventory
- jacobianTol = pyre.inventory.float("jacobian_tolerance", default=1.0e-06)
- jacobianTol.meta['tip'] = "Minimum allowable determinant of Jacobian."
+ minJacobian = pyre.inventory.float("min_jacobian", default=1.0e-06)
+ minJacobian.meta['tip'] = "Minimum allowable determinant of Jacobian."
- from CellFIAT import CellFIAT
- cell = pyre.inventory.facility("cell", factory=CellFIAT)
+ from FIATSimplex import FIATSimplex
+ cell = pyre.inventory.facility("cell", factory=FIATSimplex)
cell.meta['tip'] = "Reference cell with basis fns and quadrature rules."
@@ -66,22 +56,23 @@
Constructor.
"""
Component.__init__(self, name, facility="quadrature")
- #import pylith.feassemble.feassemble as bindings
- #self.cppHandle = bindings.Quadrature()
- self.spaceDim = 3
+ self.cppHandle = None
+ self.spaceDim = None
return
def initialize(self):
"""
Initialize C++ quadrature object.
"""
- self.cell.initialize()
-
- # Set minimum allowable determinant of Jacobian
- #self.cppHandle.jacobianTol = self.jacobianTol
+ self.cppHandle.minJacobian = self.minJacobian
- # Get basis functions, quadrature points
- #self.cppHandle.initialize(LOTS OF ARGS)
+ c = self.cell
+ c.initialize()
+ self.cppHandle.initialize(c.basis, c.basisDeriv,
+ c.quadrature.get_points(),
+ c.quadrature.get_weights(),
+ c.cellDim, c.numCorners, c.numQuadPts,
+ self.spaceDim)
return
@@ -92,12 +83,85 @@
Set members based using inventory.
"""
Component._configure(self)
- self.jacobianTol = self.inventory.jacobianTol
+ self.minJacobian = self.inventory.minJacobian
self.cell = self.inventory.cell
return
-# version
-__id__ = "$Id$"
+# ----------------------------------------------------------------------
+# Quadrature1D class
+class Quadrature1D(Quadrature):
+ """
+ Python object for integrating over 1-D finite-elements in a 1-D
+ domain using quadrature.
+ """
+ def __init__(self, name="quadrature1d"):
+ """
+ Constructor.
+ """
+ Quadrature.__init(self, name)
+ import pylith.feassemble.feassemble as bindings
+ self.cppHandle = bindings.Quadrature1D()
+ self.spaceDim = 1
+ return
+
+
+# ----------------------------------------------------------------------
+# Quadrature1Din2D class
+class Quadrature1Din2D(Quadrature):
+ """
+ Python object for integrating over 1-D finite-elements in a 2-D
+ domain using quadrature.
+ """
+
+ def __init__(self, name="quadrature1din2d"):
+ """
+ Constructor.
+ """
+ Quadrature.__init(self, name)
+ import pylith.feassemble.feassemble as bindings
+ self.cppHandle = bindings.Quadrature1Din2D()
+ self.spaceDim = 2
+ return
+
+
+# ----------------------------------------------------------------------
+# Quadrature1Din3D class
+class Quadrature1Din3D(Quadrature):
+ """
+ Python object for integrating over 1-D finite-elements in a 3-D
+ domain using quadrature.
+ """
+
+ def __init__(self, name="quadrature1din3d"):
+ """
+ Constructor.
+ """
+ Quadrature.__init(self, name)
+ import pylith.feassemble.feassemble as bindings
+ self.cppHandle = bindings.Quadrature1Din3D()
+ self.spaceDim = 3
+ return
+
+
+# ----------------------------------------------------------------------
+# Quadrature2D class
+class Quadrature2D(Quadrature):
+ """
+ Python object for integrating over 2-D finite-elements in a 2-D
+ domain using quadrature.
+ """
+
+ def __init__(self, name="quadrature2d"):
+ """
+ Constructor.
+ """
+ Quadrature.__init(self, name)
+ import pylith.feassemble.feassemble as bindings
+ self.cppHandle = bindings.Quadrature at D()
+ self.spaceDim = 2
+ return
+
+
# End of file
Modified: short/3D/PyLith/trunk/pylith/feassemble/ReferenceCell.py
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/ReferenceCell.py 2006-09-26 23:23:28 UTC (rev 4628)
+++ short/3D/PyLith/trunk/pylith/feassemble/ReferenceCell.py 2006-09-26 23:24:51 UTC (rev 4629)
@@ -34,6 +34,9 @@
self.basis = None
self.basisDeriv = None
self.quadrature = None
+ self.cellDim = None
+ self.numCorners = None
+ self.numQuadPts = None
return
def initialize(self):
More information about the cig-commits
mailing list