[cig-commits] r7080 - in short/3D/PyLith/trunk: . pylith/feassemble
unittests/pytests/feassemble
brad at geodynamics.org
brad at geodynamics.org
Wed Jun 6 12:51:08 PDT 2007
Author: brad
Date: 2007-06-06 12:51:08 -0700 (Wed, 06 Jun 2007)
New Revision: 7080
Modified:
short/3D/PyLith/trunk/TODO
short/3D/PyLith/trunk/pylith/feassemble/FIATLagrange.py
short/3D/PyLith/trunk/pylith/feassemble/FIATSimplex.py
short/3D/PyLith/trunk/unittests/pytests/feassemble/TestFIATLagrange.py
short/3D/PyLith/trunk/unittests/pytests/feassemble/TestFIATSimplex.py
Log:
Added vertices (dual basis) to FIATLagrange and FIATSimplex objects. Updated (Python) unit tests and added tests for more cells.
Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO 2007-06-06 18:04:51 UTC (rev 7079)
+++ short/3D/PyLith/trunk/TODO 2007-06-06 19:51:08 UTC (rev 7080)
@@ -4,41 +4,30 @@
1. Simple tests with analytical solutions
- b. 2-D
+ a. 2-D
i. tri3 cells
(1) axial compression
(2) shear
- (3) axialshear
ii. quad4 cells
(1) axial compression
(2) shear
- (3) axialshear
- c. 3-D
+ b. 3-D
i. tet4 cells
(1) axial compression
(2) shear
- (3) axialshear
ii. hex8 cells
(1) axial compression
(2) shear
- (3) axialshear
-2. Allow use of all elasticity constants (9 for 2-D, 36 for 3-D).
- a. Materials C++ code
- b. Integrator C++ code
- b. Material C++ unit tests
-
-3. Add dualBasis to Quadrature.
+2. Add dualBasis to Quadrature.
a. Python
- ReferenceCell
- FIATSimplex
Quadrature()
b. C++
Quadrature
c. C++ unit tests
d. Python unit tests
-4. Implement faults for kinematic source
+3. Implement faults for kinematic source
a. Creation of cohesive cells
i. Add tests for interpolated meshes.
Double check consistency in ordering of vertices (positive/negative).
@@ -83,6 +72,11 @@
constructor
initialize()
+3. Allow use of all elasticity constants (9 for 2-D, 36 for 3-D).
+ a. Materials C++ code
+ b. Integrator C++ code
+ b. Material C++ unit tests
+
4. Implement absorbing boundary conditions
======================================================================
Modified: short/3D/PyLith/trunk/pylith/feassemble/FIATLagrange.py
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/FIATLagrange.py 2007-06-06 18:04:51 UTC (rev 7079)
+++ short/3D/PyLith/trunk/pylith/feassemble/FIATLagrange.py 2007-06-06 19:51:08 UTC (rev 7080)
@@ -86,6 +86,9 @@
element = self._setupElement()
dim = self.cellDim
+ # Get coordinates of vertices (dual basis)
+ vertices = numpy.array(self._setupVertices(element))
+
# Evaluate basis functions at quadrature points
quadpts = numpy.array(quadrature.get_points())
quadwts = numpy.array(quadrature.get_weights())
@@ -101,12 +104,19 @@
self.numCorners = numBasisFns**dim
if dim == 1:
+ self.vertices = numpy.array(vertices)
self.quadPts = quadpts
self.quadWts = quadwts
self.basis = basis
self.basisDeriv = basisDeriv
else:
if dim == 2:
+ self.vertices = numpy.zeros((numBasisFns, numBasisFns, dim))
+ for q in range(numBasisFns):
+ for r in range(numBasisFns):
+ self.vertices[q][r][0] = vertices[r]
+ self.vertices[q][r][1] = vertices[q]
+
self.quadPts = numpy.zeros((numQuadPts, numQuadPts, dim))
self.quadWts = numpy.zeros((numQuadPts, numQuadPts))
self.basis = numpy.zeros((numQuadPts, numQuadPts,
@@ -124,6 +134,15 @@
self.basisDeriv[q][r][f][g][0] = basisDeriv[r][g][0]*basis[q][f]
self.basisDeriv[q][r][f][g][1] = basis[r][g]*basisDeriv[q][f][0]
elif dim == 3:
+ self.vertices = numpy.zeros((numBasisFns, numBasisFns, numBasisFns,
+ dim))
+ for q in range(numBasisFns):
+ for r in range(numBasisFns):
+ for s in range(numBasisFns):
+ self.vertices[q][r][s][0] = vertices[s]
+ self.vertices[q][r][s][1] = vertices[r]
+ self.vertices[q][r][s][2] = vertices[q]
+
self.quadPts = numpy.zeros((numQuadPts, numQuadPts,
numQuadPts, dim))
self.quadWts = numpy.zeros((numQuadPts, numQuadPts, numQuadPts))
@@ -146,19 +165,23 @@
self.basisDeriv[q][r][s][f][g][h][0] = basisDeriv[s][h][0]*basis[r][g]*basis[q][f]
self.basisDeriv[q][r][s][f][g][h][1] = basis[s][h]*basisDeriv[r][g][0]*basis[q][f]
self.basisDeriv[q][r][s][f][g][h][2] = basis[s][h]*basis[r][g]*basisDeriv[q][f][0]
+ self.vertices = numpy.reshape(self.vertices, (self.numCorners, dim))
self.quadPts = numpy.reshape(self.quadPts, (self.numQuadPts, dim))
self.quadWts = numpy.reshape(self.quadWts, (self.numQuadPts))
self.basis = numpy.reshape(self.basis, (self.numQuadPts, self.numCorners))
self.basisDeriv = numpy.reshape(self.basisDeriv, (self.numQuadPts, self.numCorners, dim))
- self._info.line("Basis (quad pts):")
+
+ self._info.line("Vertices: ")
+ self._info.line(self.vertices)
+ self._info.line("Quad pts:")
+ self._info.line(quadrature.get_points())
+ self._info.line("Quad wts:")
+ self._info.line(quadrature.get_weights())
+ self._info.line("Basis fns @ quad pts ):")
self._info.line(self.basis)
- self._info.line("Basis derivatives (quad pts):")
+ self._info.line("Basis fn derivatives @ quad pts:")
self._info.line(self.basisDeriv)
- self._info.line("Quad pts:")
- self._info.line(self.quadPts)
- self._info.line("Quad wts:")
- self._info.line(self.quadWts)
- self._info.log()
+ self._info.log()
return
@@ -198,6 +221,13 @@
return FIAT.Lagrange.Lagrange(FIAT.shapes.LINE, self.degree)
+ def _setupVertices(self, element):
+ """
+ Setup evaluation functional points for reference cell.
+ """
+ return element.Udual.pts
+
+
# FACTORIES ////////////////////////////////////////////////////////////
def reference_cell():
Modified: short/3D/PyLith/trunk/pylith/feassemble/FIATSimplex.py
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/FIATSimplex.py 2007-06-06 18:04:51 UTC (rev 7079)
+++ short/3D/PyLith/trunk/pylith/feassemble/FIATSimplex.py 2007-06-06 19:51:08 UTC (rev 7080)
@@ -86,7 +86,10 @@
"""
quadrature = self._setupQuadrature()
basisFns = self._setupBasisFns()
-
+
+ # Get coordinates of vertices (dual basis)
+ self.vertices = numpy.array(self._setupVertices(), dtype=numpy.float64)
+
# Evaluate basis functions at quadrature points
quadpts = quadrature.get_points()
basis = numpy.array(basisFns.tabulate(quadpts)).transpose()
@@ -106,14 +109,16 @@
self.numCorners = len(basisFns)
self.numQuadPts = len(quadrature.get_weights())
- self._info.line("Basis (quad pts):")
- self._info.line(self.basis)
- self._info.line("Basis derivatives (quad pts):")
- self._info.line(self.basisDeriv)
+ self._info.line("Vertices: ")
+ self._info.line(self.vertices)
self._info.line("Quad pts:")
self._info.line(quadrature.get_points())
self._info.line("Quad wts:")
self._info.line(quadrature.get_weights())
+ self._info.line("Basis fns @ quad pts ):")
+ self._info.line(self.basis)
+ self._info.line("Basis fn derivatives @ quad pts:")
+ self._info.line(self.basisDeriv)
self._info.log()
return
Modified: short/3D/PyLith/trunk/unittests/pytests/feassemble/TestFIATLagrange.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/feassemble/TestFIATLagrange.py 2007-06-06 18:04:51 UTC (rev 7079)
+++ short/3D/PyLith/trunk/unittests/pytests/feassemble/TestFIATLagrange.py 2007-06-06 19:51:08 UTC (rev 7080)
@@ -16,141 +16,340 @@
import unittest
import numpy
+from pylith.feassemble.FIATLagrange import FIATLagrange
from pylith.utils.testarray import test_double
# ----------------------------------------------------------------------
-def N0(p):
- return (1-p[0])*(1-p[1])/4.0
+class Line2(object):
-def N0p(p):
- return -(1-p[1])/4.0
+ def __init__(self):
+ """
+ Setup line2 cell.
+ """
+ vertices = numpy.array([[-1.0], [1.0]])
+ quadPts = numpy.array( [[0.0]],
+ dtype=numpy.float64 )
+ quadWts = numpy.array( [2.0], dtype=numpy.float64 )
-def N0q(p):
- return -(1-p[0])/4.0
+ # Compute basis fns and derivatives at quadrature points
+ basis = numpy.zeros( (1, 2), dtype=numpy.float64)
+ basisDeriv = numpy.zeros( (1, 2, 1), dtype=numpy.float64)
+ iQuad = 0
+ for q in quadPts:
+ basis[iQuad] = numpy.array([self.N0(q), self.N1(q)],
+ dtype=numpy.float64).reshape( (2,) )
+ deriv = numpy.array([[self.N0p(q)], [self.N1p(q)]],
+ dtype=numpy.float64)
+ basisDeriv[iQuad] = deriv.reshape((2, 1))
+ iQuad += 1
-def N1(p):
- return (1+p[0])*(1-p[1])/4.0
+ self.cellDim = 1
+ self.numCorners = len(vertices)
+ self.numQuadPts = len(quadPts)
+ self.vertices = vertices
+ self.quadPts = quadPts
+ self.quadWts = quadWts
+ self.basis = basis
+ self.basisDeriv = basisDeriv
+ return
-def N1p(p):
- return (1-p[1])/4.0
-def N1q(p):
- return -(1+p[0])/4.0
+ def N0(self, p):
+ return 0.5*(1.0-p)
-def N2(p):
- return (1-p[0])*(1+p[1])/4.0
+ def N0p(self, p):
+ return -0.5
-def N2p(p):
- return -(1+p[1])/4.0
+ def N1(self, p):
+ return 0.5*(1.0+p)
-def N2q(p):
- return (1-p[0])/4.0
+ def N1p(self, p):
+ return 0.5
-def N3(p):
- return (1+p[0])*(1+p[1])/4.0
+# ----------------------------------------------------------------------
+class Line3(object):
-def N3p(p):
- return (1+p[1])/4.0
+ def __init__(self):
+ """
+ Setup line3 cell.
+ """
+ vertices = numpy.array([[-1.0], [1.0], [0.0]])
+ quadPts = numpy.array([ [-1.0/3**0.5],
+ [+1.0/3**0.5] ])
+ quadWts = numpy.array( [1.0, 1.0])
-def N3q(p):
- return (1+p[0])/4.0
+ # Compute basis fns and derivatives at quadrature points
+ basis = numpy.zeros( (2, 3), dtype=numpy.float64)
+ basisDeriv = numpy.zeros( (2, 3, 1), dtype=numpy.float64)
+ iQuad = 0
+ for q in quadPts:
+ basis[iQuad] = numpy.array([self.N0(q), self.N1(q), self.N2(q)],
+ dtype=numpy.float64).reshape( (3,) )
+ deriv = numpy.array([[self.N0p(q)], [self.N1p(q)], [self.N2p(q)]])
+ basisDeriv[iQuad] = deriv.reshape((3, 1))
+ iQuad += 1
-def Q0(p):
- return (1-p[0])*(1-p[1])*(1-p[2])/8.0
+ self.cellDim = 1
+ self.numCorners = len(vertices)
+ self.numQuadPts = len(quadPts)
+ self.vertices = vertices
+ self.quadPts = quadPts
+ self.quadWts = quadWts
+ self.basis = basis
+ self.basisDeriv = basisDeriv
+ return
-def Q0p(p):
- return -(1-p[1])*(1-p[2])/8.0
-def Q0q(p):
- return -(1-p[0])*(1-p[2])/8.0
+ def N0(self, p):
+ return -0.5*p*(1.0-p)
-def Q0r(p):
- return -(1-p[0])*(1-p[1])/8.0
+ def N0p(self, p):
+ return 1.0*p - 0.5
-def Q1(p):
- return (1+p[0])*(1-p[1])*(1-p[2])/8.0
+ def N1(self, p):
+ return 0.5*p*(1.0+p)
-def Q1p(p):
- return (1-p[1])*(1-p[2])/8.0
+ def N1p(self, p):
+ return 1.0*p + 0.5
-def Q1q(p):
- return -(1+p[0])*(1-p[2])/8.0
+ def N2(self, p):
+ return (1.0-p*p)
-def Q1r(p):
- return -(1+p[0])*(1-p[1])/8.0
+ def N2p(self, p):
+ return -2.0*p
-def Q2(p):
- return (1-p[0])*(1+p[1])*(1-p[2])/8.0
+# ----------------------------------------------------------------------
+class Quad4(object):
-def Q2p(p):
- return -(1+p[1])*(1-p[2])/8.0
+ def __init__(self):
+ """
+ Setup quad4 cell.
+ """
+ vertices = numpy.array([[-1.0, -1.0],
+ [+1.0, -1.0],
+ [-1.0, +1.0],
+ [+1.0, +1.0]])
+ quadPts = numpy.array([ [-1.0/3**0.5, -1.0/3**0.5],
+ [+1.0/3**0.5, -1.0/3**0.5],
+ [-1.0/3**0.5, +1.0/3**0.5],
+ [+1.0/3**0.5, +1.0/3**0.5] ])
+ quadWts = numpy.array( [1.0, 1.0, 1.0, 1.0])
-def Q2q(p):
- return (1-p[0])*(1-p[2])/8.0
+ # Compute basis fns and derivatives at quadrature points
+ basis = numpy.zeros( (4, 4), dtype=numpy.float64)
+ basisDeriv = numpy.zeros( (4, 4, 2), dtype=numpy.float64)
+ iQuad = 0
+ for q in quadPts:
+ basis[iQuad] = numpy.array([self.N0(q), self.N1(q),
+ self.N2(q), self.N3(q)],
+ dtype=numpy.float64).reshape( (4,) )
+ deriv = numpy.array([[self.N0p(q), self.N0q(q)],
+ [self.N1p(q), self.N1q(q)],
+ [self.N2p(q), self.N2q(q)],
+ [self.N3p(q), self.N3q(q)]])
+ basisDeriv[iQuad] = deriv.reshape((4, 2))
+ iQuad += 1
-def Q2r(p):
- return -(1-p[0])*(1+p[1])/8.0
+ self.cellDim = 2
+ self.numCorners = len(vertices)
+ self.numQuadPts = len(quadPts)
+ self.vertices = vertices
+ self.quadPts = quadPts
+ self.quadWts = quadWts
+ self.basis = basis
+ self.basisDeriv = basisDeriv
+ return
-def Q3(p):
- return (1+p[0])*(1+p[1])*(1-p[2])/8.0
-def Q3p(p):
- return (1+p[1])*(1-p[2])/8.0
+ def N0(self, p):
+ return (1-p[0])*(1-p[1])/4.0
-def Q3q(p):
- return (1+p[0])*(1-p[2])/8.0
+ def N0p(self, p):
+ return -(1-p[1])/4.0
-def Q3r(p):
- return -(1+p[0])*(1+p[1])/8.0
+ def N0q(self, p):
+ return -(1-p[0])/4.0
-def Q4(p):
- return (1-p[0])*(1-p[1])*(1+p[2])/8.0
+ def N1(self, p):
+ return (1+p[0])*(1-p[1])/4.0
-def Q4p(p):
- return -(1-p[1])*(1+p[2])/8.0
+ def N1p(self, p):
+ return (1-p[1])/4.0
-def Q4q(p):
- return -(1-p[0])*(1+p[2])/8.0
+ def N1q(self, p):
+ return -(1+p[0])/4.0
-def Q4r(p):
- return (1-p[0])*(1-p[1])/8.0
+ def N2(self, p):
+ return (1-p[0])*(1+p[1])/4.0
-def Q5(p):
- return (1+p[0])*(1-p[1])*(1+p[2])/8.0
+ def N2p(self, p):
+ return -(1+p[1])/4.0
-def Q5p(p):
- return (1-p[1])*(1+p[2])/8.0
+ def N2q(self, p):
+ return (1-p[0])/4.0
-def Q5q(p):
- return -(1+p[0])*(1+p[2])/8.0
+ def N3(self, p):
+ return (1+p[0])*(1+p[1])/4.0
-def Q5r(p):
- return (1+p[0])*(1-p[1])/8.0
+ def N3p(self, p):
+ return (1+p[1])/4.0
-def Q6(p):
- return (1-p[0])*(1+p[1])*(1+p[2])/8.0
+ def N3q(self, p):
+ return (1+p[0])/4.0
-def Q6p(p):
- return -(1+p[1])*(1+p[2])/8.0
+# ----------------------------------------------------------------------
+class Hex8(object):
-def Q6q(p):
- return (1-p[0])*(1+p[2])/8.0
+ def __init__(self):
+ """
+ Setup hex8 cell.
+ """
+ vertices = numpy.array([[-1.0, -1.0, -1.0],
+ [+1.0, -1.0, -1.0],
+ [-1.0, +1.0, -1.0],
+ [+1.0, +1.0, -1.0],
+ [-1.0, -1.0, +1.0],
+ [+1.0, -1.0, +1.0],
+ [-1.0, +1.0, +1.0],
+ [+1.0, +1.0, +1.0]])
+ quadPts = numpy.array([ [-1.0/3**0.5, -1.0/3**0.5, -1.0/3**0.5],
+ [+1.0/3**0.5, -1.0/3**0.5, -1.0/3**0.5],
+ [-1.0/3**0.5, +1.0/3**0.5, -1.0/3**0.5],
+ [+1.0/3**0.5, +1.0/3**0.5, -1.0/3**0.5],
+ [-1.0/3**0.5, -1.0/3**0.5, +1.0/3**0.5],
+ [+1.0/3**0.5, -1.0/3**0.5, +1.0/3**0.5],
+ [-1.0/3**0.5, +1.0/3**0.5, +1.0/3**0.5],
+ [+1.0/3**0.5, +1.0/3**0.5, +1.0/3**0.5]])
+ quadWts = numpy.array( [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0])
-def Q6r(p):
- return (1-p[0])*(1+p[1])/8.0
+ # Compute basis fns and derivatives at quadrature points
+ basis = numpy.zeros( (8, 8), dtype=numpy.float64)
+ basisDeriv = numpy.zeros( (8, 8, 3), dtype=numpy.float64)
+ iQuad = 0
+ for q in quadPts:
+ basis[iQuad] = numpy.array([self.N0(q), self.N1(q),
+ self.N2(q), self.N3(q),
+ self.N4(q), self.N5(q),
+ self.N6(q), self.N7(q)],
+ dtype=numpy.float64).reshape( (8,) )
+ deriv = numpy.array([[self.N0p(q), self.N0q(q), self.N0r(q)],
+ [self.N1p(q), self.N1q(q), self.N1r(q)],
+ [self.N2p(q), self.N2q(q), self.N2r(q)],
+ [self.N3p(q), self.N3q(q), self.N3r(q)],
+ [self.N4p(q), self.N4q(q), self.N4r(q)],
+ [self.N5p(q), self.N5q(q), self.N5r(q)],
+ [self.N6p(q), self.N6q(q), self.N6r(q)],
+ [self.N7p(q), self.N7q(q), self.N7r(q)]])
+ basisDeriv[iQuad] = deriv.reshape((8, 3))
+ iQuad += 1
-def Q7(p):
- return (1+p[0])*(1+p[1])*(1+p[2])/8.0
+ self.cellDim = 3
+ self.numCorners = len(vertices)
+ self.numQuadPts = len(quadPts)
+ self.vertices = vertices
+ self.quadPts = quadPts
+ self.quadWts = quadWts
+ self.basis = basis
+ self.basisDeriv = basisDeriv
+ return
-def Q7p(p):
- return (1+p[1])*(1+p[2])/8.0
-def Q7q(p):
- return (1+p[0])*(1+p[2])/8.0
+ def N0(self, p):
+ return (1-p[0])*(1-p[1])*(1-p[2])/8.0
+
+ def N0p(self, p):
+ return -(1-p[1])*(1-p[2])/8.0
+
+ def N0q(self, p):
+ return -(1-p[0])*(1-p[2])/8.0
+
+ def N0r(self, p):
+ return -(1-p[0])*(1-p[1])/8.0
+
+ def N1(self, p):
+ return (1+p[0])*(1-p[1])*(1-p[2])/8.0
+
+ def N1p(self, p):
+ return (1-p[1])*(1-p[2])/8.0
+
+ def N1q(self, p):
+ return -(1+p[0])*(1-p[2])/8.0
+
+ def N1r(self, p):
+ return -(1+p[0])*(1-p[1])/8.0
+
+ def N2(self, p):
+ return (1-p[0])*(1+p[1])*(1-p[2])/8.0
+
+ def N2p(self, p):
+ return -(1+p[1])*(1-p[2])/8.0
+
+ def N2q(self, p):
+ return (1-p[0])*(1-p[2])/8.0
+
+ def N2r(self, p):
+ return -(1-p[0])*(1+p[1])/8.0
+
+ def N3(self, p):
+ return (1+p[0])*(1+p[1])*(1-p[2])/8.0
+
+ def N3p(self, p):
+ return (1+p[1])*(1-p[2])/8.0
+
+ def N3q(self, p):
+ return (1+p[0])*(1-p[2])/8.0
-def Q7r(p):
- return (1+p[0])*(1+p[1])/8.0
+ def N3r(self, p):
+ return -(1+p[0])*(1+p[1])/8.0
+ def N4(self, p):
+ return (1-p[0])*(1-p[1])*(1+p[2])/8.0
+
+ def N4p(self, p):
+ return -(1-p[1])*(1+p[2])/8.0
+
+ def N4q(self, p):
+ return -(1-p[0])*(1+p[2])/8.0
+
+ def N4r(self, p):
+ return (1-p[0])*(1-p[1])/8.0
+
+ def N5(self, p):
+ return (1+p[0])*(1-p[1])*(1+p[2])/8.0
+
+ def N5p(self, p):
+ return (1-p[1])*(1+p[2])/8.0
+
+ def N5q(self, p):
+ return -(1+p[0])*(1+p[2])/8.0
+
+ def N5r(self, p):
+ return (1+p[0])*(1-p[1])/8.0
+
+ def N6(self, p):
+ return (1-p[0])*(1+p[1])*(1+p[2])/8.0
+
+ def N6p(self, p):
+ return -(1+p[1])*(1+p[2])/8.0
+
+ def N6q(self, p):
+ return (1-p[0])*(1+p[2])/8.0
+
+ def N6r(self, p):
+ return (1-p[0])*(1+p[1])/8.0
+
+ def N7(self, p):
+ return (1+p[0])*(1+p[1])*(1+p[2])/8.0
+
+ def N7p(self, p):
+ return (1+p[1])*(1+p[2])/8.0
+
+ def N7q(self, p):
+ return (1+p[0])*(1+p[2])/8.0
+
+ def N7r(self, p):
+ return (1+p[0])*(1+p[1])/8.0
+
# ----------------------------------------------------------------------
class TestFIATLagrange(unittest.TestCase):
"""
@@ -158,106 +357,87 @@
"""
- def test_initialize_quad(self):
+ def test_initialize_line2(self):
"""
- Test initialize().
+ Test initialize() with line2 cell.
"""
- from pylith.feassemble.FIATLagrange import FIATLagrange
+ from pylith.feassemble.FIATSimplex import FIATSimplex
cell = FIATLagrange()
- cell.cellDim = 2
+ cell.cellDim = 1
cell.degree = 1
- cell.order = 2
- quadPtsE = numpy.array( [(-1.0/3**0.5, -1.0/3**0.5),
- (+1.0/3**0.5, -1.0/3**0.5),
- (-1.0/3**0.5, +1.0/3**0.5),
- (+1.0/3**0.5, +1.0/3**0.5)],
- dtype=numpy.float64 )
- quadWtsE = numpy.array( [1.0, 1.0, 1.0, 1.0], dtype=numpy.float64 )
+ cell.order = 1
+ cell.initialize()
- # Compute basis fns and derivatives at quadrature points
- basis = numpy.zeros( (4, 4), dtype=numpy.float64)
- basisDeriv = numpy.zeros( (4, 4, 2), dtype=numpy.float64)
- iQuad = 0
- for q in quadPtsE:
- basis[iQuad] = numpy.array([N0(q), N1(q), N2(q), N3(q)],
- dtype=numpy.float64).reshape( (4,) )
- deriv = numpy.array([[N0p(q), N0q(q)],
- [N1p(q), N1q(q)],
- [N2p(q), N2q(q)],
- [N3p(q), N3q(q)]],
- dtype=numpy.float64)
- basisDeriv[iQuad] = deriv.reshape((4, 2))
- iQuad += 1
+ cellE = Line2()
+ self._checkVals(cellE, cell)
+ return
+
+ def test_initialize_line3(self):
+ """
+ Test initialize() with line3 cell.
+ """
+ from pylith.feassemble.FIATSimplex import FIATSimplex
+ cell = FIATLagrange()
+ cell.cellDim = 1
+ cell.degree = 2
+ cell.order = 2
cell.initialize()
- # Check basic attributes
- self.assertEqual(2, cell.cellDim)
- self.assertEqual(4, cell.numCorners)
- self.assertEqual(4, cell.numQuadPts)
+ cellE = Line3()
+ self._checkVals(cellE, cell)
+ return
- # Check arrays
- test_double(self, quadPtsE, cell.quadPts)
- test_double(self, quadWtsE, cell.quadWts)
- test_double(self, basis, cell.basis)
- test_double(self, basisDeriv, cell.basisDeriv)
+ def test_initialize_quad4(self):
+ """
+ Test initialize() with quad4 cell.
+ """
+ from pylith.feassemble.FIATSimplex import FIATSimplex
+ cell = FIATLagrange()
+ cell.cellDim = 2
+ cell.degree = 1
+ cell.order = 2
+ cell.initialize()
+
+ cellE = Quad4()
+ self._checkVals(cellE, cell)
return
- def test_initialize_hex(self):
+ def test_initialize_hex8(self):
"""
- Test initialize().
+ Test initialize() with hex8 cell.
"""
- from pylith.feassemble.FIATLagrange import FIATLagrange
+ from pylith.feassemble.FIATSimplex import FIATSimplex
cell = FIATLagrange()
cell.cellDim = 3
cell.degree = 1
cell.order = 2
- quadPtsE = numpy.array( [(-1.0/3**0.5, -1.0/3**0.5, -1.0/3**0.5),
- (+1.0/3**0.5, -1.0/3**0.5, -1.0/3**0.5),
- (-1.0/3**0.5, +1.0/3**0.5, -1.0/3**0.5),
- (+1.0/3**0.5, +1.0/3**0.5, -1.0/3**0.5),
- (-1.0/3**0.5, -1.0/3**0.5, +1.0/3**0.5),
- (+1.0/3**0.5, -1.0/3**0.5, +1.0/3**0.5),
- (-1.0/3**0.5, +1.0/3**0.5, +1.0/3**0.5),
- (+1.0/3**0.5, +1.0/3**0.5, +1.0/3**0.5)],
- dtype=numpy.float64 )
- quadWtsE = numpy.array( [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], dtype=numpy.float64 )
+ cell.initialize()
- # Compute basis fns and derivatives at quadrature points
- iQuad = 0
- basis = numpy.zeros( (8, 8), dtype=numpy.float64)
- basisDeriv = numpy.zeros( (8, 8, 3), dtype=numpy.float64)
- for q in quadPtsE:
- basis[iQuad] = numpy.array([Q0(q), Q1(q), Q2(q), Q3(q), Q4(q), Q5(q), Q6(q), Q7(q)],
- dtype=numpy.float64).reshape( (8,) )
- deriv = numpy.array([[Q0p(q), Q0q(q), Q0r(q)],
- [Q1p(q), Q1q(q), Q1r(q)],
- [Q2p(q), Q2q(q), Q2r(q)],
- [Q3p(q), Q3q(q), Q3r(q)],
- [Q4p(q), Q4q(q), Q4r(q)],
- [Q5p(q), Q5q(q), Q5r(q)],
- [Q6p(q), Q6q(q), Q6r(q)],
- [Q7p(q), Q7q(q), Q7r(q)]],
- dtype=numpy.float64)
- basisDeriv[iQuad] = deriv.reshape((8, 3))
- iQuad += 1
+ cellE = Hex8()
+ self._checkVals(cellE, cell)
+ return
- cell.initialize()
+ def _checkVals(self, cellE, cell):
+ """
+ Check known values against those generated by FIATSimplex.
+ """
+
# Check basic attributes
- self.assertEqual(3, cell.cellDim)
- self.assertEqual(8, cell.numCorners)
- self.assertEqual(8, cell.numQuadPts)
+ self.assertEqual(cellE.cellDim, cell.cellDim)
+ self.assertEqual(cellE.numCorners, cell.numCorners)
+ self.assertEqual(cellE.numQuadPts, cell.numQuadPts)
# Check arrays
- test_double(self, quadPtsE, cell.quadPts)
- test_double(self, quadWtsE, cell.quadWts)
- test_double(self, basis, cell.basis)
- test_double(self, basisDeriv, cell.basisDeriv)
-
+ test_double(self, cellE.vertices, cell.vertices)
+ test_double(self, cellE.quadPts, cell.quadPts)
+ test_double(self, cellE.quadWts, cell.quadWts)
+ test_double(self, cellE.basis, cell.basis)
+ test_double(self, cellE.basisDeriv, cell.basisDeriv)
return
-# End of file
+# End of file
Modified: short/3D/PyLith/trunk/unittests/pytests/feassemble/TestFIATSimplex.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/feassemble/TestFIATSimplex.py 2007-06-06 18:04:51 UTC (rev 7079)
+++ short/3D/PyLith/trunk/unittests/pytests/feassemble/TestFIATSimplex.py 2007-06-06 19:51:08 UTC (rev 7080)
@@ -16,32 +16,261 @@
import unittest
import numpy
+from pylith.feassemble.FIATSimplex import FIATSimplex
from pylith.utils.testarray import test_double
# ----------------------------------------------------------------------
-def N0(p):
- return -0.5*p*(1.0-p)
+class Line2(object):
-def N0p(p):
- return 1.0*p - 0.5
+ def __init__(self):
+ """
+ Setup line2 cell.
+ """
+ vertices = numpy.array([[-1.0], [1.0]])
+ quadPts = numpy.array( [[0.0]],
+ dtype=numpy.float64 )
+ quadWts = numpy.array( [2.0], dtype=numpy.float64 )
-def N1(p):
- return 0.5*p*(1.0+p)
+ # Compute basis fns and derivatives at quadrature points
+ basis = numpy.zeros( (1, 2), dtype=numpy.float64)
+ basisDeriv = numpy.zeros( (1, 2, 1), dtype=numpy.float64)
+ iQuad = 0
+ for q in quadPts:
+ basis[iQuad] = numpy.array([self.N0(q), self.N1(q)],
+ dtype=numpy.float64).reshape( (2,) )
+ deriv = numpy.array([[self.N0p(q)], [self.N1p(q)]],
+ dtype=numpy.float64)
+ basisDeriv[iQuad] = deriv.reshape((2, 1))
+ iQuad += 1
-def N1p(p):
- return 1.0*p + 0.5
+ self.cellDim = 1
+ self.numCorners = len(vertices)
+ self.numQuadPts = len(quadPts)
+ self.vertices = vertices
+ self.quadPts = quadPts
+ self.quadWts = quadWts
+ self.basis = basis
+ self.basisDeriv = basisDeriv
+ return
-def N2(p):
- return (1.0-p*p)
-def N2p(p):
- return -2.0*p
+ def N0(self, p):
+ return 0.5*(1.0-p)
-def nodesRef():
- return [-1.0, 1.0, 0.0]
+ def N0p(self, p):
+ return -0.5
+ def N1(self, p):
+ return 0.5*(1.0+p)
+ def N1p(self, p):
+ return 0.5
+
# ----------------------------------------------------------------------
+class Line3(object):
+
+ def __init__(self):
+ """
+ Setup line3 cell.
+ """
+ vertices = numpy.array([[-1.0], [1.0], [0.0]])
+ quadPts = numpy.array([ [-1.0/3**0.5],
+ [+1.0/3**0.5] ])
+ quadWts = numpy.array( [1.0, 1.0])
+
+ # Compute basis fns and derivatives at quadrature points
+ basis = numpy.zeros( (2, 3), dtype=numpy.float64)
+ basisDeriv = numpy.zeros( (2, 3, 1), dtype=numpy.float64)
+ iQuad = 0
+ for q in quadPts:
+ basis[iQuad] = numpy.array([self.N0(q), self.N1(q), self.N2(q)],
+ dtype=numpy.float64).reshape( (3,) )
+ deriv = numpy.array([[self.N0p(q)], [self.N1p(q)], [self.N2p(q)]])
+ basisDeriv[iQuad] = deriv.reshape((3, 1))
+ iQuad += 1
+
+ self.cellDim = 1
+ self.numCorners = len(vertices)
+ self.numQuadPts = len(quadPts)
+ self.vertices = vertices
+ self.quadPts = quadPts
+ self.quadWts = quadWts
+ self.basis = basis
+ self.basisDeriv = basisDeriv
+ return
+
+
+ def N0(self, p):
+ return -0.5*p*(1.0-p)
+
+ def N0p(self, p):
+ return 1.0*p - 0.5
+
+ def N1(self, p):
+ return 0.5*p*(1.0+p)
+
+ def N1p(self, p):
+ return 1.0*p + 0.5
+
+ def N2(self, p):
+ return (1.0-p*p)
+
+ def N2p(self, p):
+ return -2.0*p
+
+# ----------------------------------------------------------------------
+class Tri3(object):
+
+ def __init__(self):
+ """
+ Setup tri33 cell.
+ """
+ vertices = numpy.array([[-1.0, -1.0],
+ [+1.0, -1.0],
+ [-1.0, +1.0]])
+ quadPts = numpy.array([ [-1.0/3.0, -1.0/3.0] ])
+ quadWts = numpy.array( [2.0])
+
+ # Compute basis fns and derivatives at quadrature points
+ basis = numpy.zeros( (1, 3), dtype=numpy.float64)
+ basisDeriv = numpy.zeros( (1, 3, 2), dtype=numpy.float64)
+ iQuad = 0
+ for q in quadPts:
+ basis[iQuad] = numpy.array([self.N0(q), self.N1(q), self.N2(q)],
+ dtype=numpy.float64).reshape( (3,) )
+ deriv = numpy.array([[self.N0p(q), self.N0q(q)],
+ [self.N1p(q), self.N1q(q)],
+ [self.N2p(q), self.N2q(q)]])
+ basisDeriv[iQuad] = deriv.reshape((3, 2))
+ iQuad += 1
+
+ self.cellDim = 2
+ self.numCorners = len(vertices)
+ self.numQuadPts = len(quadPts)
+ self.vertices = vertices
+ self.quadPts = quadPts
+ self.quadWts = quadWts
+ self.basis = basis
+ self.basisDeriv = basisDeriv
+ return
+
+
+ def N0(self, p):
+ return 0.5*(-p[0]-p[1])
+
+ def N0p(self, p):
+ return -0.5
+
+ def N0q(self, p):
+ return -0.5
+
+ def N1(self, p):
+ return 0.5*(1.0+p[0])
+
+ def N1p(self, p):
+ return 0.5
+
+ def N1q(self, p):
+ return 0.0
+
+ def N2(self, p):
+ return 0.5*(1.0+p[1])
+
+ def N2p(self, p):
+ return 0.0
+
+ def N2q(self, p):
+ return 0.5
+
+# ----------------------------------------------------------------------
+class Tet4(object):
+
+ def __init__(self):
+ """
+ Setup tri33 cell.
+ """
+ vertices = numpy.array([[-1.0, -1.0, -1.0],
+ [+1.0, -1.0, -1.0],
+ [-1.0, +1.0, -1.0],
+ [-1.0, -1.0, +1.0]])
+ quadPts = numpy.array([ [-1.0/2.0, -1.0/2.0, -1.0/2.0] ])
+ quadWts = numpy.array( [4.0/3.0])
+
+ # Compute basis fns and derivatives at quadrature points
+ basis = numpy.zeros( (1, 4), dtype=numpy.float64)
+ basisDeriv = numpy.zeros( (1, 4, 3), dtype=numpy.float64)
+ iQuad = 0
+ for q in quadPts:
+ basis[iQuad] = numpy.array([self.N0(q), self.N1(q),
+ self.N2(q), self.N3(q)],
+ dtype=numpy.float64).reshape( (4,) )
+ deriv = numpy.array([[self.N0p(q), self.N0q(q), self.N0r(q)],
+ [self.N1p(q), self.N1q(q), self.N1r(q)],
+ [self.N2p(q), self.N2q(q), self.N2r(q)],
+ [self.N3p(q), self.N3q(q), self.N3r(q)]])
+ basisDeriv[iQuad] = deriv.reshape((4, 3))
+ iQuad += 1
+
+ self.cellDim = 3
+ self.numCorners = len(vertices)
+ self.numQuadPts = len(quadPts)
+ self.vertices = vertices
+ self.quadPts = quadPts
+ self.quadWts = quadWts
+ self.basis = basis
+ self.basisDeriv = basisDeriv
+ return
+
+
+ def N0(self, p):
+ return 0.5*(-1-p[0]-p[1]-p[2])
+
+ def N0p(self, p):
+ return -0.5
+
+ def N0q(self, p):
+ return -0.5
+
+ def N0r(self, p):
+ return -0.5
+
+ def N1(self, p):
+ return 0.5*(1.0+p[0])
+
+ def N1p(self, p):
+ return 0.5
+
+ def N1q(self, p):
+ return 0.0
+
+ def N1r(self, p):
+ return 0.0
+
+ def N2(self, p):
+ return 0.5*(1.0+p[1])
+
+ def N2p(self, p):
+ return 0.0
+
+ def N2q(self, p):
+ return 0.5
+
+ def N2r(self, p):
+ return 0.0
+
+ def N3(self, p):
+ return 0.5*(1.0+p[2])
+
+ def N3p(self, p):
+ return 0.0
+
+ def N3q(self, p):
+ return 0.0
+
+ def N3r(self, p):
+ return 0.5
+
+# ----------------------------------------------------------------------
class TestFIATSimplex(unittest.TestCase):
"""
Unit testing of FIATSimplex object.
@@ -52,7 +281,6 @@
"""
Test _getShape().
"""
- from pylith.feassemble.FIATSimplex import FIATSimplex
cell = FIATSimplex()
import FIAT.shapes
@@ -71,45 +299,82 @@
return
- def test_initialize(self):
+ def test_initialize_line2(self):
"""
- Test initialize().
+ Test initialize() with line2 cell.
"""
- from pylith.feassemble.FIATSimplex import FIATSimplex
cell = FIATSimplex()
cell.shape = "line"
+ cell.degree = 1
+ cell.order = 1
+ cell.initialize()
+
+ cellE = Line2()
+ self._checkVals(cellE, cell)
+ return
+
+
+ def test_initialize_line3(self):
+ """
+ Test initialize() with line3 cell.
+ """
+ cell = FIATSimplex()
+ cell.shape = "line"
cell.degree = 2
cell.order = 2
- 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 )
+ cell.initialize()
- # Compute basis fns and derivatives at quadrature points
- basis = numpy.zeros( (2, 3), dtype=numpy.float64)
- basisDeriv = numpy.zeros( (2, 3, 1), dtype=numpy.float64)
- iQuad = 0
- for q in quadPtsE:
- basis[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)
- basisDeriv[iQuad] = deriv.reshape((3, 1))
- iQuad += 1
+ cellE = Line3()
+ self._checkVals(cellE, cell)
+ return
+
+ def test_initialize_tri3(self):
+ """
+ Test initialize() with tri3 cell.
+ """
+ cell = FIATSimplex()
+ cell.shape = "triangle"
+ cell.degree = 1
+ cell.order = 1
cell.initialize()
+ cellE = Tri3()
+ self._checkVals(cellE, cell)
+ return
+
+
+ def test_initialize_tet4(self):
+ """
+ Test initialize() with tet4 cell.
+ """
+ cell = FIATSimplex()
+ cell.shape = "tetrahedron"
+ cell.degree = 1
+ cell.order = 1
+ cell.initialize()
+
+ cellE = Tet4()
+ self._checkVals(cellE, cell)
+ return
+
+
+ def _checkVals(self, cellE, cell):
+ """
+ Check known values against those generated by FIATSimplex.
+ """
+
# Check basic attributes
- self.assertEqual(1, cell.cellDim)
- self.assertEqual(3, cell.numCorners)
- self.assertEqual(2, cell.numQuadPts)
+ self.assertEqual(cellE.cellDim, cell.cellDim)
+ self.assertEqual(cellE.numCorners, cell.numCorners)
+ self.assertEqual(cellE.numQuadPts, cell.numQuadPts)
# Check arrays
- test_double(self, quadPtsE, cell.quadPts)
- test_double(self, quadWtsE, cell.quadWts)
- test_double(self, basis, cell.basis)
- test_double(self, basisDeriv, cell.basisDeriv)
-
+ test_double(self, cellE.vertices, cell.vertices)
+ test_double(self, cellE.quadPts, cell.quadPts)
+ test_double(self, cellE.quadWts, cell.quadWts)
+ test_double(self, cellE.basis, cell.basis)
+ test_double(self, cellE.basisDeriv, cell.basisDeriv)
return
More information about the cig-commits
mailing list