[cig-commits] r19018 - in short/3D/PyLith/branches/v1.6-stable: pylith/feassemble unittests/pytests/feassemble
brad at geodynamics.org
brad at geodynamics.org
Wed Oct 5 16:17:28 PDT 2011
Author: brad
Date: 2011-10-05 16:17:27 -0700 (Wed, 05 Oct 2011)
New Revision: 19018
Modified:
short/3D/PyLith/branches/v1.6-stable/pylith/feassemble/FIATLagrange.py
short/3D/PyLith/branches/v1.6-stable/unittests/pytests/feassemble/TestFIATLagrange.py
Log:
Added ability to collocate quadrature points with vertices.
Modified: short/3D/PyLith/branches/v1.6-stable/pylith/feassemble/FIATLagrange.py
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/pylith/feassemble/FIATLagrange.py 2011-10-05 22:05:41 UTC (rev 19017)
+++ short/3D/PyLith/branches/v1.6-stable/pylith/feassemble/FIATLagrange.py 2011-10-05 23:17:27 UTC (rev 19018)
@@ -35,6 +35,23 @@
raise ValueError("Dimension of Lagrange element must be 1, 2, or 3.")
return dim
+
+from FIAT.quadrature import QuadratureRule
+class CollocatedQuadratureLineRule(QuadratureRule):
+ """
+ Quadrature points colocated with vertices.
+ """
+ def __init__(self, ref_el, m):
+ from FIAT.lagrange import Lagrange
+ vertices = Lagrange(ref_el, m).dual.get_nodes()
+ pts = [v.get_point_dict().keys()[0] for v in vertices]
+ npts = len(pts)
+ wts = (ref_el.volume()/npts,)*npts
+
+ QuadratureRule.__init__(self, ref_el, pts, wts)
+ return
+
+
# FIATLagrange class
class FIATLagrange(ReferenceCell):
"""
@@ -53,9 +70,10 @@
## Python object for managing FIATLagrange facilities and properties.
##
## \b Properties
- ## @li \b dimension Dimension of finite-element cell
- ## @li \b degree Degree of finite-element cell
- ## @li \b quad_order Order of quadrature rule
+ ## @li \b dimension Dimension of finite-element cell.
+ ## @li \b degree Degree of finite-element cell.
+ ## @li \b quad_order Order of quadrature rule.
+ ## @li \b collocate_quad Collocate quadrature points with vertices.
##
## \b Facilities
## @li None
@@ -72,7 +90,10 @@
order = pyre.inventory.int("quad_order", default=-1)
order.meta['tip'] = "Order of quadrature rule."
+ collocateQuad = pyre.inventory.bool("collocate_quad", default=False)
+ collocateQuad.meta['tip'] = "Collocate quadrature points with vertices."
+
# PUBLIC METHODS /////////////////////////////////////////////////////
def __init__(self, name="fiatlagrange"):
@@ -411,6 +432,7 @@
self.cellDim = self.inventory.dimension
self.degree = self.inventory.degree
self.order = self.inventory.order
+ self.collocateQuad = self.inventory.collocateQuad
if self.order == -1:
self.order = self.degree+1
@@ -462,9 +484,15 @@
"""
from FIAT.quadrature import make_quadrature
from FIAT.reference_element import default_simplex
- return make_quadrature(default_simplex(1), self.order)
+ if not self.collocateQuad:
+ q = make_quadrature(default_simplex(1), self.order)
+ else:
+ q = CollocatedQuadratureLineRule(default_simplex(1), self.order)
+ return q
+
+
def _setupElement(self):
"""
Setup the finite element for reference cell.
Modified: short/3D/PyLith/branches/v1.6-stable/unittests/pytests/feassemble/TestFIATLagrange.py
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/unittests/pytests/feassemble/TestFIATLagrange.py 2011-10-05 22:05:41 UTC (rev 19017)
+++ short/3D/PyLith/branches/v1.6-stable/unittests/pytests/feassemble/TestFIATLagrange.py 2011-10-05 23:17:27 UTC (rev 19018)
@@ -74,6 +74,40 @@
# ----------------------------------------------------------------------
+class Line2Collocated(Line2):
+
+ def __init__(self):
+ """
+ Setup line2 cell.
+ """
+ vertices = numpy.array([[-1.0], [1.0]])
+ quadPts = vertices[:]
+ quadWts = numpy.array( [1.0,1.0], dtype=numpy.float64 )
+
+ # Compute basis fns and derivatives at quadrature points
+ basis = numpy.zeros( (2, 2), dtype=numpy.float64)
+ basisDeriv = numpy.zeros( (2, 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
+
+ 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
+
+
+# ----------------------------------------------------------------------
class Line3(object):
def __init__(self):
@@ -207,6 +241,49 @@
# ----------------------------------------------------------------------
+class Quad4Collocated(Quad4):
+
+ 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, -1.0],
+ [+1.0, -1.0],
+ [-1.0, +1.0],
+ [+1.0, +1.0]])
+ quadWts = numpy.array( [1.0, 1.0, 1.0, 1.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
+
+ 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
+
+
+# ----------------------------------------------------------------------
class Quad9(object):
def __init__(self):
@@ -1014,6 +1091,27 @@
return
+ def test_initialize_line2_collocated(self):
+ """
+ Test initialize() with line2 cell.
+ """
+ cell = FIATLagrange()
+ cell.inventory.dimension = 1
+ cell.inventory.degree = 1
+ cell.inventory.order = 1
+ cell.inventory.collocateQuad = True
+ cell._configure()
+ cell._info.activate()
+ cell.initialize(spaceDim=1)
+
+
+ cellE = Line2Collocated()
+ self._checkVals(cellE, cell)
+ from pylith.feassemble.CellGeometry import GeometryLine1D
+ self.failUnless(isinstance(cell.geometry, GeometryLine1D))
+ return
+
+
def test_initialize_line3(self):
"""
Test initialize() with line3 cell.
@@ -1050,6 +1148,25 @@
return
+ def test_initialize_quad4_collocated(self):
+ """
+ Test initialize() with quad4 cell.
+ """
+ cell = FIATLagrange()
+ cell.inventory.dimension = 2
+ cell.inventory.degree = 1
+ cell.inventory.order = 1
+ cell.inventory.collocateQuad = True
+ cell._configure()
+ cell.initialize(spaceDim=2)
+
+ cellE = Quad4Collocated()
+ self._checkVals(cellE, cell)
+ from pylith.feassemble.CellGeometry import GeometryQuad2D
+ self.failUnless(isinstance(cell.geometry, GeometryQuad2D))
+ return
+
+
def test_initialize_quad9(self):
"""
Test initialize() with quad9 cell.
More information about the CIG-COMMITS
mailing list