[cig-commits] r14297 - in short/3D/PyLith/branches/pylith-swig: modulesrc/feassemble pylith/feassemble unittests/pytests/feassemble
brad at geodynamics.org
brad at geodynamics.org
Wed Mar 11 17:10:24 PDT 2009
Author: brad
Date: 2009-03-11 17:10:24 -0700 (Wed, 11 Mar 2009)
New Revision: 14297
Added:
short/3D/PyLith/branches/pylith-swig/unittests/pytests/feassemble/TestMeshQuadrature.py
Removed:
short/3D/PyLith/branches/pylith-swig/unittests/pytests/feassemble/TestQuadrature.py
Modified:
short/3D/PyLith/branches/pylith-swig/modulesrc/feassemble/Makefile.am
short/3D/PyLith/branches/pylith-swig/modulesrc/feassemble/QuadratureRefCell.i
short/3D/PyLith/branches/pylith-swig/modulesrc/feassemble/feassemble.i
short/3D/PyLith/branches/pylith-swig/pylith/feassemble/FIATLagrange.py
short/3D/PyLith/branches/pylith-swig/pylith/feassemble/FIATSimplex.py
short/3D/PyLith/branches/pylith-swig/pylith/feassemble/Quadrature.py
short/3D/PyLith/branches/pylith-swig/unittests/pytests/feassemble/testfeassemble.py
Log:
More work on updating quadrature and integrators.
Modified: short/3D/PyLith/branches/pylith-swig/modulesrc/feassemble/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-swig/modulesrc/feassemble/Makefile.am 2009-03-11 20:27:45 UTC (rev 14296)
+++ short/3D/PyLith/branches/pylith-swig/modulesrc/feassemble/Makefile.am 2009-03-12 00:10:24 UTC (rev 14297)
@@ -27,7 +27,9 @@
GeometryLine2D.i \
GeometryLine3D.i \
GeometryTri2D.i \
- GeometryTri3D.i
+ GeometryTri3D.i \
+ QuadratureRefCell.i \
+ Quadrature.i
swig_generated = \
feassemble_wrap.cxx \
Modified: short/3D/PyLith/branches/pylith-swig/modulesrc/feassemble/QuadratureRefCell.i
===================================================================
--- short/3D/PyLith/branches/pylith-swig/modulesrc/feassemble/QuadratureRefCell.i 2009-03-11 20:27:45 UTC (rev 14296)
+++ short/3D/PyLith/branches/pylith-swig/modulesrc/feassemble/QuadratureRefCell.i 2009-03-12 00:10:24 UTC (rev 14297)
@@ -65,7 +65,30 @@
* @param numQuadPts Number of quadrature points
* @param spaceDim Number of dimensions in coordinates of cell vertices
*/
+ %apply(double* IN_ARRAY2, int DIM1, int DIM2) {
+ (const double* basis,
+ const int numQuadPts,
+ const int numBasis)
+ };
+ %apply(double* IN_ARRAY3, int DIM1, int DIM2, int DIM3) {
+ (const double* basisDerivRef,
+ const int numQuadPts,
+ const int numBasis,
+ const int spaceDim)
+ };
+ %apply(double* IN_ARRAY2, int DIM1, int DIM2) {
+ (const double* quadPtsRef,
+ const int numQuadPts,
+ const int cellDim)
+ };
+ %apply(double* IN_ARRAY1, int DIM1) {
+ (const double* quadWts,
+ const int numQuadPts)
+ };
void initialize(const double* basis,
+ const int numQuadPts1,
+ const int numBasis1,
+ const double
const double* basisDerivRef,
const double* quadPtsRef,
const double* quadWts,
@@ -73,6 +96,10 @@
const int numBasis,
const int numQuadPts,
const int spaceDim);
+ %clear(const double* basis, const int numQuadPts, const int numBasis);
+ %clear(const double* basisDerivRef, const int numQuadPts, const int numBasis, const int spaceDim);
+ %clear(const double* quadPtsRef, const int numQuadPts, const int cellDim);
+ %clear(const double* quadWts, const int numQuadPts);
/** Set geometry associated with reference cell.
*
Modified: short/3D/PyLith/branches/pylith-swig/modulesrc/feassemble/feassemble.i
===================================================================
--- short/3D/PyLith/branches/pylith-swig/modulesrc/feassemble/feassemble.i 2009-03-11 20:27:45 UTC (rev 14296)
+++ short/3D/PyLith/branches/pylith-swig/modulesrc/feassemble/feassemble.i 2009-03-12 00:10:24 UTC (rev 14297)
@@ -33,6 +33,7 @@
#include "pylith/topology/Mesh.hh"
#include "pylith/topology/SubMesh.hh"
#include "pylith/feassemble/Quadrature.hh"
+#include "pylith/feassemble/Integrator.hh"
%}
@@ -73,6 +74,7 @@
%include "QuadratureRefCell.i"
%include "Quadrature.i"
+%include "Integrator.i"
// Template instatiation
%template(MeshQuadrature) pylith::feassemble::Quadrature<pylith::topology::Mesh>;
Modified: short/3D/PyLith/branches/pylith-swig/pylith/feassemble/FIATLagrange.py
===================================================================
--- short/3D/PyLith/branches/pylith-swig/pylith/feassemble/FIATLagrange.py 2009-03-11 20:27:45 UTC (rev 14296)
+++ short/3D/PyLith/branches/pylith-swig/pylith/feassemble/FIATLagrange.py 2009-03-12 00:10:24 UTC (rev 14297)
@@ -56,10 +56,6 @@
import pyre.inventory
- dimension = pyre.inventory.int("dimension", default=3,
- validator=validateDimension)
- dimension.meta['tip'] = "Dimension of finite-element cell."
-
degree = pyre.inventory.int("degree", default=1)
degree.meta['tip'] = "Degree of finite-element cell."
@@ -86,8 +82,8 @@
if self.cellDim > 0:
quadrature = self._setupQuadrature()
- element = self._setupElement()
- dim = self.cellDim
+ element = self._setupElement()
+ dim = self.cellDim
# Get coordinates of vertices (dual basis)
vertices = numpy.array(self._setupVertices(element))
@@ -668,7 +664,6 @@
import FIAT.shapes
ReferenceCell._configure(self)
- self.cellDim = self.inventory.dimension
self.degree = self.inventory.degree
self.order = self.inventory.order
Modified: short/3D/PyLith/branches/pylith-swig/pylith/feassemble/FIATSimplex.py
===================================================================
--- short/3D/PyLith/branches/pylith-swig/pylith/feassemble/FIATSimplex.py 2009-03-11 20:27:45 UTC (rev 14296)
+++ short/3D/PyLith/branches/pylith-swig/pylith/feassemble/FIATSimplex.py 2009-03-12 00:10:24 UTC (rev 14297)
@@ -68,7 +68,7 @@
degree.meta['tip'] = "Degree of finite-element cell."
order = pyre.inventory.int("quad_order", default=-1)
- order.meta['tip'] = "Order of quadrature rule."
+ order.meta['tip'] = "Order of quadrature rule [-1, order = degree]."
# PUBLIC METHODS /////////////////////////////////////////////////////
Modified: short/3D/PyLith/branches/pylith-swig/pylith/feassemble/Quadrature.py
===================================================================
--- short/3D/PyLith/branches/pylith-swig/pylith/feassemble/Quadrature.py 2009-03-11 20:27:45 UTC (rev 14296)
+++ short/3D/PyLith/branches/pylith-swig/pylith/feassemble/Quadrature.py 2009-03-12 00:10:24 UTC (rev 14297)
@@ -20,8 +20,8 @@
from pyre.components.Component import Component
# ----------------------------------------------------------------------
-# Quadrature class
-class Quadrature(Component):
+# QuadratureBase class
+class QuadratureBase(Component):
"""
Python abstract base class for integrating over finite-elements
using quadrature.
@@ -68,37 +68,20 @@
Constructor.
"""
Component.__init__(self, name, facility="quadrature")
- self.minJacobian = 1.0e-06
- self.cppHandle = None
- self.spaceDim = None
return
- def preinitialize(self):
+ def preinitialize(self, spaceDim):
"""
Setup quadrature object.
"""
- self._createCppHandle()
-
- self.cppHandle.minJacobian = self.minJacobian
- self.cppHandle.checkConditioning = self.checkConditioning
-
self._info.log("Initializing reference cell.")
cell = self.cell
- cell.initialize(self.spaceDim)
+ cell.initialize(spaceDim)
- if cell.cellDim != self.cellDim:
- raise TypeError("Dimension of reference cell '%d' does not match "
- "dimension of quadrature implementation '%d'." % \
- (cell.cellDim, self.cellDim))
-
-
self._info.log("Initializing C++ quadrature.")
- self.cppHandle.initialize(cell.basis, cell.basisDeriv,
- cell.quadPts, cell.quadWts,
- cell.cellDim, cell.numCorners, cell.numQuadPts,
- self.spaceDim)
- self.cppHandle.refGeometry = cell.geometry.cppHandle
+ self._initialize(cell)
+ self.refGeometry(cell.geometry)
return
@@ -116,18 +99,47 @@
Set members based using inventory.
"""
Component._configure(self)
- self.minJacobian = self.inventory.minJacobian
- self.checkConditioning = self.inventory.checkConditioning
+ self.minJacobian(self.inventory.minJacobian)
+ self.checkConditioning(self.inventory.checkConditioning)
self.cell = self.inventory.cell
return
- def _createCppHandle(self):
+# ----------------------------------------------------------------------
+from feassemble import MeshQuadrature as ModuleMeshQuadrature
+
+# MeshQuadrature class
+class MeshQuadrature(QuadratureBase, ModuleMeshQuadrature):
+ """
+ Python object for integrating over finite-elements using quadrature.
+ """
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="meshquadrature"):
"""
- Create handle to corresponding C++ object.
+ Constructor.
"""
- raise NotImplementedError("Please implement _createCppHandle() in " \
- "derived class.")
-
-
+ QuadratureBase.__init__(self, name)
+ ModuleMeshQuadrature.__init__(self)
+ return
+
+
+ # PRIVATE METHODS ////////////////////////////////////////////////////
+
+ def _initialize(self, cell):
+ """
+ Initialize C++ quadrature object.
+ """
+ print "cell.basis shape: ", cell.basis.shape
+ print "cell.basisDeriv shape: ", cell.basisDeriv.shape
+ print "cell.quadPts shape: ", cell.quadPts.shape
+ print "cell.quadWts shape: ", cell.quadWts.shape
+ print "cell.geometry.spaceDim: ", cell.geometry.spaceDim()
+ ModuleMeshQuadrature.initialize(self, cell.basis, cell.basisDeriv,
+ cell.quadPts, cell.quadWts,
+ cell.geometry.spaceDim())
+ return
+
+
# End of file
Copied: short/3D/PyLith/branches/pylith-swig/unittests/pytests/feassemble/TestMeshQuadrature.py (from rev 14295, short/3D/PyLith/branches/pylith-swig/unittests/pytests/feassemble/TestQuadrature.py)
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/pytests/feassemble/TestMeshQuadrature.py (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/unittests/pytests/feassemble/TestMeshQuadrature.py 2009-03-12 00:10:24 UTC (rev 14297)
@@ -0,0 +1,167 @@
+#!/usr/bin/env python
+#
+# ======================================================================
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# {LicenseText}
+#
+# ======================================================================
+#
+
+## @file unittests/pytests/feassemble/TestQuadrature.py
+
+## @brief Unit testing of Python Quadrature object.
+
+import unittest
+import numpy
+
+from pylith.feassemble.Quadrature import MeshQuadrature
+from pylith.feassemble.FIATSimplex import FIATSimplex
+from pylith.feassemble.FIATLagrange import FIATLagrange
+
+# ----------------------------------------------------------------------
+def N0(p):
+ return -0.5*p*(1.0-p)
+
+def N0p(p):
+ return -0.5*(1.0-p) + 0.5*p
+
+def N1(p):
+ return 0.5*p*(1.0+p)
+
+def N1p(p):
+ return +0.5*(1.0+p) + 0.5*p
+
+def N2(p):
+ return (1.0-p**2)
+
+def N2p(p):
+ return -2.0*p
+
+# ----------------------------------------------------------------------
+class TestMeshQuadrature(unittest.TestCase):
+ """
+ Unit testing of Python Quadrature object.
+ """
+
+ def test_minJacobian(self):
+ """
+ Test minJacobian().
+ """
+ minJacobian = 4.0e-02;
+ q = MeshQuadrature()
+ q.minJacobian(minJacobian)
+ self.assertEqual(minJacobian, q.minJacobian())
+ return
+
+
+ def test_checkConditioning(self):
+ """
+ Test checkConditioning().
+ """
+ q = MeshQuadrature()
+
+ flag = False # default
+ self.assertEqual(flag, q.checkConditioning())
+
+ flag = True
+ q.checkConditioning(flag)
+ self.assertEqual(flag, q.checkConditioning())
+
+ flag = False
+ q.checkConditioning(flag)
+ self.assertEqual(flag, q.checkConditioning())
+
+ return
+
+
+ def test_initialize(self):
+ """
+ Test initialize().
+ """
+ cell = FIATSimplex()
+ cell.inventory.shape = "line"
+ cell.inventory.degree = 2
+ cell.inventory.order = 2
+ cell._configure()
+
+ verticesE = numpy.array([ [-1.0], [1.0], [0.0] ])
+ quadPtsE = numpy.array( [[-1.0/3**0.5],
+ [+1.0/3**0.5]],
+ dtype=numpy.float64 )
+ quadWtsE = numpy.array( [1.0, 1.0], dtype=numpy.float64 )
+
+ # Compute basis functions and derivatives at quadrature points
+ basisE = numpy.zeros( (2, 3), dtype=numpy.float64)
+ basisDerivE = numpy.zeros( (2, 3, 1), dtype=numpy.float64)
+ iQuad = 0
+ for q in quadPtsE:
+ basisE[iQuad] = numpy.array([N0(q), N1(q), N2(q)],
+ dtype=numpy.float64).reshape( (3,) )
+ deriv = numpy.array([[N0p(q)], [N1p(q)], [N2p(q)]],
+ dtype=numpy.float64)
+ basisDerivE[iQuad] = deriv.reshape((3, 1))
+ iQuad += 1
+
+ quadrature = MeshQuadrature()
+ quadrature.inventory.cell = cell
+ quadrature._configure()
+
+ quadrature.preinitialize(spaceDim=1)
+ quadrature.initialize()
+
+ self.assertEqual(1, quadrature.cellDim())
+ self.assertEqual(1, quadrature.spaceDim())
+ self.assertEqual(3, quadrature.numBasis())
+ self.assertEqual(2, quadrature.numQuadPts())
+
+ from pylith.utils.testarray import test_double
+ #import pylith.feassemble.testfeassemble as testmodule
+
+ #vertices = testmodule.vertices(quadrature.cppHandle)
+ #test_double(self, verticesE, numpy.array(vertices))
+
+ #basis = testmodule.basis(quadrature.cppHandle)
+ #test_double(self, basisE, numpy.array(basis))
+
+ #basisDeriv = testmodule.basisDeriv(quadrature.cppHandle)
+ #test_double(self, basisDerivE, numpy.array(basisDeriv))
+
+ #quadWts = testmodule.quadWts(quadrature.cppHandle)
+ #test_double(self, quadWtsE, numpy.array(quadWts))
+
+ #quadPts = testmodule.quadPtsRef(quadrature.cppHandle)
+ #test_double(self, quadPtsE, numpy.array(quadPts))
+
+ from pylith.feassemble.CellGeometry import GeometryLine1D
+ self.failUnless(isinstance(quadrature.refGeometry(), GeometryLine1D))
+
+ return
+
+
+ def test_simplex1D(self):
+ """
+ Test setup of quadrature for simplex cells for a 1-D problem.
+ """
+ spaceDim = 1
+
+ cell = FIATSimplex()
+ cell.inventory.shape = "line"
+ cell.inventory.degree = 1
+ cell._configure()
+
+ quadrature = MeshQuadrature()
+ quadrature.inventory.cell = cell
+ quadrature._configure()
+
+ quadrature.preinitialize(spaceDim)
+ self.assertEqual(1, quadrature.cellDim())
+ self.assertEqual(spaceDim, quadrature.spaceDim())
+ self.assertEqual(2, quadrature.numBasis())
+
+ return
+
+
+# End of file
Deleted: short/3D/PyLith/branches/pylith-swig/unittests/pytests/feassemble/TestQuadrature.py
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/pytests/feassemble/TestQuadrature.py 2009-03-11 20:27:45 UTC (rev 14296)
+++ short/3D/PyLith/branches/pylith-swig/unittests/pytests/feassemble/TestQuadrature.py 2009-03-12 00:10:24 UTC (rev 14297)
@@ -1,173 +0,0 @@
-#!/usr/bin/env python
-#
-# ======================================================================
-#
-# Brad T. Aagaard
-# U.S. Geological Survey
-#
-# {LicenseText}
-#
-# ======================================================================
-#
-
-## @file unittests/pytests/feassemble/TestQuadrature.py
-
-## @brief Unit testing of Python Quadrature object.
-
-import unittest
-import numpy
-
-from pylith.feassemble.quadrature.Quadrature0D import Quadrature0D
-from pylith.feassemble.quadrature.Quadrature1D import Quadrature1D
-from pylith.feassemble.quadrature.Quadrature1Din2D import Quadrature1Din2D
-from pylith.feassemble.quadrature.Quadrature1Din3D import Quadrature1Din3D
-from pylith.feassemble.quadrature.Quadrature2D import Quadrature2D
-from pylith.feassemble.quadrature.Quadrature2Din3D import Quadrature2Din3D
-from pylith.feassemble.quadrature.Quadrature3D import Quadrature3D
-
-
-# ----------------------------------------------------------------------
-def N0(p):
- return -0.5*p*(1.0-p)
-
-def N0p(p):
- return -0.5*(1.0-p) + 0.5*p
-
-def N1(p):
- return 0.5*p*(1.0+p)
-
-def N1p(p):
- return +0.5*(1.0+p) + 0.5*p
-
-def N2(p):
- return (1.0-p**2)
-
-def N2p(p):
- return -2.0*p
-
-# ----------------------------------------------------------------------
-class TestQuadrature(unittest.TestCase):
- """
- Unit testing of Python Quadrature object.
- """
-
- def test_minJacobian(self):
- """
- Test minJacobian attribute.
- """
- minJacobian = 4.0e-02;
- q = Quadrature1D()
- q.minJacobian = minJacobian
- self.assertEqual(minJacobian, q.minJacobian)
- return
-
-
- def test_initialize(self):
- """
- Test initialize().
- """
- from pylith.feassemble.FIATSimplex import FIATSimplex
- cell = FIATSimplex()
- cell.shape = "line"
- cell.degree = 2
- cell.order = 2
-
- verticesE = numpy.array([ [-1.0], [1.0], [0.0] ])
- quadPtsE = numpy.array( [[-1.0/3**0.5],
- [+1.0/3**0.5]],
- dtype=numpy.float64 )
- quadWtsE = numpy.array( [1.0, 1.0], dtype=numpy.float64 )
-
- # Compute basis functions and derivatives at quadrature points
- basisE = numpy.zeros( (2, 3), dtype=numpy.float64)
- basisDerivE = numpy.zeros( (2, 3, 1), dtype=numpy.float64)
- iQuad = 0
- for q in quadPtsE:
- basisE[iQuad] = numpy.array([N0(q), N1(q), N2(q)],
- dtype=numpy.float64).reshape( (3,) )
- deriv = numpy.array([[N0p(q)], [N1p(q)], [N2p(q)]],
- dtype=numpy.float64)
- basisDerivE[iQuad] = deriv.reshape((3, 1))
- iQuad += 1
-
- quadrature = Quadrature1D()
- quadrature._configure()
- quadrature.cell = cell
-
- quadrature.preinitialize()
- quadrature.initialize()
-
- from pylith.utils.testarray import test_double
- import pylith.feassemble.testfeassemble as testmodule
-
- vertices = testmodule.vertices(quadrature.cppHandle)
- test_double(self, verticesE, numpy.array(vertices))
-
- basis = testmodule.basis(quadrature.cppHandle)
- test_double(self, basisE, numpy.array(basis))
-
- basisDeriv = testmodule.basisDeriv(quadrature.cppHandle)
- test_double(self, basisDerivE, numpy.array(basisDeriv))
-
- quadWts = testmodule.quadWts(quadrature.cppHandle)
- test_double(self, quadWtsE, numpy.array(quadWts))
-
- quadPts = testmodule.quadPtsRef(quadrature.cppHandle)
- test_double(self, quadPtsE, numpy.array(quadPts))
-
- from pylith.feassemble.geometry.GeometryLine1D import GeometryLine1D
- self.failUnless(isinstance(quadrature.cell.geometry, GeometryLine1D))
-
- return
-
-
- def test_constructors(self):
- """
- Test constructors for quadrature objects.
- """
- q = Quadrature0D()
- q._createCppHandle()
- self.assertEqual(1, q.spaceDim)
- self.assertEqual(0, q.cellDim)
- self.failIfEqual(None, q.cppHandle)
-
- q = Quadrature1D()
- q._createCppHandle()
- self.assertEqual(1, q.spaceDim)
- self.assertEqual(1, q.cellDim)
- self.failIfEqual(None, q.cppHandle)
-
- q = Quadrature1Din2D()
- q._createCppHandle()
- self.assertEqual(2, q.spaceDim)
- self.assertEqual(1, q.cellDim)
- self.failIfEqual(None, q.cppHandle)
-
- q = Quadrature1Din3D()
- q._createCppHandle()
- self.assertEqual(3, q.spaceDim)
- self.assertEqual(1, q.cellDim)
- self.failIfEqual(None, q.cppHandle)
-
- q = Quadrature2D()
- q._createCppHandle()
- self.assertEqual(2, q.spaceDim)
- self.assertEqual(2, q.cellDim)
- self.failIfEqual(None, q.cppHandle)
-
- q = Quadrature2Din3D()
- q._createCppHandle()
- self.assertEqual(3, q.spaceDim)
- self.assertEqual(2, q.cellDim)
- self.failIfEqual(None, q.cppHandle)
-
- q = Quadrature3D()
- q._createCppHandle()
- self.assertEqual(3, q.spaceDim)
- self.assertEqual(3, q.cellDim)
- self.failIfEqual(None, q.cppHandle)
-
- return
-
-
-# End of file
Modified: short/3D/PyLith/branches/pylith-swig/unittests/pytests/feassemble/testfeassemble.py
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/pytests/feassemble/testfeassemble.py 2009-03-11 20:27:45 UTC (rev 14296)
+++ short/3D/PyLith/branches/pylith-swig/unittests/pytests/feassemble/testfeassemble.py 2009-03-12 00:10:24 UTC (rev 14297)
@@ -63,8 +63,8 @@
from TestFIATLagrange import TestFIATLagrange
suite.addTest(unittest.makeSuite(TestFIATLagrange))
- #from TestQuadrature import TestQuadrature
- #suite.addTest(unittest.makeSuite(TestQuadrature))
+ from TestMeshQuadrature import TestMeshQuadrature
+ suite.addTest(unittest.makeSuite(TestMeshQuadrature))
#from TestIntegrator import TestIntegrator
#suite.addTest(unittest.makeSuite(TestIntegrator))
More information about the CIG-COMMITS
mailing list