[cig-commits] r17126 - in short/3D/PyLith/branches/v1.5-stable: libsrc/faults pylith/feassemble pylith/perf unittests/pytests/feassemble
brad at geodynamics.org
brad at geodynamics.org
Thu Aug 26 16:31:22 PDT 2010
Author: brad
Date: 2010-08-26 16:31:22 -0700 (Thu, 26 Aug 2010)
New Revision: 17126
Modified:
short/3D/PyLith/branches/v1.5-stable/libsrc/faults/CohesiveTopology.cc
short/3D/PyLith/branches/v1.5-stable/pylith/feassemble/FIATLagrange.py
short/3D/PyLith/branches/v1.5-stable/pylith/perf/Mesh.py
short/3D/PyLith/branches/v1.5-stable/unittests/pytests/feassemble/TestFIATLagrange.py
Log:
Added unit test to check FIATLagrange for 2-D cells with degree == 2 (quadratic quadrilaterals). Fixed ordering of basis stuff in FIATLagrange for 2-D. Faults will not work without additional Sieve code to handle additional vertices on faces.
Modified: short/3D/PyLith/branches/v1.5-stable/libsrc/faults/CohesiveTopology.cc
===================================================================
--- short/3D/PyLith/branches/v1.5-stable/libsrc/faults/CohesiveTopology.cc 2010-08-26 21:17:31 UTC (rev 17125)
+++ short/3D/PyLith/branches/v1.5-stable/libsrc/faults/CohesiveTopology.cc 2010-08-26 23:31:22 UTC (rev 17126)
@@ -545,7 +545,7 @@
<< cellNeighbors.size() << std::endl;
throw ALE::Exception("Invalid number of neighbors");
} // if
- } else if (numCorners == 4) {
+ } else if (numCorners == 4 || numCorners == 9) {
if (cellNeighbors.size() > 4) {
std::cout << "Cell " << *c_iter
<< " has an invalid number of neighbors "
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-26 21:17:31 UTC (rev 17125)
+++ short/3D/PyLith/branches/v1.5-stable/pylith/feassemble/FIATLagrange.py 2010-08-26 23:31:22 UTC (rev 17126)
@@ -122,33 +122,53 @@
if dim == 2:
self.vertices = numpy.zeros((self.numCorners, dim))
n = 0
- # Bottom
- for r in range(0, numBasisFns-1):
- self.vertices[n][0] = vertices[r]
+ # 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
+
+ # Edges
+ # Bottom
+ for p in range(2, numBasisFns):
+ self.vertices[n][0] = vertices[p]
self.vertices[n][1] = vertices[0]
n += 1
- # Right
- for q in range(0, numBasisFns-1):
- self.vertices[n][0] = vertices[numBasisFns-1]
+ # Right
+ for q in range(2, numBasisFns):
+ self.vertices[n][0] = vertices[1]
self.vertices[n][1] = vertices[q]
n += 1
- # Top
- for r in range(numBasisFns-1, 0, -1):
- self.vertices[n][0] = vertices[r]
- self.vertices[n][1] = vertices[numBasisFns-1]
+
+ # Top
+ for p in range(numBasisFns-1, 1, -1):
+ self.vertices[n][0] = vertices[p]
+ self.vertices[n][1] = vertices[1]
n += 1
- # Left
- for q in range(numBasisFns-1, 0, -1):
+
+ # Left
+ for r in range(numBasisFns-1, 1, -1):
self.vertices[n][0] = vertices[0]
- self.vertices[n][1] = vertices[q]
+ self.vertices[n][1] = 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]
+
+ # 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 2D function tabulation')
+
+ if not n == self.numCorners:
+ raise RuntimeError('Invalid 2D function tabulation, n is %d' % n)
self.quadPts = numpy.zeros((numQuadPts*numQuadPts, dim))
self.quadWts = numpy.zeros((numQuadPts*numQuadPts,))
@@ -156,203 +176,82 @@
numBasisFns*numBasisFns))
self.basisDeriv = numpy.zeros((numQuadPts*numQuadPts,
numBasisFns*numBasisFns, dim))
+
+ # Order of basis functions and quadrature points don't matter
n = 0
- # Bottom
- for r in range(0, numQuadPts-1):
- self.quadPts[n][0] = quadpts[r]
- self.quadPts[n][1] = quadpts[0]
- self.quadWts[n] = quadwts[r]*quadwts[0]
- m = 0
- # Bottom
- for g in range(0, numBasisFns-1):
- self.basis[n][m] = basis[r][g]*basis[0][0]
- self.basisDeriv[n][m][0] = basisDeriv[r][g][0]*basis[0][0]
- self.basisDeriv[n][m][1] = basis[r][g]*basisDeriv[0][0][0]
+ 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.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
- # Right
- for f in range(0, numBasisFns-1):
- self.basis[n][m] = basis[r][numBasisFns-1]*basis[0][f]
- self.basisDeriv[n][m][0] = basisDeriv[r][numBasisFns-1][0]*basis[0][f]
- self.basisDeriv[n][m][1] = basis[r][numBasisFns-1]*basisDeriv[0][f][0]
+ 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
- # Top
- for g in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[r][g]*basis[0][numBasisFns-1]
- self.basisDeriv[n][m][0] = basisDeriv[r][g][0]*basis[0][numBasisFns-1]
- self.basisDeriv[n][m][1] = basis[r][g]*basisDeriv[0][numBasisFns-1][0]
+ 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
- # Left
- for f in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[r][0]*basis[0][f]
- self.basisDeriv[n][m][0] = basisDeriv[r][0][0]*basis[0][f]
- self.basisDeriv[n][m][1] = basis[r][0]*basisDeriv[0][f][0]
+ 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
- # 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]
- self.basisDeriv[0][r][f][g][0] = basisDeriv[r][g][0]*basis[0][f]
- self.basisDeriv[0][r][f][g][1] = basis[r][g]*basisDeriv[0][f][0]
+
+ # Edges
+ # Bottom
+ for bp in range(2, numBasisFns):
+ 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
- if not m == self.numCorners: raise RuntimeError('Invalid 2D 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.quadWts[n] = quadwts[numQuadPts-1]*quadwts[q]
- m = 0
- # Bottom
- for g in range(0, numBasisFns-1):
- self.basis[n][m] = basis[numQuadPts-1][g]*basis[q][0]
- self.basisDeriv[n][m][0] = basisDeriv[numQuadPts-1][g][0]*basis[q][0]
- self.basisDeriv[n][m][1] = basis[numQuadPts-1][g]*basisDeriv[q][0][0]
- m += 1
- # Right
- for f in range(0, numBasisFns-1):
- self.basis[n][m] = basis[numQuadPts-1][numBasisFns-1]*basis[q][f]
- self.basisDeriv[n][m][0] = basisDeriv[numQuadPts-1][numBasisFns-1][0]*basis[q][f]
- self.basisDeriv[n][m][1] = basis[numQuadPts-1][numBasisFns-1]*basisDeriv[q][f][0]
- m += 1
- # Top
- for g in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[numQuadPts-1][g]*basis[q][numBasisFns-1]
- self.basisDeriv[n][m][0] = basisDeriv[numQuadPts-1][g][0]*basis[q][numBasisFns-1]
- self.basisDeriv[n][m][1] = basis[numQuadPts-1][g]*basisDeriv[q][numBasisFns-1][0]
- m += 1
- # Left
- for f in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[numQuadPts-1][0]*basis[q][f]
- self.basisDeriv[n][m][0] = basisDeriv[numQuadPts-1][0][0]*basis[q][f]
- self.basisDeriv[n][m][1] = basis[numQuadPts-1][0]*basisDeriv[q][f][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[0][f]
- self.basisDeriv[q][numQuadPts-1][f][g][0] = basisDeriv[numQuadPts-1][g][0]*basis[q][f]
- self.basisDeriv[q][numQuadPts-1][f][g][1] = basis[numQuadPts-1][g]*basisDeriv[q][f][0]
+
+ # 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
- if not m == self.numCorners: raise RuntimeError('Invalid 2D 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.quadWts[n] = quadwts[r]*quadwts[numQuadPts-1]
- m = 0
- # Bottom
- for g in range(0, numBasisFns-1):
- self.basis[n][m] = basis[r][g]*basis[numQuadPts-1][0]
- self.basisDeriv[n][m][0] = basisDeriv[r][g][0]*basis[numQuadPts-1][0]
- self.basisDeriv[n][m][1] = basis[r][g]*basisDeriv[numQuadPts-1][0][0]
- m += 1
- # Right
- for f in range(0, numBasisFns-1):
- self.basis[n][m] = basis[r][numBasisFns-1]*basis[numQuadPts-1][f]
- self.basisDeriv[n][m][0] = basisDeriv[r][numBasisFns-1][0]*basis[numQuadPts-1][f]
- self.basisDeriv[n][m][1] = basis[r][numBasisFns-1]*basisDeriv[numQuadPts-1][f][0]
- m += 1
- # Top
- for g in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[r][g]*basis[numQuadPts-1][numBasisFns-1]
- self.basisDeriv[n][m][0] = basisDeriv[r][g][0]*basis[numQuadPts-1][numBasisFns-1]
- self.basisDeriv[n][m][1] = basis[r][g]*basisDeriv[numQuadPts-1][numBasisFns-1][0]
- m += 1
- # Left
- for f in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[r][0]*basis[numQuadPts-1][f]
- self.basisDeriv[n][m][0] = basisDeriv[r][0][0]*basis[numQuadPts-1][f]
- self.basisDeriv[n][m][1] = basis[r][0]*basisDeriv[numQuadPts-1][f][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]
- self.basisDeriv[numQuadPts-1][r][f][g][0] = basisDeriv[r][g][0]*basis[numQuadPts-1][f]
- self.basisDeriv[numQuadPts-1][r][f][g][1] = basis[r][g]*basisDeriv[numQuadPts-1][f][0]
+
+ # 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
- if not m == self.numCorners: raise RuntimeError('Invalid 2D 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.quadWts[n] = quadwts[0]*quadwts[q]
- m = 0
- # Bottom
- for g in range(0, numBasisFns-1):
- self.basis[n][m] = basis[0][g]*basis[q][0]
- self.basisDeriv[n][m][0] = basisDeriv[0][g][0]*basis[q][0]
- self.basisDeriv[n][m][1] = basis[0][g]*basisDeriv[q][0][0]
- m += 1
- # Right
- for f in range(0, numBasisFns-1):
- self.basis[n][m] = basis[0][numBasisFns-1]*basis[q][f]
- self.basisDeriv[n][m][0] = basisDeriv[0][numBasisFns-1][0]*basis[q][f]
- self.basisDeriv[n][m][1] = basis[0][numBasisFns-1]*basisDeriv[q][f][0]
- m += 1
- # Top
- for g in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[0][g]*basis[q][numBasisFns-1]
- self.basisDeriv[n][m][0] = basisDeriv[0][g][0]*basis[q][numBasisFns-1]
- self.basisDeriv[n][m][1] = basis[0][g]*basisDeriv[q][numBasisFns-1][0]
- m += 1
- # Left
- for f in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[0][0]*basis[q][f]
- self.basisDeriv[n][m][0] = basisDeriv[0][0][0]*basis[q][f]
- self.basisDeriv[n][m][1] = basis[0][0]*basisDeriv[q][f][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[0][f]
- self.basisDeriv[q][0][f][g][0] = basisDeriv[0][g][0]*basis[q][f]
- self.basisDeriv[q][0][f][g][1] = basis[0][g]*basisDeriv[q][f][0]
+
+ # 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
- if not m == self.numCorners: raise RuntimeError('Invalid 2D 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]
- self.quadPts[n][1] = quadpts[q]
- self.quadWts[n] = quadwts[r]*quadwts[q]
- m = 0
- # Bottom
- for g in range(0, numBasisFns-1):
- self.basis[n][m] = basis[r][g]*basis[q][0]
- self.basisDeriv[n][m][0] = basisDeriv[r][g][0]*basis[q][0]
- self.basisDeriv[n][m][1] = basis[r][g]*basisDeriv[q][0][0]
- m += 1
- # Right
- for f in range(0, numBasisFns-1):
- self.basis[n][m] = basis[r][numBasisFns-1]*basis[q][f]
- self.basisDeriv[n][m][0] = basisDeriv[r][numBasisFns-1][0]*basis[q][f]
- self.basisDeriv[n][m][1] = basis[r][numBasisFns-1]*basisDeriv[q][f][0]
- m += 1
- # Top
- for g in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[r][g]*basis[q][numBasisFns-1]
- self.basisDeriv[n][m][0] = basisDeriv[r][g][0]*basis[q][numBasisFns-1]
- self.basisDeriv[n][m][1] = basis[r][g]*basisDeriv[q][numBasisFns-1][0]
- m += 1
- # Left
- for f in range(numBasisFns-1, 0, -1):
- self.basis[n][m] = basis[r][0]*basis[q][f]
- self.basisDeriv[n][m][0] = basisDeriv[r][0][0]*basis[q][f]
- self.basisDeriv[n][m][1] = basis[r][0]*basisDeriv[q][f][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]
- 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]
+
+ # 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 == self.numCorners: raise RuntimeError('Invalid 2D function tabulation')
+
+ if not m == numBasisFns**2: raise RuntimeError('Invalid 2D quadrature')
n += 1
+
if not n == self.numQuadPts: raise RuntimeError('Invalid 2D quadrature')
elif dim == 3:
self.vertices = numpy.zeros((self.numCorners, dim))
Modified: short/3D/PyLith/branches/v1.5-stable/pylith/perf/Mesh.py
===================================================================
--- short/3D/PyLith/branches/v1.5-stable/pylith/perf/Mesh.py 2010-08-26 21:17:31 UTC (rev 17125)
+++ short/3D/PyLith/branches/v1.5-stable/pylith/perf/Mesh.py 2010-08-26 23:31:22 UTC (rev 17126)
@@ -27,6 +27,7 @@
(1,3): 'line3',
(2,3): 'tri3',
(2,4): 'quad4',
+ (2,9): 'quad9',
(3,4): 'tet4',
(3,8): 'hex8'
}
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-26 21:17:31 UTC (rev 17125)
+++ short/3D/PyLith/branches/v1.5-stable/unittests/pytests/feassemble/TestFIATLagrange.py 2010-08-26 23:31:22 UTC (rev 17126)
@@ -72,6 +72,7 @@
def N1p(self, p):
return 0.5
+
# ----------------------------------------------------------------------
class Line3(object):
@@ -124,6 +125,7 @@
def N2p(self, p):
return -2.0*p
+
# ----------------------------------------------------------------------
class Quad4(object):
@@ -137,8 +139,8 @@
[-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] ])
quadWts = numpy.array( [1.0, 1.0, 1.0, 1.0])
# Compute basis fns and derivatives at quadrature points
@@ -203,7 +205,153 @@
def N3q(self, p):
return (1-p[0])/4.0
+
# ----------------------------------------------------------------------
+class Quad9(object):
+
+ def __init__(self):
+ """
+ Setup quad9 cell.
+ """
+ vertices = numpy.array([[-1.0, -1.0], # corners
+ [+1.0, -1.0],
+ [+1.0, +1.0],
+ [-1.0, +1.0],
+ [ 0.0, -1.0], # edges
+ [+1.0, 0.0],
+ [ 0.0, +1.0],
+ [-1.0, 0.0],
+ [ 0.0, 0.0], # face
+ ])
+ quadPts = numpy.array([ [-(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],
+ [ 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] ])
+ quadWts = numpy.array( [25.0/81, 40.0/81, 25.0/81,
+ 40.0/81, 64.0/81, 40.0/81,
+ 25.0/81, 40.0/81, 25.0/81])
+
+ # Compute basis fns and derivatives at quadrature points
+ basis = numpy.zeros( (9, 9), dtype=numpy.float64)
+ basisDeriv = numpy.zeros( (9, 9, 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), 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)],
+ [self.N3p(q), self.N3q(q)],
+ [self.N4p(q), self.N4q(q)],
+ [self.N5p(q), self.N5q(q)],
+ [self.N6p(q), self.N6q(q)],
+ [self.N7p(q), self.N7q(q)],
+ [self.N8p(q), self.N8q(q)]])
+ basisDeriv[iQuad] = deriv.reshape((9, 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])*(1.0-p[0])*(-0.5*p[1])*(1.0-p[1])
+
+ def N0p(self, p):
+ return (p[0]-0.5)*(-0.5*p[1])*(1.0-p[1])
+
+ def N0q(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(p[1]-0.5)
+
+ def N1(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(-0.5*p[1])*(1.0-p[1])
+
+ def N1p(self, p):
+ return (p[0]+0.5)*(-0.5*p[1])*(1.0-p[1])
+
+ def N1q(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(p[1]-0.5)
+
+ def N2(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(+0.5*p[1])*(1.0+p[1])
+
+ def N2p(self, p):
+ return (p[0]+0.5)*(+0.5*p[1])*(1.0+p[1])
+
+ def N2q(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(p[1]+0.5)
+
+ def N3(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(+0.5*p[1])*(1.0+p[1])
+
+ def N3p(self, p):
+ return (p[0]-0.5)*(+0.5*p[1])*(1.0+p[1])
+
+ def N3q(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(p[1]+0.5)
+
+ def N4(self, p):
+ return (1.0-p[0]**2)*(-0.5*p[1])*(1.0-p[1])
+
+ def N4p(self, p):
+ return (-2.0*p[0])*(-0.5*p[1])*(1.0-p[1])
+
+ def N4q(self, p):
+ return (1.0-p[0]**2)*(p[1]-0.5)
+
+ def N5(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(1.0-p[1]**2)
+
+ def N5p(self, p):
+ return (p[0]+0.5)*(1.0-p[1]**2)
+
+ def N5q(self, p):
+ return (+0.5*p[0])*(1.0+p[0])*(-2.0*p[1])
+
+ def N6(self, p):
+ return (1.0-p[0]**2)*(+0.5*p[1])*(1.0+p[1])
+
+ def N6p(self, p):
+ return (-2.0*p[0])*(+0.5*p[1])*(1.0+p[1])
+
+ def N6q(self, p):
+ return (1.0-p[0]**2)*(p[1]+0.5)
+
+ def N7(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(1.0-p[1]**2)
+
+ def N7p(self, p):
+ return (p[0]-0.5)*(1.0-p[1]**2)
+
+ def N7q(self, p):
+ return (-0.5*p[0])*(1.0-p[0])*(-2.0*p[1])
+
+ def N8(self, p):
+ return (1.0-p[0]**2)*(1.0-p[1]**2)
+
+ def N8p(self, p):
+ return (-2.0*p[0])*(1.0-p[1]**2)
+
+ def N8q(self, p):
+ return (1.0-p[0]**2)*(-2.0*p[1])
+
+
+# ----------------------------------------------------------------------
class Hex8(object):
def __init__(self):
@@ -366,11 +514,11 @@
"""
Test initialize() with line2 cell.
"""
- from pylith.feassemble.FIATSimplex import FIATSimplex
cell = FIATLagrange()
- cell.cellDim = 1
- cell.degree = 1
- cell.order = 1
+ cell.inventory.dimension = 1
+ cell.inventory.degree = 1
+ cell.inventory.order = 1
+ cell._configure()
cell.initialize(spaceDim=1)
cellE = Line2()
@@ -384,11 +532,11 @@
"""
Test initialize() with line3 cell.
"""
- from pylith.feassemble.FIATSimplex import FIATSimplex
cell = FIATLagrange()
- cell.cellDim = 1
- cell.degree = 2
- cell.order = 2
+ cell.inventory.dimension = 1
+ cell.inventory.degree = 2
+ cell.inventory.order = 2
+ cell._configure()
cell.initialize(spaceDim=2)
cellE = Line3()
@@ -402,11 +550,11 @@
"""
Test initialize() with quad4 cell.
"""
- from pylith.feassemble.FIATSimplex import FIATSimplex
cell = FIATLagrange()
- cell.cellDim = 2
- cell.degree = 1
- cell.order = 2
+ cell.inventory.dimension = 2
+ cell.inventory.degree = 1
+ cell.inventory.order = 2
+ cell._configure()
cell.initialize(spaceDim=2)
cellE = Quad4()
@@ -416,15 +564,33 @@
return
+ def test_initialize_quad9(self):
+ """
+ Test initialize() with quad9 cell.
+ """
+ cell = FIATLagrange()
+ cell.inventory.dimension = 2
+ cell.inventory.degree = 2
+ cell.inventory.order = 5
+ cell._configure()
+ cell.initialize(spaceDim=2)
+
+ cellE = Quad9()
+ self._checkVals(cellE, cell)
+ from pylith.feassemble.CellGeometry import GeometryQuad2D
+ self.failUnless(isinstance(cell.geometry, GeometryQuad2D))
+ return
+
+
def test_initialize_hex8(self):
"""
Test initialize() with hex8 cell.
"""
- from pylith.feassemble.FIATSimplex import FIATSimplex
cell = FIATLagrange()
- cell.cellDim = 3
- cell.degree = 1
- cell.order = 2
+ cell.inventory.dimension = 3
+ cell.inventory.degree = 1
+ cell.inventory.order = 2
+ cell._configure()
cell.initialize(spaceDim=3)
cellE = Hex8()
@@ -445,7 +611,7 @@
def _checkVals(self, cellE, cell):
"""
- Check known values against those generated by FIATSimplex.
+ Check known values against those generated by FIATLagrange.
"""
# Check basic attributes
More information about the CIG-COMMITS
mailing list