[cig-commits] r17128 - in short/3D/PyLith/branches/v1.5-stable: pylith/feassemble unittests/pytests/feassemble
brad at geodynamics.org
brad at geodynamics.org
Thu Aug 26 18:00:49 PDT 2010
Author: brad
Date: 2010-08-26 18:00:49 -0700 (Thu, 26 Aug 2010)
New Revision: 17128
Modified:
short/3D/PyLith/branches/v1.5-stable/pylith/feassemble/FIATLagrange.py
short/3D/PyLith/branches/v1.5-stable/unittests/pytests/feassemble/TestFIATLagrange.py
Log:
Started fixing orderings for higher order hex cells.
Modified: short/3D/PyLith/branches/v1.5-stable/pylith/feassemble/FIATLagrange.py
===================================================================
--- short/3D/PyLith/branches/v1.5-stable/pylith/feassemble/FIATLagrange.py 2010-08-27 01:00:41 UTC (rev 17127)
+++ short/3D/PyLith/branches/v1.5-stable/pylith/feassemble/FIATLagrange.py 2010-08-27 01:00:49 UTC (rev 17128)
@@ -155,9 +155,9 @@
n += 1
# Left
- for r in range(numBasisFns-1, 1, -1):
+ for q in range(numBasisFns-1, 1, -1):
self.vertices[n][0] = vertices[0]
- self.vertices[n][1] = vertices[r]
+ self.vertices[n][1] = vertices[q]
n += 1
# Face
@@ -168,7 +168,8 @@
n += 1
if not n == self.numCorners:
- raise RuntimeError('Invalid 2D function tabulation, n is %d' % n)
+ raise RuntimeError('Invalid 2-D function tabulation: '+str(n)+ \
+ ' should be '+str(self.numCorners))
self.quadPts = numpy.zeros((numQuadPts*numQuadPts, dim))
self.quadWts = numpy.zeros((numQuadPts*numQuadPts,))
@@ -177,7 +178,8 @@
self.basisDeriv = numpy.zeros((numQuadPts*numQuadPts,
numBasisFns*numBasisFns, dim))
- # Order of basis functions and quadrature points don't matter
+ # Order of quadrature points doesn't matter
+ # Order of basis functions should match vertices for isoparametric
n = 0
for q in range(0, numQuadPts):
for p in range(0, numQuadPts):
@@ -249,48 +251,189 @@
self.basisDeriv[n][m][1] = basis[p][bp]*basisDeriv[q][bq][0]
m += 1
- if not m == numBasisFns**2: raise RuntimeError('Invalid 2D quadrature')
+ if not m == numBasisFns**2:
+ raise RuntimeError('Invalid 2-D quadrature')
n += 1
- if not n == self.numQuadPts: raise RuntimeError('Invalid 2D quadrature')
+ if not n == self.numQuadPts:
+ raise RuntimeError('Invalid 2-D quadrature')
+
elif dim == 3:
self.vertices = numpy.zeros((self.numCorners, dim))
n = 0
- # Depth
- for s in range(numBasisFns):
- # Bottom
- for r in range(0, numBasisFns-1):
- self.vertices[n][0] = vertices[r]
- self.vertices[n][1] = vertices[0]
- self.vertices[n][2] = vertices[s]
+ # Corners
+ self.vertices[n][0] = vertices[0]
+ self.vertices[n][1] = vertices[0]
+ self.vertices[n][2] = vertices[0]
+ n += 1
+ self.vertices[n][0] = vertices[1]
+ self.vertices[n][1] = vertices[0]
+ self.vertices[n][2] = vertices[0]
+ n += 1
+ self.vertices[n][0] = vertices[1]
+ self.vertices[n][1] = vertices[1]
+ self.vertices[n][2] = vertices[0]
+ n += 1
+ self.vertices[n][0] = vertices[0]
+ self.vertices[n][1] = vertices[1]
+ self.vertices[n][2] = vertices[0]
+ n += 1
+ self.vertices[n][0] = vertices[0]
+ self.vertices[n][1] = vertices[0]
+ self.vertices[n][2] = vertices[1]
+ n += 1
+ self.vertices[n][0] = vertices[1]
+ self.vertices[n][1] = vertices[0]
+ self.vertices[n][2] = vertices[1]
+ n += 1
+ self.vertices[n][0] = vertices[1]
+ self.vertices[n][1] = vertices[1]
+ self.vertices[n][2] = vertices[1]
+ n += 1
+ self.vertices[n][0] = vertices[0]
+ self.vertices[n][1] = vertices[1]
+ self.vertices[n][2] = vertices[1]
+ n += 1
+
+ # Edges
+ # Bottom front
+ for p in range(2, numBasisFns):
+ self.vertices[n][0] = vertices[p]
+ self.vertices[n][1] = vertices[0]
+ self.vertices[n][2] = vertices[0]
+ n += 1
+ # Bottom right
+ for q in range(2, numBasisFns):
+ self.vertices[n][0] = vertices[1]
+ self.vertices[n][1] = vertices[q]
+ self.vertices[n][2] = vertices[0]
+ n += 1
+
+ # Bottom back
+ for p in range(numBasisFns-1, 1, -1):
+ self.vertices[n][0] = vertices[p]
+ self.vertices[n][1] = vertices[1]
+ self.vertices[n][2] = vertices[0]
+ n += 1
+
+ # Bottom left
+ for q in range(numBasisFns-1, 1, -1):
+ self.vertices[n][0] = vertices[0]
+ self.vertices[n][1] = vertices[q]
+ self.vertices[n][2] = vertices[0]
+ n += 1
+ # Middle left front
+ for r in range(2, numBasisFns):
+ self.vertices[n][0] = vertices[0]
+ self.vertices[n][1] = vertices[0]
+ self.vertices[n][2] = vertices[r]
+ n += 1
+ # Middle right front
+ for r in range(2, numBasisFns):
+ self.vertices[n][0] = vertices[1]
+ self.vertices[n][1] = vertices[0]
+ self.vertices[n][2] = vertices[r]
+ n += 1
+
+ # Middle right back
+ for r in range(2, numBasisFns):
+ self.vertices[n][0] = vertices[1]
+ self.vertices[n][1] = vertices[1]
+ self.vertices[n][2] = vertices[r]
+ n += 1
+
+ # Middle left back
+ for r in range(2, numBasisFns):
+ self.vertices[n][0] = vertices[0]
+ self.vertices[n][1] = vertices[1]
+ self.vertices[n][2] = vertices[r]
+ n += 1
+ # Top front
+ for p in range(2, numBasisFns):
+ self.vertices[n][0] = vertices[p]
+ self.vertices[n][1] = vertices[0]
+ self.vertices[n][2] = vertices[1]
+ n += 1
+ # Top right
+ for q in range(2, numBasisFns):
+ self.vertices[n][0] = vertices[1]
+ self.vertices[n][1] = vertices[q]
+ self.vertices[n][2] = vertices[1]
+ n += 1
+
+ # Top back
+ for p in range(numBasisFns-1, 1, -1):
+ self.vertices[n][0] = vertices[p]
+ self.vertices[n][1] = vertices[1]
+ self.vertices[n][2] = vertices[1]
+ n += 1
+
+ # Top left
+ for r in range(numBasisFns-1, 1, -1):
+ self.vertices[n][0] = vertices[0]
+ self.vertices[n][1] = vertices[r]
+ self.vertices[n][2] = vertices[1]
+ n += 1
+
+ # Face
+ # Interior
+ for r in range(2, numBasisFns):
+ for q in range(2, numBasisFns):
+ for p in range(2, numBasisFns):
+ self.vertices[n][0] = vertices[p]
+ self.vertices[n][1] = vertices[q]
+ self.vertices[n][2] = vertices[r]
+ n += 1
+
+ # Bottom
+ for q in range(2, numBasisFns):
+ for p in range(2, numBasisFns):
+ self.vertices[n][0] = vertices[p]
+ self.vertices[n][1] = vertices[q]
+ self.vertices[n][2] = vertices[0]
n += 1
- # Right
- for q in range(0, numBasisFns-1):
- self.vertices[n][0] = vertices[numBasisFns-1]
+ # Top
+ for q in range(2, numBasisFns):
+ for p in range(2, numBasisFns):
+ self.vertices[n][0] = vertices[p]
self.vertices[n][1] = vertices[q]
- self.vertices[n][2] = vertices[s]
+ self.vertices[n][2] = vertices[1]
n += 1
- # Top
- for r in range(numBasisFns-1, 0, -1):
- self.vertices[n][0] = vertices[r]
- self.vertices[n][1] = vertices[numBasisFns-1]
- self.vertices[n][2] = vertices[s]
- n += 1
- # Left
- for q in range(numBasisFns-1, 0, -1):
+ # Left
+ for r in range(2, numBasisFns):
+ for q in range(2, numBasisFns):
self.vertices[n][0] = vertices[0]
self.vertices[n][1] = vertices[q]
- self.vertices[n][2] = vertices[s]
+ self.vertices[n][2] = vertices[r]
n += 1
- # Interior
- for q in range(1, numBasisFns-1):
- for r in range(1, numBasisFns-1):
- self.vertices[n][0] = vertices[r]
- self.vertices[n][1] = vertices[q]
- self.vertices[n][2] = vertices[s]
- n += 1
- if not n == self.numCorners: raise RuntimeError('Invalid 3D function tabulation: '+str(n)+' should be '+str(self.numCorners))
+ # Right
+ for r in range(2, numBasisFns):
+ for q in range(2, numBasisFns):
+ self.vertices[n][0] = vertices[1]
+ self.vertices[n][1] = vertices[q]
+ self.vertices[n][2] = vertices[r]
+ n += 1
+ # Front
+ for r in range(2, numBasisFns):
+ for p in range(2, numBasisFns):
+ self.vertices[n][0] = vertices[p]
+ self.vertices[n][1] = vertices[0]
+ self.vertices[n][2] = vertices[r]
+ n += 1
+ # Back
+ for r in range(2, numBasisFns):
+ for p in range(2, numBasisFns):
+ self.vertices[n][0] = vertices[p]
+ self.vertices[n][1] = vertices[1]
+ self.vertices[n][2] = vertices[r]
+ n += 1
+ if not n == self.numCorners:
+ raise RuntimeError('Invalid 3-D function tabulation: '+str(n)+ \
+ ' should be '+str(self.numCorners))
+
+ print "VERTICES",self.vertices
+
self.quadPts = numpy.zeros((numQuadPts*numQuadPts*numQuadPts, dim))
self.quadWts = numpy.zeros((numQuadPts*numQuadPts*numQuadPts,))
self.basis = numpy.zeros((numQuadPts*numQuadPts*numQuadPts,
Modified: short/3D/PyLith/branches/v1.5-stable/unittests/pytests/feassemble/TestFIATLagrange.py
===================================================================
--- short/3D/PyLith/branches/v1.5-stable/unittests/pytests/feassemble/TestFIATLagrange.py 2010-08-27 01:00:41 UTC (rev 17127)
+++ short/3D/PyLith/branches/v1.5-stable/unittests/pytests/feassemble/TestFIATLagrange.py 2010-08-27 01:00:49 UTC (rev 17128)
@@ -505,6 +505,178 @@
return (1-p[0])*(1+p[1])/8.0
# ----------------------------------------------------------------------
+class Hex27(object):
+
+ def __init__(self):
+ """
+ Setup hex8 cell.
+ """
+ vertices = numpy.array([[-1.0, -1.0, -1.0], # Corners
+ [+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],
+ [ 0.0, -1.0, -1.0], # Bottom edges
+ [+1.0, 0.0, -1.0],
+ [ 0.0, +1.0, -1.0],
+ [-1.0, 0.0, -1.0],
+ [-1.0, -1.0, 0.0], # Middle edges
+ [+1.0, -1.0, 0.0],
+ [+1.0, +1.0, 0.0],
+ [-1.0, +1.0, 0.0],
+ [ 0.0, -1.0, +1.0], # Top edges
+ [+1.0, 0.0, +1.0],
+ [ 0.0, +1.0, +1.0],
+ [-1.0, 0.0, +1.0],
+ [ 0.0, 0.0, 0.0], # Interior
+ [ 0.0, 0.0, -1.0], # Faces
+ [ 0.0, 0.0, +1.0],
+ [-1.0, 0.0, 0.0],
+ [+1.0, 0.0, 0.0],
+ [ 0.0, -1.0, 0.0],
+ [ 0.0, +1.0, 0.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])
+
+ # 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
+
+ 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 (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 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):
"""
Unit testing of FIATLagrange object.
More information about the CIG-COMMITS
mailing list