[cig-commits] r7048 - in short/3D/PyLith/trunk: pylith/feassemble
unittests/pytests/feassemble
knepley at geodynamics.org
knepley at geodynamics.org
Sat Jun 2 09:35:16 PDT 2007
Author: knepley
Date: 2007-06-02 09:35:15 -0700 (Sat, 02 Jun 2007)
New Revision: 7048
Modified:
short/3D/PyLith/trunk/pylith/feassemble/FIATLagrange.py
short/3D/PyLith/trunk/unittests/pytests/feassemble/TestFIATLagrange.py
Log:
Test now passes for quadrilateral elements
Modified: short/3D/PyLith/trunk/pylith/feassemble/FIATLagrange.py
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/FIATLagrange.py 2007-06-02 04:30:16 UTC (rev 7047)
+++ short/3D/PyLith/trunk/pylith/feassemble/FIATLagrange.py 2007-06-02 16:35:15 UTC (rev 7048)
@@ -83,34 +83,75 @@
1-D Lagrange elements.
"""
quadrature = self._setupQuadrature()
- basisFns = self._setupBasisFns()
+ element = self._setupElement()
+ dim = self.cellDim
# Evaluate basis functions at quadrature points
- quadpts = quadrature.get_points()
- basis = numpy.array(basisFns.tabulate(quadpts)).transpose()
- self.basis = numpy.reshape(basis.flatten(), basis.shape)
+ quadpts = numpy.array(quadrature.get_points())
+ 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())
# Evaluate derivatives of basis functions at quadrature points
- import FIAT.shapes
- dim = FIAT.shapes.dimension(basisFns.base.shape)
- basisDeriv = numpy.array([basisFns.deriv_all(d).tabulate(quadpts) \
- for d in range(dim)]).transpose()
- self.basisDeriv = numpy.reshape(basisDeriv.flatten(), basisDeriv.shape)
+ basisDeriv = numpy.array([element.function_space().deriv_all(d).tabulate(quadrature.get_points()) \
+ for d in range(1)]).transpose()
- self.quadPts = numpy.array(quadrature.get_points())
- self.quadWts = numpy.array(quadrature.get_weights())
+ self.numQuadPts = numQuadPts**dim
+ self.numCorners = numBasisFns**dim
- self.numCorners = len(basisFns)
- self.numQuadPts = len(quadrature.get_weights())
-
+ if dim == 1:
+ self.quadPts = quadpts
+ self.quadWts = quadwts
+ self.basis = basis
+ self.basisDeriv = basisDeriv
+ else:
+ if dim == 2:
+ self.quadPts = numpy.zeros((numQuadPts, numQuadPts, dim))
+ self.quadWts = numpy.zeros((numQuadPts, numQuadPts))
+ self.basis = numpy.zeros((numQuadPts, numQuadPts, numBasisFns, numBasisFns))
+ self.basisDeriv = numpy.zeros((numQuadPts, numQuadPts, numBasisFns, numBasisFns, dim))
+ for q in range(numQuadPts):
+ for r in range(numQuadPts):
+ self.quadPts[q][r][0] = quadpts[q]
+ self.quadPts[q][r][1] = quadpts[r]
+ self.quadWts[q][r] = quadwts[q]*quadwts[r]
+ for f in range(numBasisFns):
+ for g in range(numBasisFns):
+ self.basis[q][r][f][g] = basis[q][f]*basis[r][g]
+ self.basisDeriv[q][r][f][g][0] = basisDeriv[q][f][0]*basis[r][g]
+ self.basisDeriv[q][r][f][g][1] = basis[q][f]*basisDeriv[r][g][0]
+ elif dim == 3:
+ self.quadPts = numpy.array((numQuadPts, numQuadPts, numQuadPts, dim))
+ self.quadWts = numpy.array((numQuadPts, numQuadPts, numQuadPts))
+ self.basis = numpy.array((numQuadPts, numQuadPts, numQuadPts, numBasisFns, numBasisFns, numBasisFns))
+ self.basisDeriv = numpy.array((numQuadPts, numQuadPts, numQuadPts, numBasisFns, numBasisFns, numBasisFns, dim))
+ for q in range(numQuadPts):
+ for r in range(numQuadPts):
+ for s in range(numQuadPts):
+ self.quadPts[q][r][s][0] = quadpts[q]
+ self.quadPts[q][r][s][1] = quadpts[r]
+ self.quadPts[q][r][s][2] = quadpts[s]
+ self.quadWts[q][r][s] = quadwts[q]*quadwts[r]*quadwts[s]
+ for f in range(numBasisFns):
+ for g in range(numBasisFns):
+ for h in range(numBasisFns):
+ self.basis[q][r][s][f][g][h] = basis[q][f]*basis[r][g]*basis[s][h]
+ self.basisDeriv[q][r][s][f][g][h][0] = basisDeriv[q][f][0]*basis[r][g]*basis[s][h]
+ self.basisDeriv[q][r][s][f][g][h][1] = basis[q][f]*basisDeriv[r][g][0]*basis[s][h]
+ self.basisDeriv[q][r][s][f][g][h][2] = basis[q][f]*basis[r][g]*basisDeriv[s][h][0]
+ 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(self.basis)
self._info.line("Basis derivatives (quad pts):")
self._info.line(self.basisDeriv)
self._info.line("Quad pts:")
- self._info.line(quadrature.get_points())
+ self._info.line(self.quadPts)
self._info.line("Quad wts:")
- self._info.line(quadrature.get_weights())
+ self._info.line(self.quadWts)
self._info.log()
return
@@ -137,20 +178,20 @@
"""
Setup quadrature rule for reference cell.
"""
- # :TODO: ADD STUFF HERE
- return
+ import FIAT.quadrature
+ import FIAT.shapes
+ return FIAT.quadrature.make_quadrature_by_degree(FIAT.shapes.LINE, self.order)
- def _setupBasisFns(self):
+ def _setupElement(self):
"""
- Setup basis functions for reference cell.
+ Setup the finite element for reference cell.
"""
- from FIAT.Lagrange import Lagrange
+ import FIAT.Lagrange
+ import FIAT.shapes
+ return FIAT.Lagrange.Lagrange(FIAT.shapes.LINE, self.degree)
- # :TODO: ADD STUFF HERE
- return
-
# FACTORIES ////////////////////////////////////////////////////////////
def reference_cell():
Modified: short/3D/PyLith/trunk/unittests/pytests/feassemble/TestFIATLagrange.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/feassemble/TestFIATLagrange.py 2007-06-02 04:30:16 UTC (rev 7047)
+++ short/3D/PyLith/trunk/unittests/pytests/feassemble/TestFIATLagrange.py 2007-06-02 16:35:15 UTC (rev 7048)
@@ -23,37 +23,37 @@
return (1-p[0])*(1-p[1])/4.0
def N0p(p):
- return -0.25
+ return -(1-p[1])/4.0
def N0q(p):
- return -0.25
+ return -(1-p[0])/4.0
def N1(p):
- return (1+p[0])*(1-p[1])/4.0
+ return (1-p[0])*(1+p[1])/4.0
def N1p(p):
- return +0.25
+ return -(1+p[1])/4.0
def N1q(p):
- return -0.25
+ return (1-p[0])/4.0
def N2(p):
- return (1+p[0])*(1+p[1])/4.0
+ return (1+p[0])*(1-p[1])/4.0
def N2p(p):
- return +0.25
+ return (1-p[1])/4.0
def N2q(p):
- return +0.25
+ return -(1+p[0])/4.0
def N3(p):
- return (1-p[0])*(1+p[1])/4.0
+ return (1+p[0])*(1+p[1])/4.0
def N3p(p):
- return -0.25
+ return (1+p[1])/4.0
def N3q(p):
- return +0.25
+ return (1+p[0])/4.0
# ----------------------------------------------------------------------
class TestFIATLagrange(unittest.TestCase):
@@ -68,11 +68,12 @@
"""
from pylith.feassemble.FIATLagrange import FIATLagrange
cell = FIATLagrange()
+ cell.cellDim = 2
cell.degree = 1
- cell.order = 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)],
dtype=numpy.float64 )
quadWtsE = numpy.array( [1.0, 1.0, 1.0, 1.0], dtype=numpy.float64 )
More information about the cig-commits
mailing list