[cig-commits] r8934 - in cs/benchmark/cigma/trunk/src: . cigma
cigmalib
luis at geodynamics.org
luis at geodynamics.org
Wed Dec 19 12:05:08 PST 2007
Author: luis
Date: 2007-12-19 12:05:08 -0800 (Wed, 19 Dec 2007)
New Revision: 8934
Added:
cs/benchmark/cigma/trunk/src/cigmalib/
cs/benchmark/cigma/trunk/src/cigmalib/Cell.py
cs/benchmark/cigma/trunk/src/cigmalib/FiatLagrange.py
cs/benchmark/cigma/trunk/src/cigmalib/FiatQuadrature.py
cs/benchmark/cigma/trunk/src/cigmalib/FiatShapes.py
cs/benchmark/cigma/trunk/src/cigmalib/FiatSimplex.py
cs/benchmark/cigma/trunk/src/cigmalib/__init__.py
cs/benchmark/cigma/trunk/src/cigmalib/api.py
Removed:
cs/benchmark/cigma/trunk/src/cigma/Cell.py
cs/benchmark/cigma/trunk/src/cigma/FiatLagrange.py
cs/benchmark/cigma/trunk/src/cigma/FiatQuadrature.py
cs/benchmark/cigma/trunk/src/cigma/FiatShapes.py
cs/benchmark/cigma/trunk/src/cigma/FiatSimplex.py
cs/benchmark/cigma/trunk/src/cigma/__init__.py
cs/benchmark/cigma/trunk/src/cigma/api.py
Log:
Renamed FIAT wrapper python module to cigmalib
Deleted: cs/benchmark/cigma/trunk/src/cigma/Cell.py
===================================================================
--- cs/benchmark/cigma/trunk/src/cigma/Cell.py 2007-12-19 20:05:00 UTC (rev 8933)
+++ cs/benchmark/cigma/trunk/src/cigma/Cell.py 2007-12-19 20:05:08 UTC (rev 8934)
@@ -1,30 +0,0 @@
-#!/usr/bin/env python
-
-class Cell(object):
- """
- Python abstract base class for managing basis functions
- and quadrature rules of a reference finite-element cell.
- """
-
- def __init__(self, name):
-
- self.name = name.lower()
-
- self.geometry = None
- self.basisVert = None # basis fns at vertices
- self.basisDerivVert = None # basis fn derivs at vertices
- self.basisQuad = None # basis fns at quad pts
- self.basisDerivQuad = None # basis fn derivs at quad points
-
- self.quadPts = None # coordinates of quad pts
- self.quadWts = None # weights of quad pts
-
- self.cellDim = None # dimension of cell
- self.numCorners = None # number of vertices in cell
- self.numQuadPts = None # number of quadrature points
-
- return
-
- def configure(self):
- pass
-
Deleted: cs/benchmark/cigma/trunk/src/cigma/FiatLagrange.py
===================================================================
--- cs/benchmark/cigma/trunk/src/cigma/FiatLagrange.py 2007-12-19 20:05:00 UTC (rev 8933)
+++ cs/benchmark/cigma/trunk/src/cigma/FiatLagrange.py 2007-12-19 20:05:08 UTC (rev 8934)
@@ -1,636 +0,0 @@
-#!/usr/bin/env python
-
-import numpy
-from Cell import Cell
-
-class FiatLagrange(Cell):
- """
- Python object for managing basis functions and quadrature rules of a
- Lagrange reference finite-element cell using FIAT.
- """
-
- def __init__(self, name, **kw):
- Cell.__init__(self, name)
- self.cellDim = kw.get('dimension', 3) # dimension of finite element cell
- self.degree = kw.get('degree', 1) # degree of finite-element cell
- self.order = kw.get('order', -1) # order of finite-element cell
- self.configure()
- return
-
- def configure(self):
- Cell.configure(self)
- if self.order == -1:
- self.order = 2*(self.degree) + 1
- return
-
- def setupGeometry(self, spaceDim):
- self.geometry = None
- pass # XXX: cell geometry object
-
- def setupQuadrature(self):
- """
- Setup quadrature rule for reference cell.
- """
- from FIAT.quadrature import make_quadrature_by_degree
- import FIAT.shapes
- return make_quadrature_by_degree(FIAT.shapes.LINE, self.order)
-
- def setupElement(self):
- """
- Setup the finite element for reference cell
- """
- from FIAT.Lagrange import Lagrange
- import FIAT.shapes
- return Lagrange(FIAT.shapes.LINE, self.degree)
-
- def setupVertices(self, element):
- """
- Setup evaluation functional points for reference cell.
- """
- return element.Udual.pts
-
- def initialize(self, spaceDim):
- """
- Initialize reference finite-element cell from a tensor
- product of 1-D Lagrange elements.
- """
- self.setupGeometry(spaceDim)
-
- assert self.cellDim > 0
-
- quadrature = self.setupQuadrature()
- element = self.setupElement()
- dim = self.cellDim
-
- # Get coordinates of vertices (dual basis)
- vertices = numpy.array(self.setupVertices(element))
-
- # Evaluate basis functions at quadrature points
- 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
- 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
-
- if dim == 1:
- #######################################################################
- self.vertices = numpy.array(vertices)
- self.quadPts = quadpts
- self.quadWts = quadwts
- self.basis = numpy.reshape(basis.flatten(), basis.shape)
- self.basisDeriv = numpy.reshape(basisDeriv.flatten(), basisDeriv.shape)
- elif dim == 2:
- #######################################################################
- v = self.vertices = numpy.zeros((self.numCorners, dim))
- n = 0
- # Bottom
- for r in range(0, numBasisFns-1):
- v[n,0] = vertices[r]
- v[n,1] = vertices[0]
- n += 1
- # Right
- for q in range(0, numBasisFns-1):
- v[n,0] = vertices[numBasisFns-1]
- v[n,1] = vertices[q]
- n += 1
- # Top
- for r in range(numBasisFns-1, 0, -1):
- v[n,0] = vertices[r]
- v[n,1] = vertices[numBasisFns-1]
- n += 1
- # Left
- for q in range(numBasisFns-1, 0, -1):
- v[n,0] = vertices[0]
- v[n,1] = vertices[q]
- # Interior
- for q in range(1, numBasisFns-1):
- for r in range(1, numBasisFns-1):
- v[n,0] = vertices[r]
- v[n,1] = vertices[q]
- n += 1
-
- if n != self.numCorners:
- raise RuntimeError("Invalid 2D function tabulation")
-
- qp = self.quadPts = numpy.zeros((self.numQuadPts, dim))
- qw = self.quadWts = numpy.zeros((self.numQuadPts,))
- bs = self.basis = numpy.zeros((self.numQuadPts, self.numCorners))
- dbs = self.basisDeriv = numpy.zeros((self.numQuadPts, self.numCorners, dim))
-
- 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]
- 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]
- 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]
- 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]
- 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]
- 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]
- 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]
- 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]
- 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]
- m += 1
- if not m == self.numCorners:
- raise RuntimeError('Invalid 2D function tabulation')
- n += 1
- if not n == self.numQuadPts:
- raise RuntimeError('Invalid 2D quadrature')
-
- 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))
- self.basis = numpy.reshape(self.basis, (self.numQuadPts, self.numCorners))
- self.basisDeriv = numpy.reshape(self.basisDeriv, (self.numQuadPts, self.numCorners, dim))
-
- 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]
- n += 1
- # Right
- for q in range(0, numBasisFns-1):
- self.vertices[n][0] = vertices[numBasisFns-1]
- self.vertices[n][1] = vertices[q]
- self.vertices[n][2] = vertices[s]
- 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):
- self.vertices[n][0] = vertices[0]
- self.vertices[n][1] = vertices[q]
- self.vertices[n][2] = vertices[s]
- 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))
-
- 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))
- self.basisDeriv = numpy.zeros((numQuadPts*numQuadPts*numQuadPts,
- numBasisFns*numBasisFns*numBasisFns,
- dim))
- 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]
- self.quadPts[n][1] = quadpts[q]
- self.quadPts[n][2] = quadpts[s]
- self.quadWts[n] = quadwts[r]*quadwts[q]*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[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')
- n += 1
-
- if not n == self.numQuadPts:
- raise RuntimeError('Invalid 3D quadrature')
-
- 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))
- self.basis = numpy.reshape(self.basis, (self.numQuadPts, self.numCorners))
- self.basisDeriv = numpy.reshape(self.basisDeriv, (self.numQuadPts, self.numCorners, dim))
-
- return
Deleted: cs/benchmark/cigma/trunk/src/cigma/FiatQuadrature.py
===================================================================
--- cs/benchmark/cigma/trunk/src/cigma/FiatQuadrature.py 2007-12-19 20:05:00 UTC (rev 8933)
+++ cs/benchmark/cigma/trunk/src/cigma/FiatQuadrature.py 2007-12-19 20:05:08 UTC (rev 8934)
@@ -1,117 +0,0 @@
-#!/usr/bin/env python
-import FIAT.shapes
-from FIAT.quadrature import make_quadrature_by_degree
-from numpy import array, zeros
-from FiatShapes import *
-
-
-def make_rule(shape):
- def qr(order):
- q = make_quadrature_by_degree(shape, order)
- x = array(q.get_points())
- w = array(q.get_weights())
- return (x,w)
- return qr
-
-line_qr = make_rule(FIAT.shapes.LINE)
-tri_qr = make_rule(FIAT.shapes.TRIANGLE)
-tet_qr = make_rule(FIAT.shapes.TETRAHEDRON)
-
-def quad_qr(order):
- qpts,qwts = line_qr(order)
- nq = len(qpts)
- x = zeros((nq*nq, 2))
- w = zeros((nq*nq,))
- n = 0
- # Bottom
- for r in range(0,nq-1):
- x[n,0] = qpts[r]
- x[n,1] = qpts[0]
- w[n] = qwts[r]*qwts[0]
- n += 1
- # Right
- for q in range(0,nq-1):
- x[n,0] = qpts[nq-1]
- x[n,1] = qpts[q]
- w[n] = qwts[nq-1]*qwts[q]
- n += 1
- # Top
- for r in range(nq-1, 0, -1):
- x[n,0] = qpts[r]
- x[n,1] = qpts[nq-1]
- w[n] = qwts[r]*qwts[nq-1]
- n += 1
- # Left
- for q in range(nq-1, 0, -1):
- x[n,0] = qpts[0]
- x[n,1] = qpts[q]
- w[n] = qwts[0]*qwts[q]
- n += 1
- # Interior
- for q in range(1, nq-1):
- for r in range(1, nq-1):
- x[n,0] = qpts[r]
- x[n,1] = qpts[q]
- w[n] = qwts[r]*qwts[q]
- n += 1
- assert (n == nq*nq*nq)
- return (x,w)
-
-
-def hex_qr(order):
- qpts,qwts = line_qr(order)
- nq = len(qpts)
- x = zeros((nq*nq*nq, 3))
- w = zeros((nq*nq*nq,))
- n = 0
- for s in range(nq):
- # Bottom
- for r in range(0,nq-1):
- x[n,0] = qpts[r]
- x[n,1] = qpts[0]
- x[n,2] = qpts[s]
- w[n] = qwts[r]*qwts[0]*qwts[s]
- n += 1
- # Right
- for q in range(0,nq-1):
- x[n,0] = qpts[nq-1]
- x[n,1] = qpts[q]
- x[n,2] = qpts[s]
- w[n] = qwts[nq-1]*qwts[q]*qwts[s]
- n += 1
- # Top
- for r in range(nq-1, 0, -1):
- x[n,0] = qpts[r]
- x[n,1] = qpts[nq-1]
- x[n,2] = qpts[s]
- w[n] = qwts[r]*qwts[nq-1]*qwts[s]
- n += 1
- # Left
- for q in range(nq-1, 0, -1):
- x[n,0] = qpts[0]
- x[n,1] = qpts[q]
- x[n,2] = qpts[s]
- w[n] = qwts[0]*qwts[q]*qwts[s]
- n += 1
- # Interior
- for q in range(1, nq-1):
- for r in range(1, nq-1):
- x[n,0] = qpts[r]
- x[n,1] = qpts[q]
- x[n,2] = qpts[s]
- w[n] = qwts[r]*qwts[q]*qwts[s]
- n += 1
- assert (n == nq*nq*nq)
- return (x,w)
-
-
-def quadrature(shape, order):
- qr = {
- LINE: line_qr,
- TRI: tri_qr,
- TETRA: tet_qr,
- QUAD: quad_qr,
- HEX: hex_qr,
- }.get(shape)
- return qr(order)
-
Deleted: cs/benchmark/cigma/trunk/src/cigma/FiatShapes.py
===================================================================
--- cs/benchmark/cigma/trunk/src/cigma/FiatShapes.py 2007-12-19 20:05:00 UTC (rev 8933)
+++ cs/benchmark/cigma/trunk/src/cigma/FiatShapes.py 2007-12-19 20:05:08 UTC (rev 8934)
@@ -1,8 +0,0 @@
-#!/usr/bin/env python
-
-LINE = 1
-TRI = 2
-TETRA = 3
-QUAD = 4
-HEX = 5
-
Deleted: cs/benchmark/cigma/trunk/src/cigma/FiatSimplex.py
===================================================================
--- cs/benchmark/cigma/trunk/src/cigma/FiatSimplex.py 2007-12-19 20:05:00 UTC (rev 8933)
+++ cs/benchmark/cigma/trunk/src/cigma/FiatSimplex.py 2007-12-19 20:05:08 UTC (rev 8934)
@@ -1,89 +0,0 @@
-#!/usr/bin/env python
-
-import numpy
-from Cell import Cell
-
-class FiatSimplex(Cell):
- """
- Python object for managing basis functions and quadrature rules
- of a simplex reference finite-element cell using FIAT.
- """
- def __init__(self, name, **kw):
- Cell.__init__(self, name)
- self.degree = kw.get('degree', 1)
- self.order = kw.get('order', -1)
- self.configure()
- return
-
- def configure(self):
- self.shape = self.getShape()
- if self.order == -1:
- self.order = self.degree
- return
-
- def getShape(self):
- import FIAT.shapes
- shape = None
- name = self.name.lower()
- if name == "tetrahedron":
- shape = FIAT.shapes.TETRAHEDRON
- elif name == "triangle":
- shape = FIAT.shapes.TRIANGLE
- elif name == "line":
- shape == FIAT.shapes.LINE
- return shape
-
- def setupGeometry(self, spacedim):
- self.geometry = None
- pass # XXX: cell geometry object
-
- def setupQuadrature(self):
- """
- Setup quadrature rule for reference cell.
- """
- import FIAT.quadrature
- return FIAT.quadrature.make_quadrature_by_degree(self.shape, self.order)
-
- def setupBasisFns(self):
- """
- Setup basis functions for reference cell
- """
- from FIAT.Lagrange import Lagrange
- return Lagrange(self.shape, self.degree).function_space()
-
- def setupVertices(self):
- """
- Setup evaulation functional points for reference
- """
- from FIAT.Lagrange import Lagrange
- return Lagrange(self.shape, self.degree).Udual.pts
-
- def initialize(self, spaceDim):
- self.setupGeometry(spaceDim)
-
- quadrature = self.setupQuadrature()
- basisFns = self.setupBasisFns()
-
- # Get coordinates of vertices (dual basis)
- self.vertices = numpy.array(self.setupVertices(), dtype=numpy.float64)
-
- # 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)
-
- # 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)
-
- self.quadPts = numpy.array(quadrature.get_points())
- self.quadWts = numpy.array(quadrature.get_weights())
-
- self.cellDim = dim
- self.numCorners = len(basisFns)
- self.numQuadPts = len(quadrature.get_weights())
-
- return
Deleted: cs/benchmark/cigma/trunk/src/cigma/__init__.py
===================================================================
Deleted: cs/benchmark/cigma/trunk/src/cigma/api.py
===================================================================
--- cs/benchmark/cigma/trunk/src/cigma/api.py 2007-12-19 20:05:00 UTC (rev 8933)
+++ cs/benchmark/cigma/trunk/src/cigma/api.py 2007-12-19 20:05:08 UTC (rev 8934)
@@ -1,11 +0,0 @@
-import numpy
-
-def quadrature(shape, order):
- x = numpy.zeros((10,3))
- w = numpy.zeros((10,))
- return (x,w)
-
-def tabulate(shape, degree, x):
- tab = numpy.zeros((10,8))
- return tab
-
Copied: cs/benchmark/cigma/trunk/src/cigmalib/Cell.py (from rev 8933, cs/benchmark/cigma/trunk/src/cigma/Cell.py)
===================================================================
--- cs/benchmark/cigma/trunk/src/cigma/Cell.py 2007-12-19 20:05:00 UTC (rev 8933)
+++ cs/benchmark/cigma/trunk/src/cigmalib/Cell.py 2007-12-19 20:05:08 UTC (rev 8934)
@@ -0,0 +1,30 @@
+#!/usr/bin/env python
+
+class Cell(object):
+ """
+ Python abstract base class for managing basis functions
+ and quadrature rules of a reference finite-element cell.
+ """
+
+ def __init__(self, name):
+
+ self.name = name.lower()
+
+ self.geometry = None
+ self.basisVert = None # basis fns at vertices
+ self.basisDerivVert = None # basis fn derivs at vertices
+ self.basisQuad = None # basis fns at quad pts
+ self.basisDerivQuad = None # basis fn derivs at quad points
+
+ self.quadPts = None # coordinates of quad pts
+ self.quadWts = None # weights of quad pts
+
+ self.cellDim = None # dimension of cell
+ self.numCorners = None # number of vertices in cell
+ self.numQuadPts = None # number of quadrature points
+
+ return
+
+ def configure(self):
+ pass
+
Copied: cs/benchmark/cigma/trunk/src/cigmalib/FiatLagrange.py (from rev 8933, cs/benchmark/cigma/trunk/src/cigma/FiatLagrange.py)
===================================================================
--- cs/benchmark/cigma/trunk/src/cigma/FiatLagrange.py 2007-12-19 20:05:00 UTC (rev 8933)
+++ cs/benchmark/cigma/trunk/src/cigmalib/FiatLagrange.py 2007-12-19 20:05:08 UTC (rev 8934)
@@ -0,0 +1,636 @@
+#!/usr/bin/env python
+
+import numpy
+from Cell import Cell
+
+class FiatLagrange(Cell):
+ """
+ Python object for managing basis functions and quadrature rules of a
+ Lagrange reference finite-element cell using FIAT.
+ """
+
+ def __init__(self, name, **kw):
+ Cell.__init__(self, name)
+ self.cellDim = kw.get('dimension', 3) # dimension of finite element cell
+ self.degree = kw.get('degree', 1) # degree of finite-element cell
+ self.order = kw.get('order', -1) # order of finite-element cell
+ self.configure()
+ return
+
+ def configure(self):
+ Cell.configure(self)
+ if self.order == -1:
+ self.order = 2*(self.degree) + 1
+ return
+
+ def setupGeometry(self, spaceDim):
+ self.geometry = None
+ pass # XXX: cell geometry object
+
+ def setupQuadrature(self):
+ """
+ Setup quadrature rule for reference cell.
+ """
+ from FIAT.quadrature import make_quadrature_by_degree
+ import FIAT.shapes
+ return make_quadrature_by_degree(FIAT.shapes.LINE, self.order)
+
+ def setupElement(self):
+ """
+ Setup the finite element for reference cell
+ """
+ from FIAT.Lagrange import Lagrange
+ import FIAT.shapes
+ return Lagrange(FIAT.shapes.LINE, self.degree)
+
+ def setupVertices(self, element):
+ """
+ Setup evaluation functional points for reference cell.
+ """
+ return element.Udual.pts
+
+ def initialize(self, spaceDim):
+ """
+ Initialize reference finite-element cell from a tensor
+ product of 1-D Lagrange elements.
+ """
+ self.setupGeometry(spaceDim)
+
+ assert self.cellDim > 0
+
+ quadrature = self.setupQuadrature()
+ element = self.setupElement()
+ dim = self.cellDim
+
+ # Get coordinates of vertices (dual basis)
+ vertices = numpy.array(self.setupVertices(element))
+
+ # Evaluate basis functions at quadrature points
+ 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
+ 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
+
+ if dim == 1:
+ #######################################################################
+ self.vertices = numpy.array(vertices)
+ self.quadPts = quadpts
+ self.quadWts = quadwts
+ self.basis = numpy.reshape(basis.flatten(), basis.shape)
+ self.basisDeriv = numpy.reshape(basisDeriv.flatten(), basisDeriv.shape)
+ elif dim == 2:
+ #######################################################################
+ v = self.vertices = numpy.zeros((self.numCorners, dim))
+ n = 0
+ # Bottom
+ for r in range(0, numBasisFns-1):
+ v[n,0] = vertices[r]
+ v[n,1] = vertices[0]
+ n += 1
+ # Right
+ for q in range(0, numBasisFns-1):
+ v[n,0] = vertices[numBasisFns-1]
+ v[n,1] = vertices[q]
+ n += 1
+ # Top
+ for r in range(numBasisFns-1, 0, -1):
+ v[n,0] = vertices[r]
+ v[n,1] = vertices[numBasisFns-1]
+ n += 1
+ # Left
+ for q in range(numBasisFns-1, 0, -1):
+ v[n,0] = vertices[0]
+ v[n,1] = vertices[q]
+ # Interior
+ for q in range(1, numBasisFns-1):
+ for r in range(1, numBasisFns-1):
+ v[n,0] = vertices[r]
+ v[n,1] = vertices[q]
+ n += 1
+
+ if n != self.numCorners:
+ raise RuntimeError("Invalid 2D function tabulation")
+
+ qp = self.quadPts = numpy.zeros((self.numQuadPts, dim))
+ qw = self.quadWts = numpy.zeros((self.numQuadPts,))
+ bs = self.basis = numpy.zeros((self.numQuadPts, self.numCorners))
+ dbs = self.basisDeriv = numpy.zeros((self.numQuadPts, self.numCorners, dim))
+
+ 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]
+ 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]
+ 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]
+ 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]
+ 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]
+ 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]
+ 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]
+ 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]
+ 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]
+ m += 1
+ if not m == self.numCorners:
+ raise RuntimeError('Invalid 2D function tabulation')
+ n += 1
+ if not n == self.numQuadPts:
+ raise RuntimeError('Invalid 2D quadrature')
+
+ 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))
+ self.basis = numpy.reshape(self.basis, (self.numQuadPts, self.numCorners))
+ self.basisDeriv = numpy.reshape(self.basisDeriv, (self.numQuadPts, self.numCorners, dim))
+
+ 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]
+ n += 1
+ # Right
+ for q in range(0, numBasisFns-1):
+ self.vertices[n][0] = vertices[numBasisFns-1]
+ self.vertices[n][1] = vertices[q]
+ self.vertices[n][2] = vertices[s]
+ 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):
+ self.vertices[n][0] = vertices[0]
+ self.vertices[n][1] = vertices[q]
+ self.vertices[n][2] = vertices[s]
+ 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))
+
+ 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))
+ self.basisDeriv = numpy.zeros((numQuadPts*numQuadPts*numQuadPts,
+ numBasisFns*numBasisFns*numBasisFns,
+ dim))
+ 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]
+ self.quadPts[n][1] = quadpts[q]
+ self.quadPts[n][2] = quadpts[s]
+ self.quadWts[n] = quadwts[r]*quadwts[q]*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[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')
+ n += 1
+
+ if not n == self.numQuadPts:
+ raise RuntimeError('Invalid 3D quadrature')
+
+ 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))
+ self.basis = numpy.reshape(self.basis, (self.numQuadPts, self.numCorners))
+ self.basisDeriv = numpy.reshape(self.basisDeriv, (self.numQuadPts, self.numCorners, dim))
+
+ return
Copied: cs/benchmark/cigma/trunk/src/cigmalib/FiatQuadrature.py (from rev 8933, cs/benchmark/cigma/trunk/src/cigma/FiatQuadrature.py)
===================================================================
--- cs/benchmark/cigma/trunk/src/cigma/FiatQuadrature.py 2007-12-19 20:05:00 UTC (rev 8933)
+++ cs/benchmark/cigma/trunk/src/cigmalib/FiatQuadrature.py 2007-12-19 20:05:08 UTC (rev 8934)
@@ -0,0 +1,117 @@
+#!/usr/bin/env python
+import FIAT.shapes
+from FIAT.quadrature import make_quadrature_by_degree
+from numpy import array, zeros
+from FiatShapes import *
+
+
+def make_rule(shape):
+ def qr(order):
+ q = make_quadrature_by_degree(shape, order)
+ x = array(q.get_points())
+ w = array(q.get_weights())
+ return (x,w)
+ return qr
+
+line_qr = make_rule(FIAT.shapes.LINE)
+tri_qr = make_rule(FIAT.shapes.TRIANGLE)
+tet_qr = make_rule(FIAT.shapes.TETRAHEDRON)
+
+def quad_qr(order):
+ qpts,qwts = line_qr(order)
+ nq = len(qpts)
+ x = zeros((nq*nq, 2))
+ w = zeros((nq*nq,))
+ n = 0
+ # Bottom
+ for r in range(0,nq-1):
+ x[n,0] = qpts[r]
+ x[n,1] = qpts[0]
+ w[n] = qwts[r]*qwts[0]
+ n += 1
+ # Right
+ for q in range(0,nq-1):
+ x[n,0] = qpts[nq-1]
+ x[n,1] = qpts[q]
+ w[n] = qwts[nq-1]*qwts[q]
+ n += 1
+ # Top
+ for r in range(nq-1, 0, -1):
+ x[n,0] = qpts[r]
+ x[n,1] = qpts[nq-1]
+ w[n] = qwts[r]*qwts[nq-1]
+ n += 1
+ # Left
+ for q in range(nq-1, 0, -1):
+ x[n,0] = qpts[0]
+ x[n,1] = qpts[q]
+ w[n] = qwts[0]*qwts[q]
+ n += 1
+ # Interior
+ for q in range(1, nq-1):
+ for r in range(1, nq-1):
+ x[n,0] = qpts[r]
+ x[n,1] = qpts[q]
+ w[n] = qwts[r]*qwts[q]
+ n += 1
+ assert (n == nq*nq*nq)
+ return (x,w)
+
+
+def hex_qr(order):
+ qpts,qwts = line_qr(order)
+ nq = len(qpts)
+ x = zeros((nq*nq*nq, 3))
+ w = zeros((nq*nq*nq,))
+ n = 0
+ for s in range(nq):
+ # Bottom
+ for r in range(0,nq-1):
+ x[n,0] = qpts[r]
+ x[n,1] = qpts[0]
+ x[n,2] = qpts[s]
+ w[n] = qwts[r]*qwts[0]*qwts[s]
+ n += 1
+ # Right
+ for q in range(0,nq-1):
+ x[n,0] = qpts[nq-1]
+ x[n,1] = qpts[q]
+ x[n,2] = qpts[s]
+ w[n] = qwts[nq-1]*qwts[q]*qwts[s]
+ n += 1
+ # Top
+ for r in range(nq-1, 0, -1):
+ x[n,0] = qpts[r]
+ x[n,1] = qpts[nq-1]
+ x[n,2] = qpts[s]
+ w[n] = qwts[r]*qwts[nq-1]*qwts[s]
+ n += 1
+ # Left
+ for q in range(nq-1, 0, -1):
+ x[n,0] = qpts[0]
+ x[n,1] = qpts[q]
+ x[n,2] = qpts[s]
+ w[n] = qwts[0]*qwts[q]*qwts[s]
+ n += 1
+ # Interior
+ for q in range(1, nq-1):
+ for r in range(1, nq-1):
+ x[n,0] = qpts[r]
+ x[n,1] = qpts[q]
+ x[n,2] = qpts[s]
+ w[n] = qwts[r]*qwts[q]*qwts[s]
+ n += 1
+ assert (n == nq*nq*nq)
+ return (x,w)
+
+
+def quadrature(shape, order):
+ qr = {
+ LINE: line_qr,
+ TRI: tri_qr,
+ TETRA: tet_qr,
+ QUAD: quad_qr,
+ HEX: hex_qr,
+ }.get(shape)
+ return qr(order)
+
Copied: cs/benchmark/cigma/trunk/src/cigmalib/FiatShapes.py (from rev 8933, cs/benchmark/cigma/trunk/src/cigma/FiatShapes.py)
===================================================================
--- cs/benchmark/cigma/trunk/src/cigma/FiatShapes.py 2007-12-19 20:05:00 UTC (rev 8933)
+++ cs/benchmark/cigma/trunk/src/cigmalib/FiatShapes.py 2007-12-19 20:05:08 UTC (rev 8934)
@@ -0,0 +1,8 @@
+#!/usr/bin/env python
+
+LINE = 1
+TRI = 2
+TETRA = 3
+QUAD = 4
+HEX = 5
+
Copied: cs/benchmark/cigma/trunk/src/cigmalib/FiatSimplex.py (from rev 8933, cs/benchmark/cigma/trunk/src/cigma/FiatSimplex.py)
===================================================================
--- cs/benchmark/cigma/trunk/src/cigma/FiatSimplex.py 2007-12-19 20:05:00 UTC (rev 8933)
+++ cs/benchmark/cigma/trunk/src/cigmalib/FiatSimplex.py 2007-12-19 20:05:08 UTC (rev 8934)
@@ -0,0 +1,89 @@
+#!/usr/bin/env python
+
+import numpy
+from Cell import Cell
+
+class FiatSimplex(Cell):
+ """
+ Python object for managing basis functions and quadrature rules
+ of a simplex reference finite-element cell using FIAT.
+ """
+ def __init__(self, name, **kw):
+ Cell.__init__(self, name)
+ self.degree = kw.get('degree', 1)
+ self.order = kw.get('order', -1)
+ self.configure()
+ return
+
+ def configure(self):
+ self.shape = self.getShape()
+ if self.order == -1:
+ self.order = self.degree
+ return
+
+ def getShape(self):
+ import FIAT.shapes
+ shape = None
+ name = self.name.lower()
+ if name == "tetrahedron":
+ shape = FIAT.shapes.TETRAHEDRON
+ elif name == "triangle":
+ shape = FIAT.shapes.TRIANGLE
+ elif name == "line":
+ shape == FIAT.shapes.LINE
+ return shape
+
+ def setupGeometry(self, spacedim):
+ self.geometry = None
+ pass # XXX: cell geometry object
+
+ def setupQuadrature(self):
+ """
+ Setup quadrature rule for reference cell.
+ """
+ import FIAT.quadrature
+ return FIAT.quadrature.make_quadrature_by_degree(self.shape, self.order)
+
+ def setupBasisFns(self):
+ """
+ Setup basis functions for reference cell
+ """
+ from FIAT.Lagrange import Lagrange
+ return Lagrange(self.shape, self.degree).function_space()
+
+ def setupVertices(self):
+ """
+ Setup evaulation functional points for reference
+ """
+ from FIAT.Lagrange import Lagrange
+ return Lagrange(self.shape, self.degree).Udual.pts
+
+ def initialize(self, spaceDim):
+ self.setupGeometry(spaceDim)
+
+ quadrature = self.setupQuadrature()
+ basisFns = self.setupBasisFns()
+
+ # Get coordinates of vertices (dual basis)
+ self.vertices = numpy.array(self.setupVertices(), dtype=numpy.float64)
+
+ # 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)
+
+ # 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)
+
+ self.quadPts = numpy.array(quadrature.get_points())
+ self.quadWts = numpy.array(quadrature.get_weights())
+
+ self.cellDim = dim
+ self.numCorners = len(basisFns)
+ self.numQuadPts = len(quadrature.get_weights())
+
+ return
Copied: cs/benchmark/cigma/trunk/src/cigmalib/__init__.py (from rev 8933, cs/benchmark/cigma/trunk/src/cigma/__init__.py)
===================================================================
Copied: cs/benchmark/cigma/trunk/src/cigmalib/api.py (from rev 8933, cs/benchmark/cigma/trunk/src/cigma/api.py)
===================================================================
--- cs/benchmark/cigma/trunk/src/cigma/api.py 2007-12-19 20:05:00 UTC (rev 8933)
+++ cs/benchmark/cigma/trunk/src/cigmalib/api.py 2007-12-19 20:05:08 UTC (rev 8934)
@@ -0,0 +1,13 @@
+import numpy
+
+def quadrature(shape, order):
+ x = numpy.zeros((10,3))
+ w = numpy.zeros((10,))
+ return (x,w)
+
+def tabulate(shape, degree, x):
+ tab = numpy.zeros((10,8))
+ return tab
+
+def test(a,b):
+ return 42
More information about the cig-commits
mailing list