[cig-commits] r17135 - in short/3D/PyLith/branches/v1.5-stable: pylith/feassemble unittests/pytests/feassemble
brad at geodynamics.org
brad at geodynamics.org
Fri Aug 27 14:02:28 PDT 2010
Author: brad
Date: 2010-08-27 14:02:28 -0700 (Fri, 27 Aug 2010)
New Revision: 17135
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:
Finished fixing ordering of vertices in FIATLagrange. Added unit test for hex27 cell.
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 16:52:31 UTC (rev 17134)
+++ short/3D/PyLith/branches/v1.5-stable/pylith/feassemble/FIATLagrange.py 2010-08-27 21:02:28 UTC (rev 17135)
@@ -103,14 +103,14 @@
quadwts = numpy.array(quadrature.get_weights())
numQuadPts = len(quadpts)
basis = numpy.array(element.function_space().tabulate(quadrature.get_points())).transpose()
- numBasisFns = len(element.function_space())
+ numBasis = len(element.function_space())
# Evaluate derivatives of basis functions at quadrature points
basisDeriv = numpy.array([element.function_space().deriv_all(d).tabulate(quadrature.get_points()) \
for d in range(1)]).transpose()
self.numQuadPts = numQuadPts**dim
- self.numCorners = numBasisFns**dim
+ self.numCorners = numBasis**dim
if dim == 1:
self.vertices = numpy.array(vertices)
@@ -120,63 +120,47 @@
self.basisDeriv = numpy.reshape(basisDeriv.flatten(), basisDeriv.shape)
else:
if dim == 2:
- self.vertices = numpy.zeros((self.numCorners, dim))
- n = 0
+ # Set order of vertices and basis functions.
# Corners
- self.vertices[n][0] = vertices[0]
- self.vertices[n][1] = vertices[0]
- n += 1
- self.vertices[n][0] = vertices[1]
- self.vertices[n][1] = vertices[0]
- n += 1
- self.vertices[n][0] = vertices[1]
- self.vertices[n][1] = vertices[1]
- n += 1
- self.vertices[n][0] = vertices[0]
- self.vertices[n][1] = vertices[1]
- n += 1
-
+ vertexOrder = [(0,0), (1,0), (1,1), (0,1)]
# Edges
# Bottom
- for p in range(2, numBasisFns):
- self.vertices[n][0] = vertices[p]
- self.vertices[n][1] = vertices[0]
- n += 1
+ p = numpy.arange(2, numBasis, dtype=numpy.int32)
+ q = numpy.zeros(numBasis-2, dtype=numpy.int32)
+ vertexOrder += zip(p,q)
# Right
- for q in range(2, numBasisFns):
- self.vertices[n][0] = vertices[1]
- self.vertices[n][1] = vertices[q]
- n += 1
-
+ p = numpy.ones(numBasis-2, dtype=numpy.int32)
+ q = numpy.arange(2, numBasis, dtype=numpy.int32)
+ vertexOrder += zip(p,q)
# Top
- for p in range(numBasisFns-1, 1, -1):
- self.vertices[n][0] = vertices[p]
- self.vertices[n][1] = vertices[1]
- n += 1
-
+ p = numpy.arange(numBasis-1, 1, step=-1, dtype=numpy.int32)
+ q = numpy.ones(numBasis-2, dtype=numpy.int32)
+ vertexOrder += zip(p,q)
# Left
- for q in range(numBasisFns-1, 1, -1):
- self.vertices[n][0] = vertices[0]
+ p = numpy.zeros(numBasis-2, dtype=numpy.int32)
+ q = numpy.arange(numBasis-1, 1, step=-1, dtype=numpy.int32)
+ vertexOrder += zip(p,q)
+ # Face
+ p = numpy.arange(2, numBasis, dtype=numpy.int32)
+ q = numpy.arange(2, numBasis, dtype=numpy.int32)
+ vertexOrder += zip(p,q)
+
+ self.vertices = numpy.zeros((self.numCorners, dim))
+ n = 0
+ for (p,q) in vertexOrder:
+ self.vertices[n][0] = vertices[p]
self.vertices[n][1] = vertices[q]
n += 1
-
- # Face
- for q in range(2, numBasisFns):
- for p in range(2, numBasisFns):
- self.vertices[n][0] = vertices[p]
- self.vertices[n][1] = vertices[q]
- n += 1
-
if not n == self.numCorners:
- raise RuntimeError('Invalid 2-D function tabulation: '+str(n)+ \
- ' should be '+str(self.numCorners))
+ raise RuntimeError('Invalid 2-D vertex ordering: '+str(n)+ \
+ ' should be '+str(self.numCorners))
self.quadPts = numpy.zeros((numQuadPts*numQuadPts, dim))
self.quadWts = numpy.zeros((numQuadPts*numQuadPts,))
self.basis = numpy.zeros((numQuadPts*numQuadPts,
- numBasisFns*numBasisFns))
+ numBasis*numBasis))
self.basisDeriv = numpy.zeros((numQuadPts*numQuadPts,
- numBasisFns*numBasisFns, dim))
+ numBasis*numBasis, dim))
# Order of quadrature points doesn't matter
# Order of basis functions should match vertices for isoparametric
@@ -188,494 +172,168 @@
self.quadWts[n] = quadwts[p]*quadwts[q]
m = 0
- # Corners
- bp = 0; bq = 0
- self.basis[n][m] = basis[p][bp]*basis[q][bq]
- self.basisDeriv[n][m][0] = basisDeriv[p][bp][0]*basis[q][bq]
- self.basisDeriv[n][m][1] = basis[p][bp]*basisDeriv[q][bq][0]
- m += 1
- bp = 1; bq = 0
- self.basis[n][m] = basis[p][bp]*basis[q][bq]
- self.basisDeriv[n][m][0] = basisDeriv[p][bp][0]*basis[q][bq]
- self.basisDeriv[n][m][1] = basis[p][bp]*basisDeriv[q][bq][0]
- m += 1
- bp = 1; bq = 1
- self.basis[n][m] = basis[p][bp]*basis[q][bq]
- self.basisDeriv[n][m][0] = basisDeriv[p][bp][0]*basis[q][bq]
- self.basisDeriv[n][m][1] = basis[p][bp]*basisDeriv[q][bq][0]
- m += 1
- bp = 0; bq = 1
- self.basis[n][m] = basis[p][bp]*basis[q][bq]
- self.basisDeriv[n][m][0] = basisDeriv[p][bp][0]*basis[q][bq]
- self.basisDeriv[n][m][1] = basis[p][bp]*basisDeriv[q][bq][0]
- m += 1
-
- # Edges
- # Bottom
- for bp in range(2, numBasisFns):
- bq = 0
+ for (bp,bq) in vertexOrder:
self.basis[n][m] = basis[p][bp]*basis[q][bq]
self.basisDeriv[n][m][0] = basisDeriv[p][bp][0]*basis[q][bq]
self.basisDeriv[n][m][1] = basis[p][bp]*basisDeriv[q][bq][0]
m += 1
-
- # Right
- for bq in range(2, numBasisFns):
- bp = 1
- self.basis[n][m] = basis[p][bp]*basis[q][bq]
- self.basisDeriv[n][m][0] = basisDeriv[p][bp][0]*basis[q][bq]
- self.basisDeriv[n][m][1] = basis[p][bp]*basisDeriv[q][bq][0]
- m += 1
-
- # Top
- for bp in range(numBasisFns-1, 1, -1):
- bq = 1
- self.basis[n][m] = basis[p][bp]*basis[q][bq]
- self.basisDeriv[n][m][0] = basisDeriv[p][bp][0]*basis[q][bq]
- self.basisDeriv[n][m][1] = basis[p][bp]*basisDeriv[q][bq][0]
- m += 1
-
- # Left
- for bq in range(numBasisFns-1, 1, -1):
- bp = 0
- self.basis[n][m] = basis[p][bp]*basis[q][bq]
- self.basisDeriv[n][m][0] = basisDeriv[p][bp][0]*basis[q][bq]
- self.basisDeriv[n][m][1] = basis[p][bp]*basisDeriv[q][bq][0]
- m += 1
-
- # Face
- for bq in range(2, numBasisFns):
- for bp in range(2, numBasisFns):
- self.basis[n][m] = basis[p][bp]*basis[q][bq]
- self.basisDeriv[n][m][0] = basisDeriv[p][bp][0]*basis[q][bq]
- self.basisDeriv[n][m][1] = basis[p][bp]*basisDeriv[q][bq][0]
- m += 1
-
- if not m == numBasisFns**2:
- raise RuntimeError('Invalid 2-D quadrature')
+ if not m == self.numCorners:
+ raise RuntimeError('Invalid 2-D basis tabulation: '+str(m)+ \
+ ' should be '+str(self.numCorners))
n += 1
-
if not n == self.numQuadPts:
- raise RuntimeError('Invalid 2-D quadrature')
+ raise RuntimeError('Invalid 2-D quadrature: '+str(n)+ \
+ ' should be '+str(self.numQuadPts))
elif dim == 3:
- self.vertices = numpy.zeros((self.numCorners, dim))
- n = 0
+ # Set order of vertices and basis functions.
# 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
-
+ vertexOrder = [(0,0,0), (1,0,0), (1,1,0), (0,1,0),
+ (0,0,1), (1,0,1), (1,1,1), (0,1,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
+ p = numpy.arange(2, numBasis, dtype=numpy.int32)
+ q = numpy.zeros(numBasis-2, dtype=numpy.int32)
+ r = numpy.zeros(numBasis-2, dtype=numpy.int32)
+ vertexOrder += zip(p,q,r)
# 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
-
+ p = numpy.ones(numBasis-2, dtype=numpy.int32)
+ q = numpy.arange(2, numBasis, dtype=numpy.int32)
+ r = numpy.zeros(numBasis-2, dtype=numpy.int32)
+ vertexOrder += zip(p,q,r)
# 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
-
+ p = numpy.arange(numBasis-1, 1, step=-1, dtype=numpy.int32)
+ q = numpy.ones(numBasis-2, dtype=numpy.int32)
+ r = numpy.zeros(numBasis-2, dtype=numpy.int32)
+ vertexOrder += zip(p,q,r)
# 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
+ p = numpy.zeros(numBasis-2, dtype=numpy.int32)
+ q = numpy.arange(numBasis-1, 1, step=-1, dtype=numpy.int32)
+ r = numpy.zeros(numBasis-2, dtype=numpy.int32)
+ vertexOrder += zip(p,q,r)
# 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
+ p = numpy.zeros(numBasis-2, dtype=numpy.int32)
+ q = numpy.zeros(numBasis-2, dtype=numpy.int32)
+ r = numpy.arange(2, numBasis, dtype=numpy.int32)
+ vertexOrder += zip(p,q,r)
# 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
-
+ p = numpy.ones(numBasis-2, dtype=numpy.int32)
+ q = numpy.zeros(numBasis-2, dtype=numpy.int32)
+ r = numpy.arange(2, numBasis, dtype=numpy.int32)
+ vertexOrder += zip(p,q,r)
# 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
-
+ p = numpy.ones(numBasis-2, dtype=numpy.int32)
+ q = numpy.ones(numBasis-2, dtype=numpy.int32)
+ r = numpy.arange(2, numBasis, dtype=numpy.int32)
+ vertexOrder += zip(p,q,r)
# 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
+ p = numpy.zeros(numBasis-2, dtype=numpy.int32)
+ q = numpy.ones(numBasis-2, dtype=numpy.int32)
+ r = numpy.arange(2, numBasis, dtype=numpy.int32)
+ vertexOrder += zip(p,q,r)
# 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
+ p = numpy.arange(2, numBasis, dtype=numpy.int32)
+ q = numpy.zeros(numBasis-2, dtype=numpy.int32)
+ r = numpy.ones(numBasis-2, dtype=numpy.int32)
+ vertexOrder += zip(p,q,r)
# 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
-
+ p = numpy.ones(numBasis-2, dtype=numpy.int32)
+ q = numpy.arange(2, numBasis, dtype=numpy.int32)
+ r = numpy.ones(numBasis-2, dtype=numpy.int32)
+ vertexOrder += zip(p,q,r)
# 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
-
+ p = numpy.arange(numBasis-1, 1, step=-1, dtype=numpy.int32)
+ q = numpy.ones(numBasis-2, dtype=numpy.int32)
+ r = numpy.ones(numBasis-2, dtype=numpy.int32)
+ vertexOrder += zip(p,q,r)
# 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
-
+ p = numpy.zeros(numBasis-2, dtype=numpy.int32)
+ q = numpy.arange(numBasis-1, 1, step=-1, dtype=numpy.int32)
+ r = numpy.ones(numBasis-2, dtype=numpy.int32)
+ vertexOrder += zip(p,q,r)
# 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
+ ip = numpy.arange(2, numBasis, dtype=numpy.int32)
+ p = numpy.tile(ip, (1, (numBasis-2)*(numBasis-2)))
+ iq = numpy.arange(2, numBasis, dtype=numpy.int32)
+ q = numpy.tile(iq, ((numBasis-2), numBasis-2)).transpose()
+ ir = numpy.arange(2, numBasis, dtype=numpy.int32)
+ r = numpy.tile(ir, ((numBasis-2)*(numBasis-2), 1)).transpose()
+ vertexOrder += zip(p.ravel(),q.ravel(),r.ravel())
+
+ # Bottom / Top
+ ip = numpy.arange(2, numBasis, dtype=numpy.int32)
+ p = numpy.tile(ip, (1, 2*(numBasis-2)))
+ iq = numpy.arange(2, numBasis, dtype=numpy.int32)
+ q = numpy.tile(iq, (2, numBasis-2)).transpose()
+ ir = numpy.arange(0, 2, dtype=numpy.int32)
+ r = numpy.tile(ir, ((numBasis-2)*(numBasis-2), 1)).transpose()
+ vertexOrder += zip(p.ravel(),q.ravel(),r.ravel())
+ # Left / Right
+ ip = numpy.arange(0, 2, dtype=numpy.int32)
+ p = numpy.tile(ip, ((numBasis-2)*(numBasis-2), 1)).transpose()
+ iq = numpy.arange(2, numBasis, dtype=numpy.int32)
+ q = numpy.tile(iq, (1, 2*(numBasis-2)))
+ ir = numpy.arange(2, numBasis, dtype=numpy.int32)
+ r = numpy.tile(ir, (2, numBasis-2)).transpose()
+ vertexOrder += zip(p.ravel(),q.ravel(),r.ravel())
+ # Front / Back
+ ip = numpy.arange(2, numBasis, dtype=numpy.int32)
+ p = numpy.tile(ip, (1, 2*(numBasis-2)))
+ iq = numpy.arange(0, 2, dtype=numpy.int32)
+ q = numpy.tile(iq, ((numBasis-2)*(numBasis-2), 1)).transpose()
+ ir = numpy.arange(2, numBasis, dtype=numpy.int32)
+ r = numpy.tile(ir, (2, numBasis-2)).transpose()
+ vertexOrder += zip(p.ravel(),q.ravel(),r.ravel())
- # 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
- # 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[1]
- n += 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[r]
- n += 1
- # 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
-
+ self.vertices = numpy.zeros((self.numCorners, dim))
+ n = 0
+ for (p,q,r) in vertexOrder:
+ self.vertices[n][0] = vertices[p]
+ self.vertices[n][1] = vertices[q]
+ self.vertices[n][2] = vertices[r]
+ n += 1
+ print "VERTEX ORDER",vertexOrder
+ print "VERTICES",self.vertices
if not n == self.numCorners:
- raise RuntimeError('Invalid 3-D function tabulation: '+str(n)+ \
- ' should be '+str(self.numCorners))
+ raise RuntimeError('Invalid 3-D vertex ordering: '+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,
- numBasisFns*numBasisFns*numBasisFns))
+ numBasis*numBasis*numBasis))
self.basisDeriv = numpy.zeros((numQuadPts*numQuadPts*numQuadPts,
- numBasisFns*numBasisFns*numBasisFns,
+ numBasis*numBasis*numBasis,
dim))
+
+ # Order of quadrature points doesn't matter
+ # Order of basis functions should match vertices for isoparametric
n = 0
- # Depth
- for s in range(numQuadPts):
- # Bottom
- for r in range(0, numQuadPts-1):
- self.quadPts[n][0] = quadpts[r]
- self.quadPts[n][1] = quadpts[0]
- self.quadPts[n][2] = quadpts[s]
- self.quadWts[n] = quadwts[r]*quadwts[0]*quadwts[s]
- m = 0
- for h in range(numBasisFns):
- # Bottom
- for g in range(0, numBasisFns-1):
- self.basis[n][m] = basis[r][g]*basis[0][0]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[r][g][0]*basis[0][0]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[r][g]*basisDeriv[0][0][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[r][g]*basis[0][0]*basisDeriv[s][h][0]
- m += 1
- # Right
- for f in range(0, numBasisFns-1):
- self.basis[n][m] = basis[r][numBasisFns-1]*basis[0][f]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[r][numBasisFns-1][0]*basis[0][f]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[r][numBasisFns-1]*basisDeriv[0][f][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[r][numBasisFns-1]*basis[0][f]*basisDeriv[s][h][0]
- m += 1
- # Top
- for g in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[r][g]*basis[0][numBasisFns-1]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[r][g][0]*basis[0][numBasisFns-1]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[r][g]*basisDeriv[0][numBasisFns-1][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[r][g]*basis[0][numBasisFns-1]*basisDeriv[s][h][0]
- m += 1
- # Left
- for f in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[r][0]*basis[0][f]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[r][0][0]*basis[0][f]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[r][0]*basisDeriv[0][f][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[r][0]*basis[0][f]*basisDeriv[s][h][0]
- m += 1
- # Interior
- for f in range(1, numBasisFns-1):
- for g in range(1, numBasisFns-1):
- self.basis[n][m] = basis[r][g]*basis[0][f]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[r][g][0]*basis[0][f]*basis[s][h]
- self.basisDeriv[m][m][1] = basis[r][g]*basisDeriv[0][f][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[r][g]*basis[0][f]*basisDeriv[s][h][0]
- m += 1
- if not m == self.numCorners: raise RuntimeError('Invalid 3D function tabulation')
- n += 1
- # Right
- for q in range(0, numQuadPts-1):
- self.quadPts[n][0] = quadpts[numQuadPts-1]
- self.quadPts[n][1] = quadpts[q]
- self.quadPts[n][2] = quadpts[s]
- self.quadWts[n] = quadwts[numQuadPts-1]*quadwts[q]*quadwts[s]
- m = 0
- for h in range(numBasisFns):
- # Bottom
- for g in range(0, numBasisFns-1):
- self.basis[n][m] = basis[numQuadPts-1][g]*basis[q][0]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[numQuadPts-1][g][0]*basis[q][0]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[numQuadPts-1][g]*basisDeriv[q][0][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[numQuadPts-1][g]*basis[q][0]*basisDeriv[s][h][0]
- m += 1
- # Right
- for f in range(0, numBasisFns-1):
- self.basis[n][m] = basis[numQuadPts-1][numBasisFns-1]*basis[q][f]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[numQuadPts-1][numBasisFns-1][0]*basis[q][f]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[numQuadPts-1][numBasisFns-1]*basisDeriv[q][f][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[numQuadPts-1][numBasisFns-1]*basis[q][f]*basisDeriv[s][h][0]
- m += 1
- # Top
- for g in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[numQuadPts-1][g]*basis[q][numBasisFns-1]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[numQuadPts-1][g][0]*basis[q][numBasisFns-1]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[numQuadPts-1][g]*basisDeriv[q][numBasisFns-1][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[numQuadPts-1][g]*basis[q][numBasisFns-1]*basisDeriv[s][h][0]
- m += 1
- # Left
- for f in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[numQuadPts-1][0]*basis[q][f]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[numQuadPts-1][0][0]*basis[q][f]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[numQuadPts-1][0]*basisDeriv[q][f][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[numQuadPts-1][0]*basis[q][f]*basisDeriv[s][h][0]
- m += 1
- # Interior
- for f in range(1, numBasisFns-1):
- for g in range(1, numBasisFns-1):
- self.basis[n][m] = basis[numQuadPts-1][g]*basis[q][f]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[numQuadPts-1][g][0]*basis[q][f]*basis[s][h]
- self.basisDeriv[m][m][1] = basis[numQuadPts-1][g]*basisDeriv[q][f][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[numQuadPts-1][g]*basis[q][f]*basisDeriv[s][h][0]
- m += 1
- if not m == self.numCorners: raise RuntimeError('Invalid 3D function tabulation')
- n += 1
- # Top
- for r in range(numQuadPts-1, 0, -1):
- self.quadPts[n][0] = quadpts[r]
- self.quadPts[n][1] = quadpts[numQuadPts-1]
- self.quadPts[n][2] = quadpts[s]
- self.quadWts[n] = quadwts[r]*quadwts[numQuadPts-1]*quadwts[s]
- m = 0
- for h in range(numBasisFns):
- # Bottom
- for g in range(0, numBasisFns-1):
- self.basis[n][m] = basis[r][g]*basis[numQuadPts-1][0]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[r][g][0]*basis[numQuadPts-1][0]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[r][g]*basisDeriv[numQuadPts-1][0][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[r][g]*basis[numQuadPts-1][0]*basisDeriv[s][h][0]
- m += 1
- # Right
- for f in range(0, numBasisFns-1):
- self.basis[n][m] = basis[r][numBasisFns-1]*basis[numQuadPts-1][f]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[r][numBasisFns-1][0]*basis[numQuadPts-1][f]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[r][numBasisFns-1]*basisDeriv[numQuadPts-1][f][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[r][numBasisFns-1]*basis[numQuadPts-1][f]*basisDeriv[s][h][0]
- m += 1
- # Top
- for g in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[r][g]*basis[numQuadPts-1][numBasisFns-1]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[r][g][0]*basis[numQuadPts-1][numBasisFns-1]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[r][g]*basisDeriv[numQuadPts-1][numBasisFns-1][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[r][g]*basis[numQuadPts-1][numBasisFns-1]*basisDeriv[s][h][0]
- m += 1
- # Left
- for f in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[r][0]*basis[numQuadPts-1][f]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[r][0][0]*basis[numQuadPts-1][f]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[r][0]*basisDeriv[numQuadPts-1][f][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[r][0]*basis[numQuadPts-1][f]*basisDeriv[s][h][0]
- m += 1
- # Interior
- for f in range(1, numBasisFns-1):
- for g in range(1, numBasisFns-1):
- self.basis[n][m] = basis[r][g]*basis[numQuadPts-1][f]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[r][g][0]*basis[numQuadPts-1][f]*basis[s][h]
- self.basisDeriv[m][m][1] = basis[r][g]*basisDeriv[numQuadPts-1][f][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[r][g]*basis[numQuadPts-1][f]*basisDeriv[s][h][0]
- m += 1
- if not m == self.numCorners: raise RuntimeError('Invalid 3D function tabulation')
- n += 1
- # Left
- for q in range(numQuadPts-1, 0, -1):
- self.quadPts[n][0] = quadpts[0]
- self.quadPts[n][1] = quadpts[q]
- self.quadPts[n][2] = quadpts[s]
- self.quadWts[n] = quadwts[0]*quadwts[q]*quadwts[s]
- m = 0
- for h in range(numBasisFns):
- # Bottom
- for g in range(0, numBasisFns-1):
- self.basis[n][m] = basis[0][g]*basis[q][0]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[0][g][0]*basis[q][0]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[0][g]*basisDeriv[q][0][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[0][g]*basis[q][0]*basisDeriv[s][h][0]
- m += 1
- # Right
- for f in range(0, numBasisFns-1):
- self.basis[n][m] = basis[0][numBasisFns-1]*basis[q][f]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[0][numBasisFns-1][0]*basis[q][f]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[0][numBasisFns-1]*basisDeriv[q][f][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[0][numBasisFns-1]*basis[q][f]*basisDeriv[s][h][0]
- m += 1
- # Top
- for g in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[0][g]*basis[q][numBasisFns-1]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[0][g][0]*basis[q][numBasisFns-1]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[0][g]*basisDeriv[q][numBasisFns-1][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[0][g]*basis[q][numBasisFns-1]*basisDeriv[s][h][0]
- m += 1
- # Left
- for f in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[0][0]*basis[q][f]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[0][0][0]*basis[q][f]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[0][0]*basisDeriv[q][f][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[0][0]*basis[q][f]*basisDeriv[s][h][0]
- m += 1
- # Interior
- for f in range(1, numBasisFns-1):
- for g in range(1, numBasisFns-1):
- self.basis[n][m] = basis[0][g]*basis[q][f]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[0][g][0]*basis[q][f]*basis[s][h]
- self.basisDeriv[m][m][1] = basis[0][g]*basisDeriv[q][f][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[0][g]*basis[q][f]*basisDeriv[s][h][0]
- m += 1
- if not m == self.numCorners: raise RuntimeError('Invalid 3D function tabulation')
- n += 1
- # Interior
- for q in range(1, numQuadPts-1):
- for r in range(1, numQuadPts-1):
- self.quadPts[n][0] = quadpts[r]
+ for r in range(0, numQuadPts):
+ for q in range(0, numQuadPts):
+ for p in range(0, numQuadPts):
+ self.quadPts[n][0] = quadpts[p]
self.quadPts[n][1] = quadpts[q]
- self.quadPts[n][2] = quadpts[s]
- self.quadWts[n] = quadwts[r]*quadwts[q]*quadwts[s]
+ self.quadPts[n][2] = quadpts[r]
+ self.quadWts[n] = quadwts[p]*quadwts[q]*quadwts[r]
+
m = 0
- for h in range(numBasisFns):
- # Bottom
- for g in range(0, numBasisFns-1):
- self.basis[n][m] = basis[r][g]*basis[q][0]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[r][g][0]*basis[q][0]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[r][g]*basisDeriv[q][0][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[r][g]*basis[q][0]*basisDeriv[s][h][0]
- m += 1
- # Right
- for f in range(0, numBasisFns-1):
- self.basis[n][m] = basis[r][numBasisFns-1]*basis[q][f]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[r][numBasisFns-1][0]*basis[q][f]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[r][numBasisFns-1]*basisDeriv[q][f][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[r][numBasisFns-1]*basis[q][f]*basisDeriv[s][h][0]
- m += 1
- # Top
- for g in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[r][g]*basis[q][numBasisFns-1]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[r][g][0]*basis[q][numBasisFns-1]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[r][g]*basisDeriv[q][numBasisFns-1][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[r][g]*basis[q][numBasisFns-1]*basisDeriv[s][h][0]
- m += 1
- # Left
- for f in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[r][0]*basis[q][f]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[r][0][0]*basis[q][f]*basis[s][h]
- self.basisDeriv[n][m][1] = basis[r][0]*basisDeriv[q][f][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[r][0]*basis[q][f]*basisDeriv[s][h][0]
- m += 1
- # Interior
- for f in range(1, numBasisFns-1):
- for g in range(1, numBasisFns-1):
- self.basis[n][m] = basis[r][g]*basis[q][f]*basis[s][h]
- self.basisDeriv[n][m][0] = basisDeriv[r][g][0]*basis[q][f]*basis[s][h]
- self.basisDeriv[m][m][1] = basis[r][g]*basisDeriv[q][f][0]*basis[s][h]
- self.basisDeriv[n][m][2] = basis[r][g]*basis[q][f]*basisDeriv[s][h][0]
- m += 1
- if not m == self.numCorners: raise RuntimeError('Invalid 3D function tabulation')
+ for (bp,bq,br) in vertexOrder:
+ self.basis[n][m] = basis[p][bp]*basis[q][bq]*basis[r][br]
+ self.basisDeriv[n][m][0] = basisDeriv[p][bp][0]*basis[q][bq]*basis[r][br]
+ self.basisDeriv[n][m][1] = basis[p][bp]*basisDeriv[q][bq][0]*basis[r][br]
+ self.basisDeriv[n][m][2] = basis[p][bp]*basis[q][bq]*basisDeriv[r][br][0]
+ m += 1
+
+ if not m == self.numCorners:
+ raise RuntimeError('Invalid 3-D basis tabulation: '+str(m)+ \
+ ' should be '+str(self.numCorners))
n += 1
- if not n == self.numQuadPts: raise RuntimeError('Invalid 3D quadrature')
+ if not n == self.numQuadPts:
+ raise RuntimeError('Invalid 3-D quadrature: '+str(n)+ \
+ ' should be '+str(self.numQuadPts))
+
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))
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 16:52:31 UTC (rev 17134)
+++ short/3D/PyLith/branches/v1.5-stable/unittests/pytests/feassemble/TestFIATLagrange.py 2010-08-27 21:02:28 UTC (rev 17135)
@@ -245,7 +245,6 @@
self.N3(q), self.N4(q), self.N5(q),
self.N6(q), self.N7(q), self.N8(q)],
dtype=numpy.float64).reshape( (9,) )
- print basis[iQuad]
deriv = numpy.array([[self.N0p(q), self.N0q(q)],
[self.N1p(q), self.N1q(q)],
[self.N2p(q), self.N2q(q)],
@@ -368,12 +367,12 @@
[-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],
- [-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
@@ -504,6 +503,7 @@
def N7r(self, p):
return (1-p[0])*(1+p[1])/8.0
+
# ----------------------------------------------------------------------
class Hex27(object):
@@ -538,26 +538,58 @@
[+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])
+ quadPts = numpy.array([ [-(3.0/5)**0.5, -(3.0/5)**0.5, -(3.0/5)**0.5],
+ [ 0.0, -(3.0/5)**0.5, -(3.0/5)**0.5],
+ [+(3.0/5)**0.5, -(3.0/5)**0.5, -(3.0/5)**0.5],
+ [-(3.0/5)**0.5, 0.0, -(3.0/5)**0.5],
+ [ 0.0, 0.0, -(3.0/5)**0.5],
+ [+(3.0/5)**0.5, 0.0, -(3.0/5)**0.5],
+ [-(3.0/5)**0.5, +(3.0/5)**0.5, -(3.0/5)**0.5],
+ [ 0.0, +(3.0/5)**0.5, -(3.0/5)**0.5],
+ [+(3.0/5)**0.5, +(3.0/5)**0.5, -(3.0/5)**0.5],
+ [-(3.0/5)**0.5, -(3.0/5)**0.5, 0.0],
+ [ 0.0, -(3.0/5)**0.5, 0.0],
+ [+(3.0/5)**0.5, -(3.0/5)**0.5, 0.0],
+ [-(3.0/5)**0.5, 0.0, 0.0],
+ [ 0.0, 0.0, 0.0],
+ [+(3.0/5)**0.5, 0.0, 0.0],
+ [-(3.0/5)**0.5, +(3.0/5)**0.5, 0.0],
+ [ 0.0, +(3.0/5)**0.5, 0.0],
+ [+(3.0/5)**0.5, +(3.0/5)**0.5, 0.0],
+ [-(3.0/5)**0.5, -(3.0/5)**0.5, +(3.0/5)**0.5],
+ [ 0.0, -(3.0/5)**0.5, +(3.0/5)**0.5],
+ [+(3.0/5)**0.5, -(3.0/5)**0.5, +(3.0/5)**0.5],
+ [-(3.0/5)**0.5, 0.0, +(3.0/5)**0.5],
+ [ 0.0, 0.0, +(3.0/5)**0.5],
+ [+(3.0/5)**0.5, 0.0, +(3.0/5)**0.5],
+ [-(3.0/5)**0.5, +(3.0/5)**0.5, +(3.0/5)**0.5],
+ [ 0.0, +(3.0/5)**0.5, +(3.0/5)**0.5],
+ [+(3.0/5)**0.5, +(3.0/5)**0.5, +(3.0/5)**0.5] ])
+ quadWts = numpy.array( [125.0/729, 200.0/729, 125.0/729,
+ 200.0/729, 320.0/729, 200.0/729,
+ 125.0/729, 200.0/729, 125.0/729,
+ 200.0/729, 320.0/729, 200.0/729,
+ 320.0/729, 512.0/729, 320.0/729,
+ 200.0/729, 320.0/729, 200.0/729,
+ 125.0/729, 200.0/729, 125.0/729,
+ 200.0/729, 320.0/729, 200.0/729,
+ 125.0/729, 200.0/729, 125.0/729])
# 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)
+ basis = numpy.zeros( (27, 27), dtype=numpy.float64)
+ basisDeriv = numpy.zeros( (27, 27, 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,) )
+ basis[iQuad] = numpy.array([self.N0(q), self.N1(q), self.N2(q), self.N3(q), # Corners
+ self.N4(q), self.N5(q), self.N6(q), self.N7(q),
+ self.N8(q), self.N9(q), self.N10(q), self.N11(q), # Edges
+ self.N12(q), self.N13(q), self.N14(q), self.N15(q),
+ self.N16(q), self.N17(q), self.N18(q), self.N19(q),
+ self.N20(q), # Interior
+ self.N21(q), self.N22(q), # Faces
+ self.N23(q), self.N24(q),
+ self.N25(q), self.N26(q)],
+ dtype=numpy.float64).reshape( (27,) )
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)],
@@ -565,8 +597,27 @@
[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))
+ [self.N7p(q), self.N7q(q), self.N7r(q)],
+ [self.N8p(q), self.N8q(q), self.N8r(q)],
+ [self.N9p(q), self.N9q(q), self.N9r(q)],
+ [self.N10p(q), self.N10q(q), self.N10r(q)],
+ [self.N11p(q), self.N11q(q), self.N11r(q)],
+ [self.N12p(q), self.N12q(q), self.N12r(q)],
+ [self.N13p(q), self.N13q(q), self.N13r(q)],
+ [self.N14p(q), self.N14q(q), self.N14r(q)],
+ [self.N15p(q), self.N15q(q), self.N15r(q)],
+ [self.N16p(q), self.N16q(q), self.N16r(q)],
+ [self.N17p(q), self.N17q(q), self.N17r(q)],
+ [self.N18p(q), self.N18q(q), self.N18r(q)],
+ [self.N19p(q), self.N19q(q), self.N19r(q)],
+ [self.N20p(q), self.N20q(q), self.N20r(q)],
+ [self.N21p(q), self.N21q(q), self.N21r(q)],
+ [self.N22p(q), self.N22q(q), self.N22r(q)],
+ [self.N23p(q), self.N23q(q), self.N23r(q)],
+ [self.N24p(q), self.N24q(q), self.N24r(q)],
+ [self.N25p(q), self.N25q(q), self.N25r(q)],
+ [self.N26p(q), self.N26q(q), self.N26r(q)]])
+ basisDeriv[iQuad] = deriv.reshape((27, 3))
iQuad += 1
self.cellDim = 3
@@ -581,101 +632,356 @@
def N0(self, p):
- return (1-p[0])*(1-p[1])*(1-p[2])/8.0
-
+ return (-0.5*p[0])*(1.0-p[0])*(-0.5*p[1])*(1.0-p[1])*(-0.5*p[2])*(1.0-p[2])
+
def N0p(self, p):
- return -(1-p[1])*(1-p[2])/8.0
-
+ return (p[0]-0.5)*(-0.5*p[1])*(1.0-p[1])*(-0.5*p[2])*(1.0-p[2])
+
def N0q(self, p):
- return -(1-p[0])*(1-p[2])/8.0
-
+ return (-0.5*p[0])*(1.0-p[0])*(p[1]-0.5)*(-0.5*p[2])*(1.0-p[2])
+
def N0r(self, p):
- return -(1-p[0])*(1-p[1])/8.0
-
+ return (-0.5*p[0])*(1.0-p[0])*(-0.5*p[1])*(1.0-p[1])*(p[2]-0.5)
+
+
def N1(self, p):
- return (1+p[0])*(1-p[1])*(1-p[2])/8.0
-
+ return (+0.5*p[0])*(1.0+p[0])*(-0.5*p[1])*(1.0-p[1])*(-0.5*p[2])*(1.0-p[2])
+
def N1p(self, p):
- return (1-p[1])*(1-p[2])/8.0
-
+ return (p[0]+0.5)*(-0.5*p[1])*(1.0-p[1])*(-0.5*p[2])*(1.0-p[2])
+
def N1q(self, p):
- return -(1+p[0])*(1-p[2])/8.0
-
+ return (+0.5*p[0])*(1.0+p[0])*(p[1]-0.5)*(-0.5*p[2])*(1.0-p[2])
+
def N1r(self, p):
- return -(1+p[0])*(1-p[1])/8.0
-
+ return (+0.5*p[0])*(1.0+p[0])*(-0.5*p[1])*(1.0-p[1])*(p[2]-0.5)
+
+
def N2(self, p):
- return (1+p[0])*(1+p[1])*(1-p[2])/8.0
-
+ return (+0.5*p[0])*(1.0+p[0])*(+0.5*p[1])*(1.0+p[1])*(-0.5*p[2])*(1.0-p[2])
+
def N2p(self, p):
- return (1+p[1])*(1-p[2])/8.0
-
+ return (p[0]+0.5)*(+0.5*p[1])*(1.0+p[1])*(-0.5*p[2])*(1.0-p[2])
+
def N2q(self, p):
- return (1+p[0])*(1-p[2])/8.0
+ return (+0.5*p[0])*(1.0+p[0])*(p[1]+0.5)*(-0.5*p[2])*(1.0-p[2])
def N2r(self, p):
- return -(1+p[0])*(1+p[1])/8.0
+ return (+0.5*p[0])*(1.0+p[0])*(+0.5*p[1])*(1.0+p[1])*(p[2]-0.5)
+
def N3(self, p):
- return (1-p[0])*(1+p[1])*(1-p[2])/8.0
-
+ return (-0.5*p[0])*(1.0-p[0])*(+0.5*p[1])*(1.0+p[1])*(-0.5*p[2])*(1.0-p[2])
+
def N3p(self, p):
- return -(1+p[1])*(1-p[2])/8.0
-
+ return (p[0]-0.5)*(+0.5*p[1])*(1.0+p[1])*(-0.5*p[2])*(1.0-p[2])
+
def N3q(self, p):
- return (1-p[0])*(1-p[2])/8.0
-
+ return (-0.5*p[0])*(1.0-p[0])*(p[1]+0.5)*(-0.5*p[2])*(1.0-p[2])
+
def N3r(self, p):
- return -(1-p[0])*(1+p[1])/8.0
-
+ return (-0.5*p[0])*(1.0-p[0])*(+0.5*p[1])*(1.0+p[1])*(p[2]-0.5)
+
+
def N4(self, p):
- return (1-p[0])*(1-p[1])*(1+p[2])/8.0
+ return (-0.5*p[0])*(1.0-p[0])*(-0.5*p[1])*(1.0-p[1])*(+0.5*p[2])*(1.0+p[2])
def N4p(self, p):
- return -(1-p[1])*(1+p[2])/8.0
-
+ return (p[0]-0.5)*(-0.5*p[1])*(1.0-p[1])*(+0.5*p[2])*(1.0+p[2])
+
def N4q(self, p):
- return -(1-p[0])*(1+p[2])/8.0
-
+ return (-0.5*p[0])*(1.0-p[0])*(p[1]-0.5)*(+0.5*p[2])*(1.0+p[2])
+
def N4r(self, p):
- return (1-p[0])*(1-p[1])/8.0
-
+ return (-0.5*p[0])*(1.0-p[0])*(-0.5*p[1])*(1.0-p[1])*(p[2]+0.5)
+
+
def N5(self, p):
- return (1+p[0])*(1-p[1])*(1+p[2])/8.0
-
+ return (+0.5*p[0])*(1.0+p[0])*(-0.5*p[1])*(1.0-p[1])*(+0.5*p[2])*(1.0+p[2])
+
def N5p(self, p):
- return (1-p[1])*(1+p[2])/8.0
-
+ return (p[0]+0.5)*(-0.5*p[1])*(1.0-p[1])*(+0.5*p[2])*(1.0+p[2])
+
def N5q(self, p):
- return -(1+p[0])*(1+p[2])/8.0
-
+ return (+0.5*p[0])*(1.0+p[0])*(p[1]-0.5)*(+0.5*p[2])*(1.0+p[2])
+
def N5r(self, p):
- return (1+p[0])*(1-p[1])/8.0
-
+ return (+0.5*p[0])*(1.0+p[0])*(-0.5*p[1])*(1.0-p[1])*(p[2]+0.5)
+
+
def N6(self, p):
- return (1+p[0])*(1+p[1])*(1+p[2])/8.0
-
+ return (+0.5*p[0])*(1.0+p[0])*(+0.5*p[1])*(1.0+p[1])*(+0.5*p[2])*(1.0+p[2])
+
def N6p(self, p):
- return (1+p[1])*(1+p[2])/8.0
-
+ return (p[0]+0.5)*(+0.5*p[1])*(1.0+p[1])*(+0.5*p[2])*(1.0+p[2])
+
def N6q(self, p):
- return (1+p[0])*(1+p[2])/8.0
-
+ return (+0.5*p[0])*(1.0+p[0])*(p[1]+0.5)*(+0.5*p[2])*(1.0+p[2])
+
def N6r(self, p):
- return (1+p[0])*(1+p[1])/8.0
+ return (+0.5*p[0])*(1.0+p[0])*(+0.5*p[1])*(1.0+p[1])*(p[2]+0.5)
+
def N7(self, p):
- return (1-p[0])*(1+p[1])*(1+p[2])/8.0
-
+ return (-0.5*p[0])*(1.0-p[0])*(+0.5*p[1])*(1.0+p[1])*(+0.5*p[2])*(1.0+p[2])
+
def N7p(self, p):
- return -(1+p[1])*(1+p[2])/8.0
-
+ return (p[0]-0.5)*(+0.5*p[1])*(1.0+p[1])*(+0.5*p[2])*(1.0+p[2])
+
def N7q(self, p):
- return (1-p[0])*(1+p[2])/8.0
-
+ return (-0.5*p[0])*(1.0-p[0])*(p[1]+0.5)*(+0.5*p[2])*(1.0+p[2])
+
def N7r(self, p):
- return (1-p[0])*(1+p[1])/8.0
-
+ return (-0.5*p[0])*(1.0-p[0])*(+0.5*p[1])*(1.0+p[1])*(p[2]+0.5)
+
+
+ def N8(self, p):
+ return (1.0-p[0]**2)*(-0.5*p[1])*(1.0-p[1])*(-0.5*p[2])*(1.0-p[2])
+
+ def N8p(self, p):
+ return (-2.0*p[0])*(-0.5*p[1])*(1.0-p[1])*(-0.5*p[2])*(1.0-p[2])
+
+ def N8q(self, p):
+ return (1.0-p[0]**2)*(p[1]-0.5)*(-0.5*p[2])*(1.0-p[2])
+
+ def N8r(self, p):
+ return (1.0-p[0]**2)*(-0.5*p[1])*(1.0-p[1])*(p[2]-0.5)
+
+
+ def N9(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(1.0-p[1]**2)*(-0.5*p[2])*(1.0-p[2])
+
+ def N9p(self, p):
+ return (p[0]+0.5)*(1.0-p[1]**2)*(-0.5*p[2])*(1.0-p[2])
+
+ def N9q(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(-2.0*p[1])*(-0.5*p[2])*(1.0-p[2])
+
+ def N9r(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(1.0-p[1]**2)*(p[2]-0.5)
+
+
+ def N10(self, p):
+ return (1.0-p[0]**2)*(+0.5*p[1])*(1.0+p[1])*(-0.5*p[2])*(1.0-p[2])
+
+ def N10p(self, p):
+ return (-2.0*p[0])*(+0.5*p[1])*(1.0+p[1])*(-0.5*p[2])*(1.0-p[2])
+
+ def N10q(self, p):
+ return (1.0-p[0]**2)*(p[1]+0.5)*(-0.5*p[2])*(1.0-p[2])
+
+ def N10r(self, p):
+ return (1.0-p[0]**2)*(+0.5*p[1])*(1.0+p[1])*(p[2]-0.5)
+
+
+ def N11(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(1.0-p[1]**2)*(-0.5*p[2])*(1.0-p[2])
+
+ def N11p(self, p):
+ return (p[0]-0.5)*(1.0-p[1]**2)*(-0.5*p[2])*(1.0-p[2])
+
+ def N11q(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(-2.0*p[1])*(-0.5*p[2])*(1.0-p[2])
+
+ def N11r(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(1.0-p[1]**2)*(p[2]-0.5)
+
+
+ def N12(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(-0.5*p[1])*(1.0-p[1])*(1.0-p[2]**2)
+
+ def N12p(self, p):
+ return (p[0]-0.5)*(-0.5*p[1])*(1.0-p[1])*(1.0-p[2]**2)
+
+ def N12q(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(p[1]-0.5)*(1.0-p[2]**2)
+
+ def N12r(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(-0.5*p[1])*(1.0-p[1])*(-2.0*p[2])
+
+
+ def N13(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(-0.5*p[1])*(1.0-p[1])*(1.0-p[2]**2)
+
+ def N13p(self, p):
+ return (p[0]+0.5)*(-0.5*p[1])*(1.0-p[1])*(1.0-p[2]**2)
+
+ def N13q(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(p[1]-0.5)*(1.0-p[2]**2)
+
+ def N13r(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(-0.5*p[1])*(1.0-p[1])*(-2.0*p[2])
+
+
+ def N14(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(+0.5*p[1])*(1.0+p[1])*(1.0-p[2]**2)
+
+ def N14p(self, p):
+ return (p[0]+0.5)*(+0.5*p[1])*(1.0+p[1])*(1.0-p[2]**2)
+
+ def N14q(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(p[1]+0.5)*(1.0-p[2]**2)
+
+ def N14r(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(+0.5*p[1])*(1.0+p[1])*(-2.0*p[2])
+
+
+ def N15(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(+0.5*p[1])*(1.0+p[1])*(1.0-p[2]**2)
+
+ def N15p(self, p):
+ return (p[0]-0.5)*(+0.5*p[1])*(1.0+p[1])*(1.0-p[2]**2)
+
+ def N15q(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(p[1]+0.5)*(1.0-p[2]**2)
+
+ def N15r(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(+0.5*p[1])*(1.0+p[1])*(-2.0*p[2])
+
+
+ def N16(self, p):
+ return (1.0-p[0]**2)*(-0.5*p[1])*(1.0-p[1])*(+0.5*p[2])*(1.0+p[2])
+
+ def N16p(self, p):
+ return (-2.0*p[0])*(-0.5*p[1])*(1.0-p[1])*(+0.5*p[2])*(1.0+p[2])
+
+ def N16q(self, p):
+ return (1.0-p[0]**2)*(p[1]-0.5)*(+0.5*p[2])*(1.0+p[2])
+
+ def N16r(self, p):
+ return (1.0-p[0]**2)*(-0.5*p[1])*(1.0-p[1])*(p[2]+0.5)
+
+
+ def N17(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(1.0-p[1]**2)*(+0.5*p[2])*(1.0+p[2])
+
+ def N17p(self, p):
+ return (p[0]+0.5)*(1.0-p[1]**2)*(+0.5*p[2])*(1.0+p[2])
+
+ def N17q(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(-2.0*p[1])*(+0.5*p[2])*(1.0+p[2])
+
+ def N17r(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(1.0-p[1]**2)*(p[2]+0.5)
+
+
+ def N18(self, p):
+ return (1.0-p[0]**2)*(+0.5*p[1])*(1.0+p[1])*(+0.5*p[2])*(1.0+p[2])
+
+ def N18p(self, p):
+ return (-2.0*p[0])*(+0.5*p[1])*(1.0+p[1])*(+0.5*p[2])*(1.0+p[2])
+
+ def N18q(self, p):
+ return (1.0-p[0]**2)*(p[1]+0.5)*(+0.5*p[2])*(1.0+p[2])
+
+ def N18r(self, p):
+ return (1.0-p[0]**2)*(+0.5*p[1])*(1.0+p[1])*(p[2]+0.5)
+
+
+ def N19(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(1.0-p[1]**2)*(+0.5*p[2])*(1.0+p[2])
+
+ def N19p(self, p):
+ return (p[0]-0.5)*(1.0-p[1]**2)*(+0.5*p[2])*(1.0+p[2])
+
+ def N19q(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(-2.0*p[1])*(+0.5*p[2])*(1.0+p[2])
+
+ def N19r(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(1.0-p[1]**2)*(p[2]+0.5)
+
+
+ def N20(self, p):
+ return (1.0-p[0]**2)*(1.0-p[1]**2)*(1.0-p[2]**2)
+
+ def N20p(self, p):
+ return (-2.0*p[0])*(1.0-p[1]**2)*(1.0-p[2]**2)
+
+ def N20q(self, p):
+ return (1.0-p[0]**2)*(-2.0*p[1])*(1.0-p[2]**2)
+
+ def N20r(self, p):
+ return (1.0-p[0]**2)*(1.0-p[1]**2)*(-2.0*p[2])
+
+
+ def N21(self, p):
+ return (1.0-p[0]**2)*(1.0-p[1]**2)*(-0.5*p[2])*(1.0-p[2])
+
+ def N21p(self, p):
+ return (-2.0*p[0])*(1.0-p[1]**2)*(-0.5*p[2])*(1.0-p[2])
+
+ def N21q(self, p):
+ return (1.0-p[0]**2)*(-2.0*p[1])*(-0.5*p[2])*(1.0-p[2])
+
+ def N21r(self, p):
+ return (1.0-p[0]**2)*(1.0-p[1]**2)*(p[2]-0.5)
+
+
+ def N22(self, p):
+ return (1.0-p[0]**2)*(1.0-p[1]**2)*(+0.5*p[2])*(1.0+p[2])
+
+ def N22p(self, p):
+ return (-2.0*p[0])*(1.0-p[1]**2)*(+0.5*p[2])*(1.0+p[2])
+
+ def N22q(self, p):
+ return (1.0-p[0]**2)*(-2.0*p[1])*(+0.5*p[2])*(1.0+p[2])
+
+ def N22r(self, p):
+ return (1.0-p[0]**2)*(1.0-p[1]**2)*(p[2]+0.5)
+
+
+ def N23(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(1.0-p[1]**2)*(1.0-p[2]**2)
+
+ def N23p(self, p):
+ return (p[0]-0.5)*(1.0-p[1]**2)*(1.0-p[2]**2)
+
+ def N23q(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(-2.0*p[1])*(1.0-p[2]**2)
+
+ def N23r(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(1.0-p[1]**2)*(-2.0*p[2])
+
+
+ def N24(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(1.0-p[1]**2)*(1.0-p[2]**2)
+
+ def N24p(self, p):
+ return (p[0]+0.5)*(1.0-p[1]**2)*(1.0-p[2]**2)
+
+ def N24q(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(-2.0*p[1])*(1.0-p[2]**2)
+
+ def N24r(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(1.0-p[1]**2)*(-2.0*p[2])
+
+
+ def N25(self, p):
+ return (1.0-p[0]**2)*(-0.5*p[1])*(1.0-p[1])*(1.0-p[2]**2)
+
+ def N25p(self, p):
+ return (-2.0*p[0])*(-0.5*p[1])*(1.0-p[1])*(1.0-p[2]**2)
+
+ def N25q(self, p):
+ return (1.0-p[0]**2)*(p[1]-0.5)*(1.0-p[2]**2)
+
+ def N25r(self, p):
+ return (1.0-p[0]**2)*(-0.5*p[1])*(1.0-p[1])*(-2.0*p[2])
+
+
+ def N26(self, p):
+ return (1.0-p[0]**2)*(+0.5*p[1])*(1.0+p[1])*(1.0-p[2]**2)
+
+ def N26p(self, p):
+ return (-2.0*p[0])*(+0.5*p[1])*(1.0+p[1])*(1.0-p[2]**2)
+
+ def N26q(self, p):
+ return (1.0-p[0]**2)*(p[1]+0.5)*(1.0-p[2]**2)
+
+ def N26r(self, p):
+ return (1.0-p[0]**2)*(+0.5*p[1])*(1.0+p[1])*(-2.0*p[2])
+
+
# ----------------------------------------------------------------------
class TestFIATLagrange(unittest.TestCase):
"""
@@ -772,6 +1078,24 @@
return
+ def test_initialize_hex27(self):
+ """
+ Test initialize() with hex27 cell.
+ """
+ cell = FIATLagrange()
+ cell.inventory.dimension = 3
+ cell.inventory.degree = 2
+ cell.inventory.order = 5
+ cell._configure()
+ cell.initialize(spaceDim=3)
+
+ cellE = Hex27()
+ self._checkVals(cellE, cell)
+ from pylith.feassemble.CellGeometry import GeometryHex3D
+ self.failUnless(isinstance(cell.geometry, GeometryHex3D))
+ return
+
+
def test_factory(self):
"""
Test factory method.
More information about the CIG-COMMITS
mailing list