[cig-commits] r14086 - in cs/cigma/trunk: . examples examples/citcomcu examples/elements examples/fiat examples/gale examples/laplace1 examples/meshes examples/meshes/tetgen examples/meshes/triangle examples/test examples/utils/tetgen examples/utils/triangle utils
luis at geodynamics.org
luis at geodynamics.org
Wed Feb 18 08:14:42 PST 2009
Author: luis
Date: 2009-02-18 08:14:41 -0800 (Wed, 18 Feb 2009)
New Revision: 14086
Added:
cs/cigma/trunk/examples/citcomcu/
cs/cigma/trunk/examples/citcomcu/README
cs/cigma/trunk/examples/elements/
cs/cigma/trunk/examples/elements/README
cs/cigma/trunk/examples/elements/check-elements.sh
cs/cigma/trunk/examples/elements/reference-cells.py
cs/cigma/trunk/examples/fiat/
cs/cigma/trunk/examples/fiat/FiatCell.py
cs/cigma/trunk/examples/fiat/FiatLagrange.py
cs/cigma/trunk/examples/fiat/FiatQuadrature.py
cs/cigma/trunk/examples/fiat/FiatShapes.py
cs/cigma/trunk/examples/fiat/FiatSimplex.py
cs/cigma/trunk/examples/fiat/integration-rules.py
cs/cigma/trunk/examples/gale/
cs/cigma/trunk/examples/gale/make-vtk-files.sh
cs/cigma/trunk/examples/laplace1/compare-cube.sh
cs/cigma/trunk/examples/laplace1/compare-square.sh
cs/cigma/trunk/examples/laplace1/make-res-vtk.sh
cs/cigma/trunk/examples/meshes/
cs/cigma/trunk/examples/meshes/tetgen/
cs/cigma/trunk/examples/meshes/tetgen/README
cs/cigma/trunk/examples/meshes/tetgen/cube.py
cs/cigma/trunk/examples/meshes/tetgen/cube.sh
cs/cigma/trunk/examples/meshes/tetgen/cube1.poly
cs/cigma/trunk/examples/meshes/tetgen/cube2.poly
cs/cigma/trunk/examples/meshes/triangle/
cs/cigma/trunk/examples/meshes/triangle/README
cs/cigma/trunk/examples/meshes/triangle/square.py
cs/cigma/trunk/examples/meshes/triangle/square.sh
cs/cigma/trunk/examples/meshes/triangle/square1.poly
cs/cigma/trunk/examples/meshes/triangle/square2.poly
cs/cigma/trunk/examples/test/
cs/cigma/trunk/examples/test/compare-cube.sh
cs/cigma/trunk/examples/test/compare-square.sh
cs/cigma/trunk/examples/test/exact-integrals.txt
cs/cigma/trunk/utils/
cs/cigma/trunk/utils/README
cs/cigma/trunk/utils/create-mesh-from-tetgen-files.py
cs/cigma/trunk/utils/create-tetgen-vol-from-mesh-volume.py
cs/cigma/trunk/utils/power-plot.py
cs/cigma/trunk/utils/remove-hdf5-node.py
Removed:
cs/cigma/trunk/examples/utils/tetgen/README
cs/cigma/trunk/examples/utils/tetgen/cube.sh
cs/cigma/trunk/examples/utils/tetgen/cube1.poly
cs/cigma/trunk/examples/utils/tetgen/cube2.poly
cs/cigma/trunk/examples/utils/triangle/README
cs/cigma/trunk/examples/utils/triangle/square.py
cs/cigma/trunk/examples/utils/triangle/square.sh
cs/cigma/trunk/examples/utils/triangle/square1.poly
cs/cigma/trunk/examples/utils/triangle/square2.poly
Log:
Updated examples and utils
Added: cs/cigma/trunk/examples/citcomcu/README
===================================================================
Added: cs/cigma/trunk/examples/elements/README
===================================================================
Added: cs/cigma/trunk/examples/elements/check-elements.sh
===================================================================
--- cs/cigma/trunk/examples/elements/check-elements.sh (rev 0)
+++ cs/cigma/trunk/examples/elements/check-elements.sh 2009-02-18 16:14:41 UTC (rev 14086)
@@ -0,0 +1,27 @@
+#!/bin/bash
+
+# show list of elements
+cigma element-info
+
+# show information on element
+cigma element-info hex8
+cigma element-info tet4
+cigma element-info quad4
+cigma element-info tri3
+
+#
+cigma mesh-info reference-cells.h5:/hex8/mesh
+cigma mesh-info reference-cells.h5:/hex8/mesh --query-point 0.5,0.5
+
+# for N0, value at (0,0,0) should be 1, and at the other vertices 0
+cigma function-info reference-cells.h5:/tet4/dofs/N0 --query-point 0,0,0
+cigma function-info reference-cells.h5:/tet4/dofs/N0 --query-point 1,0,0
+cigma function-info reference-cells.h5:/tet4/dofs/N0 --query-point 0,1,0
+cigma function-info reference-cells.h5:/tet4/dofs/N0 --query-point 0,0,1
+
+# value should always be one, regardless of which point we pick in the cell
+cigma function-info reference-cells.h5:/tri3/dofs/one --query-point 0.1,0.2
+cigma function-info reference-cells.h5:/quad4/dofs/one --query-point 0.1,0.2
+cigma function-info reference-cells.h5:/tet4/dofs/one --query-point 0.1,0.2,0.3
+cigma function-info reference-cells.h5:/hex8/dofs/one --query-point 0.1,0.2,0.3
+
Added: cs/cigma/trunk/examples/elements/reference-cells.py
===================================================================
--- cs/cigma/trunk/examples/elements/reference-cells.py (rev 0)
+++ cs/cigma/trunk/examples/elements/reference-cells.py 2009-02-18 16:14:41 UTC (rev 14086)
@@ -0,0 +1,113 @@
+#!/usr/bin/env python
+
+import numpy
+import tables
+
+cells = ['tri3', 'tri6', 'quad4', 'tet4', 'tet10', 'hex8']
+
+celldim = {
+ 'tri3' : 2,
+ 'tri6' : 2,
+ 'quad4': 2,
+ 'tet4' : 3,
+ 'tet10': 3,
+ 'hex8' : 3,
+}
+
+coords = {
+ 'tri3': [[0.0, 0.0],
+ [1.0, 0.0],
+ [0.0, 1.0]],
+
+ 'tri6': [[0.0, 0.0],
+ [1.0, 0.0],
+ [0.0, 1.0],
+ [0.0, 1.0],
+ [0.5, 0.5],
+ [0.0, 0.5],
+ [0.5, 0.0]],
+
+ 'quad4': [[-1.0, -1.0],
+ [+1.0, -1.0],
+ [+1.0, +1.0],
+ [-1.0, +1.0]],
+
+ 'tet4': [[0.0, 0.0, 0.0],
+ [1.0, 0.0, 0.0],
+ [0.0, 1.0, 0.0],
+ [0.0, 0.0, 1.0]],
+
+ 'tet10': [[0.0, 0.0, 0.0],
+ [1.0, 0.0, 0.0],
+ [0.0, 1.0, 0.0],
+ [0.0, 0.0, 1.0],
+ [0.5, 0.5, 0.0],
+ [0.0, 0.5, 0.5],
+ [0.5, 0.0, 0.5],
+ [0.0, 0.0, 0.5],
+ [0.0, 0.5, 0.0],
+ [0.5, 0.0, 0.0]],
+
+ 'hex8': [[-1.0, -1.0, -1.0],
+ [+1.0, -1.0, -1.0],
+ [+1.0, +1.0, -1.0],
+ [-1.0, +1.0, -1.0],
+ [-1.0, -1.0, +1.0],
+ [+1.0, -1.0, +1.0],
+ [+1.0, +1.0, +1.0],
+ [-1.0, +1.0, +1.0]]
+}
+
+connect = {
+ 'tri3' : [range(3)],
+ 'tri6' : [range(6)],
+ 'quad4': [range(4)],
+ 'tet4' : [range(4)],
+ 'tet10': [range(10)],
+ 'hex8' : [range(8)],
+}
+
+def make_shape(nno):
+ def N(i):
+ n = numpy.zeros(nno)
+ n[i] = 1
+ return n
+ return N
+
+def main():
+ """Create a file containing reference cell meshes"""
+ h5 = tables.openFile('reference-cells.h5', 'w')
+
+ for cell in cells:
+
+ cellgroup = h5.createGroup('/', cell)
+
+ mesh = h5.createGroup(cellgroup, 'mesh')
+ mesh._v_attrs.CellType = cell
+
+ coordinates = numpy.array(coords[cell], dtype=float)
+ connectivity = numpy.array(connect[cell], dtype=int)
+ h5.createArray(mesh, 'coordinates', coordinates)
+ h5.createArray(mesh, 'connectivity', connectivity)
+
+ dofs = h5.createGroup(cellgroup, 'dofs')
+
+ nno = coordinates.shape[0]
+ meshloc = '/%s/mesh' % cell
+
+ N = make_shape(nno)
+ for i in xrange(nno):
+ d = h5.createArray(dofs, 'N%d' % i, N(i))
+ d._v_attrs.MeshLocation = meshloc
+
+ d = h5.createArray(dofs, 'one', numpy.ones(nno))
+ d._v_attrs.MeshLocation = meshloc
+
+
+
+ print 'Wrote', h5.filename
+ h5.close()
+
+if __name__ == '__main__':
+ main()
+
Property changes on: cs/cigma/trunk/examples/elements/reference-cells.py
___________________________________________________________________
Name: svn:executable
+ *
Added: cs/cigma/trunk/examples/fiat/FiatCell.py
===================================================================
--- cs/cigma/trunk/examples/fiat/FiatCell.py (rev 0)
+++ cs/cigma/trunk/examples/fiat/FiatCell.py 2009-02-18 16:14:41 UTC (rev 14086)
@@ -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
+
Added: cs/cigma/trunk/examples/fiat/FiatLagrange.py
===================================================================
--- cs/cigma/trunk/examples/fiat/FiatLagrange.py (rev 0)
+++ cs/cigma/trunk/examples/fiat/FiatLagrange.py 2009-02-18 16:14:41 UTC (rev 14086)
@@ -0,0 +1,637 @@
+#!/usr/bin/env python
+
+import numpy
+
+from FiatCell 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
Added: cs/cigma/trunk/examples/fiat/FiatQuadrature.py
===================================================================
--- cs/cigma/trunk/examples/fiat/FiatQuadrature.py (rev 0)
+++ cs/cigma/trunk/examples/fiat/FiatQuadrature.py 2009-02-18 16:14:41 UTC (rev 14086)
@@ -0,0 +1,118 @@
+#!/usr/bin/env python
+
+import FIAT.shapes
+
+from numpy import array, zeros
+from FIAT.quadrature import make_quadrature_by_degree
+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)
+ 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,
+ TRIANGLE: tri_qr,
+ QUADRANGLE: quad_qr,
+ TETRAHEDRON: tet_qr,
+ HEXAHEDRON: hex_qr,
+ }.get(shape)
+ return qr(order)
+
Added: cs/cigma/trunk/examples/fiat/FiatShapes.py
===================================================================
--- cs/cigma/trunk/examples/fiat/FiatShapes.py (rev 0)
+++ cs/cigma/trunk/examples/fiat/FiatShapes.py 2009-02-18 16:14:41 UTC (rev 14086)
@@ -0,0 +1,8 @@
+#!/usr/bin/env python
+
+LINE = 12
+TRIANGLE = 23
+QUADRANGLE = 24
+TETRAHEDRON = 34
+HEXAHEDRON = 38
+
Added: cs/cigma/trunk/examples/fiat/FiatSimplex.py
===================================================================
--- cs/cigma/trunk/examples/fiat/FiatSimplex.py (rev 0)
+++ cs/cigma/trunk/examples/fiat/FiatSimplex.py 2009-02-18 16:14:41 UTC (rev 14086)
@@ -0,0 +1,89 @@
+#!/usr/bin/env python
+
+import numpy
+from FiatCell 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
Added: cs/cigma/trunk/examples/fiat/integration-rules.py
===================================================================
--- cs/cigma/trunk/examples/fiat/integration-rules.py (rev 0)
+++ cs/cigma/trunk/examples/fiat/integration-rules.py 2009-02-18 16:14:41 UTC (rev 14086)
@@ -0,0 +1,53 @@
+#!/usr/bin/env python
+"""
+Write out a number of quadrature rules into an HDF5 file
+"""
+
+import tables
+from FiatQuadrature import *
+
+
+def transform_rule(qr):
+ def new_qr(order):
+ x,w = qr(order)
+ x += 1
+ x *= 0.5
+ w *= 0.5 ** x.shape[1]
+ return (x,w)
+ return new_qr
+
+fp = tables.openFile("integration-rules.h5", "w")
+
+cells = ['line', 'tri3', 'quad4', 'tet4', 'hex8']
+
+orders_by_cell = {
+ 'line' : range(0,31),
+ 'tri3' : range(0,31),
+ 'quad4': range(2,31),
+ 'tet4' : range(0,31),
+ 'hex8' : range(2,31),
+}
+
+rules_by_cell = {
+ 'line' : line_qr,
+ 'tri3' : transform_rule(tri_qr),
+ 'quad4': quad_qr,
+ 'tet4' : transform_rule(tet_qr),
+ 'hex8' : hex_qr,
+}
+
+for cell in cells:
+ rule = rules_by_cell[cell]
+ orders = orders_by_cell[cell]
+ cellGroup = fp.createGroup('/', cell)
+ for order in orders:
+ x,w = rule(order)
+ o = "order%d" % order
+ g = fp.createGroup(cellGroup, o)
+ g._v_attrs.CellType = cell
+ fp.createArray(g, 'points', x)
+ fp.createArray(g, 'weights', w)
+ print "Created /%s/%s" % (cell,o)
+
+fp.close()
+
Property changes on: cs/cigma/trunk/examples/fiat/integration-rules.py
___________________________________________________________________
Name: svn:executable
+ *
Added: cs/cigma/trunk/examples/gale/make-vtk-files.sh
===================================================================
--- cs/cigma/trunk/examples/gale/make-vtk-files.sh (rev 0)
+++ cs/cigma/trunk/examples/gale/make-vtk-files.sh 2009-02-18 16:14:41 UTC (rev 14086)
@@ -0,0 +1,10 @@
+#!/bin/bash
+set -x
+
+vtk-residuals -m p256.vts -i circular_inclusion.h5:/pressure_256_064 -o circular_inclusion_256_064.vtk
+vtk-residuals -m p256.vts -i circular_inclusion.h5:/pressure_256_128 -o circular_inclusion_256_128.vtk
+
+vtk-residuals -m p64.vts -i circular_inclusion.h5:/pressure_064 -o circular_inclusion_064.vtk
+vtk-residuals -m p128.vts -i circular_inclusion.h5:/pressure_128 -o circular_inclusion_128.vtk
+vtk-residuals -m p256.vts -i circular_inclusion.h5:/pressure_256 -o circular_inclusion_256.vtk
+
Property changes on: cs/cigma/trunk/examples/gale/make-vtk-files.sh
___________________________________________________________________
Name: svn:executable
+ *
Added: cs/cigma/trunk/examples/laplace1/compare-cube.sh
===================================================================
--- cs/cigma/trunk/examples/laplace1/compare-cube.sh (rev 0)
+++ cs/cigma/trunk/examples/laplace1/compare-cube.sh 2009-02-18 16:14:41 UTC (rev 14086)
@@ -0,0 +1,8 @@
+#!/bin/bash
+set -x
+
+cigma compare cube6.vtk cube5.vtk -o res.h5:/cube_6_5 -v
+cigma compare cube6.vtk cube4.vtk -o res.h5:/cube_6_4 -v
+cigma compare cube6.vtk cube3.vtk -o res.h5:/cube_6_3 -v
+cigma compare cube6.vtk cube2.vtk -o res.h5:/cube_6_2 -v
+
Property changes on: cs/cigma/trunk/examples/laplace1/compare-cube.sh
___________________________________________________________________
Name: svn:executable
+ *
Added: cs/cigma/trunk/examples/laplace1/compare-square.sh
===================================================================
--- cs/cigma/trunk/examples/laplace1/compare-square.sh (rev 0)
+++ cs/cigma/trunk/examples/laplace1/compare-square.sh 2009-02-18 16:14:41 UTC (rev 14086)
@@ -0,0 +1,8 @@
+#!/bin/bash
+set -x
+
+cigma compare square6.vtk square5.vtk -o res.h5:/square_6_5 -v
+cigma compare square6.vtk square4.vtk -o res.h5:/square_6_4 -v
+cigma compare square6.vtk square3.vtk -o res.h5:/square_6_3 -v
+cigma compare square6.vtk square2.vtk -o res.h5:/square_6_2 -v
+
Property changes on: cs/cigma/trunk/examples/laplace1/compare-square.sh
___________________________________________________________________
Name: svn:executable
+ *
Added: cs/cigma/trunk/examples/laplace1/make-res-vtk.sh
===================================================================
--- cs/cigma/trunk/examples/laplace1/make-res-vtk.sh (rev 0)
+++ cs/cigma/trunk/examples/laplace1/make-res-vtk.sh 2009-02-18 16:14:41 UTC (rev 14086)
@@ -0,0 +1,16 @@
+#!/bin/bash
+set -x
+
+#/cube_6_2 Dataset {262144, 1}
+#/cube_6_3 Dataset {262144, 1}
+#/cube_6_4 Dataset {262144, 1}
+#/cube_6_5 Dataset {262144, 1}
+#/square_6_2 Dataset {4096, 1}
+#/square_6_3 Dataset {4096, 1}
+#/square_6_4 Dataset {4096, 1}
+#/square_6_5 Dataset {4096, 1}
+
+for i in 2 3 4 5; do
+ vtk-residuals -m cube6.vtk -i res.h5:/cube_6_${i} -o res_cube_6_${i}.vtk
+ vtk-residuals -m square6.vtk -i res.h5:/square_6_${i} -o res_square_6_${i}.vtk
+done
Property changes on: cs/cigma/trunk/examples/laplace1/make-res-vtk.sh
___________________________________________________________________
Name: svn:executable
+ *
Copied: cs/cigma/trunk/examples/meshes/tetgen/README (from rev 14085, cs/cigma/trunk/examples/utils/tetgen/README)
===================================================================
--- cs/cigma/trunk/examples/meshes/tetgen/README (rev 0)
+++ cs/cigma/trunk/examples/meshes/tetgen/README 2009-02-18 16:14:41 UTC (rev 14086)
@@ -0,0 +1,4 @@
+
+TetGen: A Quality Tetrahedral Mesh Generator and Three-Dimensional Delaunay Triangulator
+http://tetgen.berlios.de/
+
Copied: cs/cigma/trunk/examples/meshes/tetgen/cube.py (from rev 14085, cs/cigma/trunk/examples/utils/triangle/square.py)
===================================================================
--- cs/cigma/trunk/examples/meshes/tetgen/cube.py (rev 0)
+++ cs/cigma/trunk/examples/meshes/tetgen/cube.py 2009-02-18 16:14:41 UTC (rev 14086)
@@ -0,0 +1,67 @@
+#!/usr/bin/env python
+
+from __future__ import with_statement
+
+import numpy
+import tables
+
+def read_coords(filename):
+ coords = None
+ with open(filename, 'r') as fp:
+ # read header
+ line = fp.readline()
+ cols = line.split()
+ nno, ndim, natt, nbdr = map(int, cols)
+ coords = numpy.zeros((nno, ndim), dtype=float)
+ # read nodes
+ for n in xrange(nno):
+ line = fp.readline()
+ cols = line.split()
+ i = int(cols[0]) - 1
+ for j in xrange(ndim):
+ coords[i,j] = float(cols[1+j])
+ return coords
+
+def read_connect(filename):
+ connect = None
+ with open(filename, 'r') as fp:
+ # read header
+ line = fp.readline()
+ cols = line.split()
+ nel, ndofs, natt = map(int, cols)
+ connect = numpy.zeros((nel, ndofs), dtype=int)
+ # read elements
+ for e in xrange(nel):
+ line = fp.readline()
+ cols = line.split()
+ i = int(cols[0]) - 1
+ for j in xrange(ndofs):
+ connect[i,j] = int(cols[1+j])
+ connect -= 1
+ return connect
+
+def main():
+ h5 = tables.openFile('cube.h5', 'w')
+
+ levels = range(1,7)
+ celltypes = {'cube1': 'tet4', 'cube2': 'tet10'}
+
+ for prefix in ('cube1', 'cube2'):
+ root = h5.createGroup('/', celltypes[prefix])
+ for r in levels:
+ coords = read_coords('%s.%d.node' % (prefix, r))
+ connect = read_connect('%s.%d.ele' % (prefix, r))
+ group = h5.createGroup(root, 'mesh%d' % r)
+ h5.createArray(group, 'coordinates', coords)
+ h5.createArray(group, 'connectivity', connect)
+ group._v_attrs.CellType = celltypes[prefix]
+
+ print 'Wrote', h5.filename
+ h5.close()
+ return 0
+
+if __name__ == '__main__':
+ import sys
+ ret = main()
+ sys.exit(ret)
+
Added: cs/cigma/trunk/examples/meshes/tetgen/cube.sh
===================================================================
--- cs/cigma/trunk/examples/meshes/tetgen/cube.sh (rev 0)
+++ cs/cigma/trunk/examples/meshes/tetgen/cube.sh 2009-02-18 16:14:41 UTC (rev 14086)
@@ -0,0 +1,23 @@
+#!/bin/bash
+
+# refine cube1
+o="-fen"
+./tetgen $o -p cube1
+./tetgen $o -r -a0.5000 cube1.1
+./tetgen $o -r -a0.2894 cube1.2
+./tetgen $o -r -a0.1674 cube1.3
+./tetgen $o -r -a0.0969 cube1.4
+./tetgen $o -r -a0.0560 cube1.5
+./tetgen $o -r -a0.0187 cube1.6
+./tetgen $o -r -a0.0108 cube1.7
+
+# refine cube2
+o="-fen -o2"
+./tetgen $o -p cube2
+./tetgen $o -r -a0.5000 cube1.1
+./tetgen $o -r -a0.2894 cube1.2
+./tetgen $o -r -a0.1674 cube1.3
+./tetgen $o -r -a0.0969 cube1.4
+./tetgen $o -r -a0.0560 cube1.5
+./tetgen $o -r -a0.0187 cube1.6
+./tetgen $o -r -a0.0108 cube1.7
Property changes on: cs/cigma/trunk/examples/meshes/tetgen/cube.sh
___________________________________________________________________
Name: svn:executable
+ *
Copied: cs/cigma/trunk/examples/meshes/tetgen/cube1.poly (from rev 14085, cs/cigma/trunk/examples/utils/tetgen/cube1.poly)
===================================================================
--- cs/cigma/trunk/examples/meshes/tetgen/cube1.poly (rev 0)
+++ cs/cigma/trunk/examples/meshes/tetgen/cube1.poly 2009-02-18 16:14:41 UTC (rev 14086)
@@ -0,0 +1,36 @@
+# Part 1 - node list
+# node count, 3 dim, no attribute, no boundary marker
+8 3 0 0
+# Node index, node coordinates
+1 -1.0 -1.0 -1.0
+2 +1.0 -1.0 -1.0
+3 +1.0 +1.0 -1.0
+4 -1.0 +1.0 -1.0
+5 -1.0 -1.0 +1.0
+6 +1.0 -1.0 +1.0
+7 +1.0 +1.0 +1.0
+8 -1.0 +1.0 +1.0
+
+# Part 2 - facet list
+# facet count, no boundary marker
+6 0
+# facets
+1 # 1 polygon, no hole, no boundary marker
+4 1 2 3 4 # front
+1
+4 5 6 7 8 # back
+1
+4 1 2 6 5 # bottom
+1
+4 2 3 7 6 # right
+1
+4 3 4 8 7 # top
+1
+4 4 1 5 8 # left
+
+# Part 3 - hole list
+0 # no hole
+
+# Part 4 - region list
+0 # no region
+
Copied: cs/cigma/trunk/examples/meshes/tetgen/cube2.poly (from rev 14085, cs/cigma/trunk/examples/utils/tetgen/cube2.poly)
===================================================================
--- cs/cigma/trunk/examples/meshes/tetgen/cube2.poly (rev 0)
+++ cs/cigma/trunk/examples/meshes/tetgen/cube2.poly 2009-02-18 16:14:41 UTC (rev 14086)
@@ -0,0 +1,36 @@
+# Part 1 - node list
+# node count, 3 dim, no attribute, no boundary marker
+8 3 0 0
+# Node index, node coordinates
+1 -1.0 -1.0 -1.0
+2 +1.0 -1.0 -1.0
+3 +1.0 +1.0 -1.0
+4 -1.0 +1.0 -1.0
+5 -1.0 -1.0 +1.0
+6 +1.0 -1.0 +1.0
+7 +1.0 +1.0 +1.0
+8 -1.0 +1.0 +1.0
+
+# Part 2 - facet list
+# facet count, no boundary marker
+6 0
+# facets
+1 # 1 polygon, no hole, no boundary marker
+4 1 2 3 4 # front
+1
+4 5 6 7 8 # back
+1
+4 1 2 6 5 # bottom
+1
+4 2 3 7 6 # right
+1
+4 3 4 8 7 # top
+1
+4 4 1 5 8 # left
+
+# Part 3 - hole list
+0 # no hole
+
+# Part 4 - region list
+0 # no region
+
Copied: cs/cigma/trunk/examples/meshes/triangle/README (from rev 14085, cs/cigma/trunk/examples/utils/triangle/README)
===================================================================
--- cs/cigma/trunk/examples/meshes/triangle/README (rev 0)
+++ cs/cigma/trunk/examples/meshes/triangle/README 2009-02-18 16:14:41 UTC (rev 14086)
@@ -0,0 +1,4 @@
+
+Triangle: A Two-Dimensional Quality Mesh Generator and Delaunay Triangulator
+http://www.cs.cmu.edu/~quake/triangle.html
+
Copied: cs/cigma/trunk/examples/meshes/triangle/square.py (from rev 14085, cs/cigma/trunk/examples/utils/triangle/square.py)
===================================================================
--- cs/cigma/trunk/examples/meshes/triangle/square.py (rev 0)
+++ cs/cigma/trunk/examples/meshes/triangle/square.py 2009-02-18 16:14:41 UTC (rev 14086)
@@ -0,0 +1,69 @@
+#!/usr/bin/env python
+
+from __future__ import with_statement
+
+import numpy
+import tables
+
+def read_coords(filename):
+ coords = None
+ with open(filename, 'r') as fp:
+ # read header
+ line = fp.readline()
+ cols = line.split()
+ nno, ndim, natt, nbdr = map(int, cols)
+ coords = numpy.zeros((nno, ndim), dtype=float)
+ # read nodes
+ for n in xrange(nno):
+ line = fp.readline()
+ cols = line.split()
+ i = int(cols[0]) - 1
+ for j in xrange(ndim):
+ coords[i,j] = float(cols[1+j])
+ return coords
+
+def read_connect(filename):
+ connect = None
+ with open(filename, 'r') as fp:
+ # read header
+ line = fp.readline()
+ cols = line.split()
+ nel, ndofs, natt = map(int, cols)
+ connect = numpy.zeros((nel, ndofs), dtype=int)
+ # read elements
+ for e in xrange(nel):
+ line = fp.readline()
+ cols = line.split()
+ i = int(cols[0]) - 1
+ for j in xrange(ndofs):
+ connect[i,j] = int(cols[1+j])
+ connect -= 1
+ return connect
+
+def main():
+
+ h5 = tables.openFile('square.h5', 'w')
+
+ levels = range(1,7);
+ celltypes = {'square1': 'tri3', 'square2': 'tri6'}
+
+ for prefix in ('square1', 'square2'):
+ root = h5.createGroup('/', celltypes[prefix])
+ for r in levels:
+ coords = read_coords('%s.%d.node' % (prefix, r))
+ connect = read_connect('%s.%d.ele' % (prefix, r))
+ group = h5.createGroup(root, 'mesh%d' % r)
+ h5.createArray(group, 'coordinates', coords)
+ h5.createArray(group, 'connectivity', connect)
+ group._v_attrs.CellType = celltypes[prefix]
+
+ print 'Wrote', h5.filename
+ h5.close()
+
+ return 0
+
+if __name__ == '__main__':
+ import sys
+ ret = main()
+ sys.exit(ret)
+
Copied: cs/cigma/trunk/examples/meshes/triangle/square.sh (from rev 14085, cs/cigma/trunk/examples/utils/triangle/square.sh)
===================================================================
--- cs/cigma/trunk/examples/meshes/triangle/square.sh (rev 0)
+++ cs/cigma/trunk/examples/meshes/triangle/square.sh 2009-02-18 16:14:41 UTC (rev 14086)
@@ -0,0 +1,20 @@
+#!/bin/bash
+
+# refine square1
+o="-evn"
+./triangle $o -p square1
+./triangle $o -a0.30 square1.1
+./triangle $o -a0.20 square1.2
+./triangle $o -a0.10 square1.3
+./triangle $o -a0.05 square1.4
+./triangle $o -a0.02 square1.5
+
+# refine square2
+o="-evn -o2"
+./triangle $o -p square2
+./triangle $o -a0.30 square2.1
+./triangle $o -a0.20 square2.2
+./triangle $o -a0.10 square2.3
+./triangle $o -a0.05 square2.4
+./triangle $o -a0.02 square2.5
+
Copied: cs/cigma/trunk/examples/meshes/triangle/square1.poly (from rev 14085, cs/cigma/trunk/examples/utils/triangle/square1.poly)
===================================================================
--- cs/cigma/trunk/examples/meshes/triangle/square1.poly (rev 0)
+++ cs/cigma/trunk/examples/meshes/triangle/square1.poly 2009-02-18 16:14:41 UTC (rev 14086)
@@ -0,0 +1,15 @@
+# A rectangle with four points in 2D, no attributes, no boundary marker
+4 2 0 0
+# List of the rectangle vertices
+1 -1.0 -1.0
+2 +1.0 -1.0
+3 +1.0 +1.0
+4 -1.0 +1.0
+# Four segments, no boundary markers
+4 0
+1 1 2
+2 2 3
+3 3 4
+4 4 1
+# No holes
+0
Copied: cs/cigma/trunk/examples/meshes/triangle/square2.poly (from rev 14085, cs/cigma/trunk/examples/utils/triangle/square2.poly)
===================================================================
--- cs/cigma/trunk/examples/meshes/triangle/square2.poly (rev 0)
+++ cs/cigma/trunk/examples/meshes/triangle/square2.poly 2009-02-18 16:14:41 UTC (rev 14086)
@@ -0,0 +1,15 @@
+# A rectangle with four points in 2D, no attributes, no boundary marker
+4 2 0 0
+# List of the rectangle vertices
+1 -1.0 -1.0
+2 +1.0 -1.0
+3 +1.0 +1.0
+4 -1.0 +1.0
+# Four segments, no boundary markers
+4 0
+1 1 2
+2 2 3
+3 3 4
+4 4 1
+# No holes
+0
Added: cs/cigma/trunk/examples/test/compare-cube.sh
===================================================================
--- cs/cigma/trunk/examples/test/compare-cube.sh (rev 0)
+++ cs/cigma/trunk/examples/test/compare-cube.sh 2009-02-18 16:14:41 UTC (rev 14086)
@@ -0,0 +1,18 @@
+#!/bin/bash
+
+# Comparing test.cube and zero should give (L2)^2 = sqrt(1160/27) = 6.55461386833450037082
+
+echo "Result should be 6.55461386833450037082"
+echo
+
+cigma compare test.cube zero -m cube.h5:/tet4/mesh1 -v
+
+cigma compare test.cube zero -m cube.h5:/tet4/mesh2 -v
+
+cigma compare \
+ -a test.cube \
+ -b zero \
+ -m cube.h5:/tet4/mesh3 \
+ -r integration-rules.h5:/tet4/order3 \
+ -v
+
Property changes on: cs/cigma/trunk/examples/test/compare-cube.sh
___________________________________________________________________
Name: svn:executable
+ *
Added: cs/cigma/trunk/examples/test/compare-square.sh
===================================================================
--- cs/cigma/trunk/examples/test/compare-square.sh (rev 0)
+++ cs/cigma/trunk/examples/test/compare-square.sh 2009-02-18 16:14:41 UTC (rev 14086)
@@ -0,0 +1,18 @@
+#!/bin/bash
+
+# Comparing test.square and zero should give (L2)^2 = 2
+
+echo "Result should be 2"
+echo
+
+cigma compare test.square zero -m square.h5:/tri3/mesh1 -v
+
+cigma compare test.square zero -m square.h5:/tri3/mesh2 -v
+
+cigma compare \
+ -a test.square \
+ -b zero \
+ -m square.h5:/tri3/mesh3 \
+ -r integration-rules.h5:/tri3/order3 \
+ -v
+
Property changes on: cs/cigma/trunk/examples/test/compare-square.sh
___________________________________________________________________
Name: svn:executable
+ *
Added: cs/cigma/trunk/examples/test/exact-integrals.txt
===================================================================
--- cs/cigma/trunk/examples/test/exact-integrals.txt (rev 0)
+++ cs/cigma/trunk/examples/test/exact-integrals.txt 2009-02-18 16:14:41 UTC (rev 14086)
@@ -0,0 +1,26 @@
+
+# First, we define the following Maple functions
+Box2 := (f,a,b) -> int(int(f^2, x=a..b), y=a..b);
+Box3 := (f,a,b) -> int(int(int(f^2, x=a..b), y=a..b), z=a..b);
+SBox3 := (f,L) -> int(int(int(f^2, x=0..L), y=0..L), z=-L..0);
+
+Box3(1+x*y*z, 0, 1) = 139/108
+Box3(1, 0, 1) = 1
+
+Box2(1+x*y, 0, 1) = 29/18
+Box2(1, 0, 1) = 1
+
+
+Box3(1+x*y*z, -1, 1) = 224/27
+Box3(1, -1, 1) = 8
+
+Box2(1+x*y, -1, 1) = 40/9
+Box2(1, -1, 1) = 4
+
+SBox3(1+x*y*z, 24) = 97796961792
+SBox3(1, 24) = 13824
+
+SBox(1+x*y*z, 240000) = 9.7844723711952224256000013824000000000 * 10^37
+SBox(1, 24000) = 1.3824 * 10^13
+
+
Deleted: cs/cigma/trunk/examples/utils/tetgen/README
===================================================================
--- cs/cigma/trunk/examples/utils/tetgen/README 2009-02-18 16:14:40 UTC (rev 14085)
+++ cs/cigma/trunk/examples/utils/tetgen/README 2009-02-18 16:14:41 UTC (rev 14086)
@@ -1,4 +0,0 @@
-
-TetGen: A Quality Tetrahedral Mesh Generator and Three-Dimensional Delaunay Triangulator
-http://tetgen.berlios.de/
-
Deleted: cs/cigma/trunk/examples/utils/tetgen/cube.sh
===================================================================
--- cs/cigma/trunk/examples/utils/tetgen/cube.sh 2009-02-18 16:14:40 UTC (rev 14085)
+++ cs/cigma/trunk/examples/utils/tetgen/cube.sh 2009-02-18 16:14:41 UTC (rev 14086)
@@ -1,21 +0,0 @@
-#!/bin/bash
-
-# refine cube1
-o="-fen"
-./tetgen $o -p cube1
-./tetgen $o -r -a0.500 cube1.1
-./tetgen $o -r -a0.250 cube1.2
-./tetgen $o -r -a0.200 cube1.3
-./tetgen $o -r -a0.125 cube1.4
-./tetgen $o -r -a0.075 cube1.5
-./tetgen $o -r -a0.010 cube1.6
-
-# refine cube2
-o="-fen -o2"
-./tetgen $o -p cube2
-./tetgen $o -r -a0.500 cube2.1
-./tetgen $o -r -a0.250 cube2.2
-./tetgen $o -r -a0.200 cube2.3
-./tetgen $o -r -a0.125 cube2.4
-./tetgen $o -r -a0.075 cube2.5
-./tetgen $o -r -a0.010 cube2.6
Deleted: cs/cigma/trunk/examples/utils/tetgen/cube1.poly
===================================================================
--- cs/cigma/trunk/examples/utils/tetgen/cube1.poly 2009-02-18 16:14:40 UTC (rev 14085)
+++ cs/cigma/trunk/examples/utils/tetgen/cube1.poly 2009-02-18 16:14:41 UTC (rev 14086)
@@ -1,36 +0,0 @@
-# Part 1 - node list
-# node count, 3 dim, no attribute, no boundary marker
-8 3 0 0
-# Node index, node coordinates
-1 0.0 0.0 0.0
-2 2.0 0.0 0.0
-3 2.0 2.0 0.0
-4 0.0 2.0 0.0
-5 0.0 0.0 2.0
-6 2.0 0.0 2.0
-7 2.0 2.0 2.0
-8 0.0 2.0 2.0
-
-# Part 2 - facet list
-# facet count, no boundary marker
-6 0
-# facets
-1 # 1 polygon, no hole, no boundary marker
-4 1 2 3 4 # front
-1
-4 5 6 7 8 # back
-1
-4 1 2 6 5 # bottom
-1
-4 2 3 7 6 # right
-1
-4 3 4 8 7 # top
-1
-4 4 1 5 8 # left
-
-# Part 3 - hole list
-0 # no hole
-
-# Part 4 - region list
-0 # no region
-
Deleted: cs/cigma/trunk/examples/utils/tetgen/cube2.poly
===================================================================
--- cs/cigma/trunk/examples/utils/tetgen/cube2.poly 2009-02-18 16:14:40 UTC (rev 14085)
+++ cs/cigma/trunk/examples/utils/tetgen/cube2.poly 2009-02-18 16:14:41 UTC (rev 14086)
@@ -1,36 +0,0 @@
-# Part 1 - node list
-# node count, 3 dim, no attribute, no boundary marker
-8 3 0 0
-# Node index, node coordinates
-1 0.0 0.0 0.0
-2 2.0 0.0 0.0
-3 2.0 2.0 0.0
-4 0.0 2.0 0.0
-5 0.0 0.0 2.0
-6 2.0 0.0 2.0
-7 2.0 2.0 2.0
-8 0.0 2.0 2.0
-
-# Part 2 - facet list
-# facet count, no boundary marker
-6 0
-# facets
-1 # 1 polygon, no hole, no boundary marker
-4 1 2 3 4 # front
-1
-4 5 6 7 8 # back
-1
-4 1 2 6 5 # bottom
-1
-4 2 3 7 6 # right
-1
-4 3 4 8 7 # top
-1
-4 4 1 5 8 # left
-
-# Part 3 - hole list
-0 # no hole
-
-# Part 4 - region list
-0 # no region
-
Deleted: cs/cigma/trunk/examples/utils/triangle/README
===================================================================
--- cs/cigma/trunk/examples/utils/triangle/README 2009-02-18 16:14:40 UTC (rev 14085)
+++ cs/cigma/trunk/examples/utils/triangle/README 2009-02-18 16:14:41 UTC (rev 14086)
@@ -1,4 +0,0 @@
-
-Triangle: A Two-Dimensional Quality Mesh Generator and Delaunay Triangulator
-http://www.cs.cmu.edu/~quake/triangle.html
-
Deleted: cs/cigma/trunk/examples/utils/triangle/square.py
===================================================================
--- cs/cigma/trunk/examples/utils/triangle/square.py 2009-02-18 16:14:40 UTC (rev 14085)
+++ cs/cigma/trunk/examples/utils/triangle/square.py 2009-02-18 16:14:41 UTC (rev 14086)
@@ -1,69 +0,0 @@
-#!/usr/bin/env python
-
-from __future__ import with_statement
-
-import numpy
-import tables
-
-def read_coords(filename):
- coords = None
- with open(filename, 'r') as fp:
- # read header
- line = fp.readline()
- cols = line.split()
- nno, ndim, natt, nbdr = map(int, cols)
- coords = numpy.zeros((nno, ndim), dtype=float)
- # read nodes
- for n in xrange(nno):
- line = fp.readline()
- cols = line.split()
- i = int(cols[0]) - 1
- for j in xrange(ndim):
- coords[i,j] = float(cols[1+j])
- return coords
-
-def read_connect(filename):
- connect = None
- with open(filename, 'r') as fp:
- # read header
- line = fp.readline()
- cols = line.split()
- nel, ndofs, natt = map(int, cols)
- connect = numpy.zeros((nel, ndofs), dtype=int)
- # read elements
- for e in xrange(nel):
- line = fp.readline()
- cols = line.split()
- i = int(cols[0]) - 1
- for j in xrange(ndofs):
- connect[i,j] = int(cols[1+j])
- connect -= 1
- return connect
-
-def main():
-
- h5 = tables.openFile('square.h5', 'w')
-
- levels = range(1,7);
- celltypes = {'square1': 'tri3', 'square2': 'tri6'}
-
- for prefix in ('square1', 'square2'):
- root = h5.createGroup('/', celltypes[prefix])
- for r in levels:
- coords = read_coords('%s.%d.node' % (prefix, r))
- connect = read_connect('%s.%d.ele' % (prefix, r))
- group = h5.createGroup(root, 'mesh%d' % r)
- h5.createArray(group, 'coordinates', coords)
- h5.createArray(group, 'connectivity', connect)
- group._v_attrs.CellType = celltypes[prefix]
-
- print 'Wrote', h5.filename
- h5.close()
-
- return 0
-
-if __name__ == '__main__':
- import sys
- ret = main()
- sys.exit(ret)
-
Deleted: cs/cigma/trunk/examples/utils/triangle/square.sh
===================================================================
--- cs/cigma/trunk/examples/utils/triangle/square.sh 2009-02-18 16:14:40 UTC (rev 14085)
+++ cs/cigma/trunk/examples/utils/triangle/square.sh 2009-02-18 16:14:41 UTC (rev 14086)
@@ -1,20 +0,0 @@
-#!/bin/bash
-
-# refine square1
-o="-evn"
-./triangle $o -p square1
-./triangle $o -a0.30 square1.1
-./triangle $o -a0.20 square1.2
-./triangle $o -a0.10 square1.3
-./triangle $o -a0.05 square1.4
-./triangle $o -a0.02 square1.5
-
-# refine square2
-o="-evn -o2"
-./triangle $o -p square2
-./triangle $o -a0.30 square2.1
-./triangle $o -a0.20 square2.2
-./triangle $o -a0.10 square2.3
-./triangle $o -a0.05 square2.4
-./triangle $o -a0.02 square2.5
-
Deleted: cs/cigma/trunk/examples/utils/triangle/square1.poly
===================================================================
--- cs/cigma/trunk/examples/utils/triangle/square1.poly 2009-02-18 16:14:40 UTC (rev 14085)
+++ cs/cigma/trunk/examples/utils/triangle/square1.poly 2009-02-18 16:14:41 UTC (rev 14086)
@@ -1,15 +0,0 @@
-# A rectangle with four points in 2D, no attributes, no boundary marker
-4 2 0 0
-# List of the rectangle vertices
-1 -1.0 -1.0
-2 +1.0 -1.0
-3 +1.0 +1.0
-4 -1.0 +1.0
-# Four segments, no boundary markers
-4 0
-1 1 2
-2 2 3
-3 3 4
-4 4 1
-# No holes
-0
Deleted: cs/cigma/trunk/examples/utils/triangle/square2.poly
===================================================================
--- cs/cigma/trunk/examples/utils/triangle/square2.poly 2009-02-18 16:14:40 UTC (rev 14085)
+++ cs/cigma/trunk/examples/utils/triangle/square2.poly 2009-02-18 16:14:41 UTC (rev 14086)
@@ -1,15 +0,0 @@
-# A rectangle with four points in 2D, no attributes, no boundary marker
-4 2 0 0
-# List of the rectangle vertices
-1 -1.0 -1.0
-2 +1.0 -1.0
-3 +1.0 +1.0
-4 -1.0 +1.0
-# Four segments, no boundary markers
-4 0
-1 1 2
-2 2 3
-3 3 4
-4 4 1
-# No holes
-0
Added: cs/cigma/trunk/utils/README
===================================================================
Added: cs/cigma/trunk/utils/create-mesh-from-tetgen-files.py
===================================================================
--- cs/cigma/trunk/utils/create-mesh-from-tetgen-files.py (rev 0)
+++ cs/cigma/trunk/utils/create-mesh-from-tetgen-files.py 2009-02-18 16:14:41 UTC (rev 14086)
@@ -0,0 +1,141 @@
+#!/usr/bin/env python
+
+from __future__ import with_statement
+
+import numpy
+import tables
+
+
+def read_coords(filename):
+ coords = None
+ with open(filename, 'r') as fp:
+ # read header
+ line = fp.readline()
+ cols = line.split()
+ nno, ndim, natt, nbdr = map(int, cols)
+ coords = numpy.zeros((nno, ndim), dtype=float)
+ # read nodes
+ for n in xrange(nno):
+ line = fp.readline()
+ cols = line.split()
+ i = int(cols[0]) - 1
+ for j in xrange(ndim):
+ coords[i,j] = float(cols[1+j])
+ return coords
+
+
+def read_connect(filename):
+ connect = None
+ with open(filename, 'r') as fp:
+ # read header
+ line = fp.readline()
+ cols = line.split()
+ nel, ndofs, natt = map(int, cols)
+ connect = numpy.zeros((nel, ndofs), dtype=int)
+ # read elements
+ for e in xrange(nel):
+ line = fp.readline()
+ cols = line.split()
+ i = int(cols[0]) - 1
+ for j in xrange(ndofs):
+ connect[i,j] = int(cols[1+j])
+ connect -= 1
+ return connect
+
+
+def open_group(h5, location):
+ try:
+ g = h5.getNode(location)
+ except tables.NoSuchNodeError, e:
+ from os.path import dirname, basename
+ parents = dirname(location)
+ group = basename(location)
+ g = h5.createGroup(parents, group, createparents=True)
+ return g
+
+def create_array(h5, location, name, value):
+ try:
+ h5.createArray(location, name, value)
+ except:
+ h5.removeNode(location, name)
+ h5.createArray(location, name, value)
+
+def write_mesh(filename, location, coords, connect, celltype):
+ h5 = tables.openFile(filename, 'a')
+ root = open_group(h5, location)
+ create_array(h5, root, 'coordinates', coords)
+ create_array(h5, root, 'connectivity', connect)
+ root._v_attrs.CellType = celltype
+ h5.close()
+
+
+def main():
+
+ from os.path import exists, splitext
+ from optparse import OptionParser
+
+ parser = OptionParser()
+
+ parser.add_option("-i", "--input",
+ dest="filename", help="Input prefix of tetgen files")
+
+ parser.add_option("-c", "--cell",
+ dest="cell", default='tet4', help="Cell type of mesh (default tet4)")
+
+ parser.add_option("-o", "--output",
+ dest="output", help="Output path for mesh datasets")
+
+ parser.add_option("-v", "--verbose",
+ action="store_true", dest="verbose")
+
+ (options, args) = parser.parse_args()
+
+ if not options.filename:
+ print 'Input file prefix not specified (use -i option)'
+ return 1
+
+ if not options.output:
+ print 'Output path not specified (use -o option)'
+ return 1
+
+ if options.verbose:
+ print "options =", options
+ print "args = ", args
+
+ cell = options.cell
+ prefix = options.filename
+ verbose = options.verbose
+
+ nodfile = '%s.node' % prefix
+ if not exists(nodfile):
+ print "File '%s' does not exist!" % nodfile
+ return 1
+
+ elefile = '%s.ele' % prefix
+ if not exists(elefile):
+ print "File '%s' does not exist!" % elefile
+ return 1
+
+ outfile = '%s.h5' % prefix
+ location = '/'
+ if options.output:
+ o = options.output.split(':')
+ outfile = o[0]
+ if len(o) > 1:
+ location = o[1]
+
+ coords = read_coords(nodfile)
+ connect = read_connect(elefile)
+
+ write_mesh(outfile, location, coords, connect, cell)
+
+ print "Wrote '%s:%s'" % (outfile, location)
+
+ return 0
+
+
+if __name__ == '__main__':
+ import sys
+ ret = main()
+ sys.exit(ret)
+
Property changes on: cs/cigma/trunk/utils/create-mesh-from-tetgen-files.py
___________________________________________________________________
Name: svn:executable
+ *
Added: cs/cigma/trunk/utils/create-tetgen-vol-from-mesh-volume.py
===================================================================
--- cs/cigma/trunk/utils/create-tetgen-vol-from-mesh-volume.py (rev 0)
+++ cs/cigma/trunk/utils/create-tetgen-vol-from-mesh-volume.py 2009-02-18 16:14:41 UTC (rev 14086)
@@ -0,0 +1,86 @@
+#!/usr/bin/env python
+
+import tables
+import numpy
+
+def main():
+
+ from optparse import OptionParser
+
+ parser = OptionParser()
+
+ parser.add_option("-i", "--input",
+ dest="volpath", help="Input HDF5 path to cell volumes")
+
+ parser.add_option("-r", "--residuals",
+ dest="respath", help="Input path to residuals dataset")
+
+ parser.add_option("-o", "--output",
+ dest="output", help="Output file for tetgen .vol file")
+
+ parser.add_option("-v", "--verbose",
+ action="store_true", dest="verbose")
+
+ (options, args) = parser.parse_args()
+
+ if options.verbose:
+ print options
+ print args
+
+ if not options.volpath:
+ print 'No cell volumes specified (use --input option)'
+ return 1
+
+ if not options.output:
+ print 'No output file specified (use --output option)'
+ return 1
+
+ vloc = '/'
+ if options.volpath:
+ o = options.volpath.split(':')
+ volfile = o[0]
+ if len(o) > 1:
+ vloc = o[1]
+
+ h5 = tables.openFile(volfile, 'r')
+ volume = h5.getNode(vloc)[:]
+ h5.close()
+
+ residuals = None
+ if options.respath:
+ o = options.respath.split(':')
+ resfile = o[0]
+ resloc = '/epsilon'
+ if len(o) > 1:
+ resloc = o[1]
+ h5 = tables.openFile(resfile, 'r')
+ residuals = h5.getNode(resloc)[:]
+ h5.close()
+
+ if volume.shape != residuals.shape:
+ print "Shape of residuals %r does not match shape of volume %r" % (residuals.shape, volume.shape)
+ return 1
+
+ #
+ # adjust volume array based on residuals
+ #
+ factor = 1.2 ** 3
+ threshold = 0.25
+ R = threshold * numpy.max(residuals)
+ volume[residuals[:] >= R ] /= factor
+ volume[residuals[:] < R] = -1.0
+
+ outfile = options.output
+ fp = open(outfile, 'w')
+ nel = volume.shape[0]
+ fp.write('%d\n' % nel)
+ for e in xrange(nel):
+ fp.write('%d %g\n' % (e+1, volume[e]))
+
+ print 'Wrote %s' % outfile
+
+ return 0
+
+
+if __name__ == '__main__':
+ main()
Property changes on: cs/cigma/trunk/utils/create-tetgen-vol-from-mesh-volume.py
___________________________________________________________________
Name: svn:executable
+ *
Added: cs/cigma/trunk/utils/power-plot.py
===================================================================
--- cs/cigma/trunk/utils/power-plot.py (rev 0)
+++ cs/cigma/trunk/utils/power-plot.py 2009-02-18 16:14:41 UTC (rev 14086)
@@ -0,0 +1,78 @@
+#!/usr/bin/env python
+
+import numpy
+import matplotlib
+
+from pylab import *
+from scipy.optimize import leastsq
+
+def power_fit(X, Y, **kw):
+
+ assert len(X) == len(Y)
+
+ # parametric function
+ fp = lambda v,x: v[0] * x ** v[1]
+
+ # error function
+ e = lambda v,x,y: (fp(v,x) - y)
+
+ # initial parameter value
+ v0 = (1,1)
+
+ # find optimal parameters
+ (v, success) = leastsq(e, v0, args=(X,Y), maxfev=1000)
+
+ # plot points and best-fit curve
+ print "Best fit: err =", v[0], "* h ^", v[1]
+
+ x = numpy.array(X)
+ y = numpy.array(Y)
+
+ plot_fn = kw.get('plot', loglog)
+ plot_fn(x, y, 'ro')
+ plot_fn(x, fp(v,x))
+ title(r"$\varepsilon\ =\ %f\ h^{%f}$" % (v[0],v[1]))
+
+ show()
+ #savefig('power-fit.png')
+
+ return
+
+
+def read_xy(filename):
+ X = []
+ Y = []
+ fp = open(filename,'r')
+ n = int(fp.readline())
+ for i in xrange(n):
+ line = fp.readline()
+ cols = line.split()
+ x,y = map(float, cols)
+ X.append(x)
+ Y.append(y)
+ fp.close()
+ return X,Y
+
+def main():
+
+ if False:
+ pass
+
+ #X = [ 0.8660254, 0.4330127, 0.21650635, 0.10825318 ]
+ #Y = [ 0.19513078406, 0.0543204762424, 0.0138572906973, 0.00304031950347]
+
+ # laplace square - weighed L2
+ #X = [ 0.088388348, 0.1767767, 0.35355339, 0.70710678 ]
+ #Y = [ 0.0441115235417, 0.202105799913, 0.788022943776, 2.64975185253 ]
+
+ import sys
+ if len(sys.argv) != 2:
+ print 'Usage: power-plot.py XYFILE'
+ sys.exit(1)
+
+ print 'Reading %s' % sys.argv[1]
+ X,Y = read_xy(sys.argv[1])
+ power_fit(X,Y)
+
+if __name__ == '__main__':
+ main()
Property changes on: cs/cigma/trunk/utils/power-plot.py
___________________________________________________________________
Name: svn:executable
+ *
Added: cs/cigma/trunk/utils/remove-hdf5-node.py
===================================================================
--- cs/cigma/trunk/utils/remove-hdf5-node.py (rev 0)
+++ cs/cigma/trunk/utils/remove-hdf5-node.py 2009-02-18 16:14:41 UTC (rev 14086)
@@ -0,0 +1,47 @@
+#!/usr/bin/env python
+
+import tables
+
+def main():
+
+ from optparse import OptionParser
+
+ parser = OptionParser()
+ parser.add_option("-f", "--file", dest="filename")
+ parser.add_option("-n", "--node", dest="node")
+ parser.add_option("-r", "--recursive", action="store_true", dest="recursive")
+ parser.add_option("-v", "--verbose", action="store_true", dest="verbose")
+
+ (options, args) = parser.parse_args()
+
+ if options.verbose:
+ print options
+ print args
+
+ if not options.filename:
+ print "Filename not specified (use --file option)"
+ return 1
+
+ if not options.node:
+ print "No node specified (use --node option)"
+ return 1
+
+ filename = options.filename
+ node = options.node
+ recursive = options.recursive
+
+ h5 = tables.openFile(filename, 'a')
+
+ if options.recursive:
+ h5.removeNode(node, recursive=True)
+ else:
+ h5.removeNode(node)
+
+ print "Removed %s from %s" % (node, filename)
+
+ h5.close()
+
+ return 0
+
+if __name__ == '__main__':
+ main()
Property changes on: cs/cigma/trunk/utils/remove-hdf5-node.py
___________________________________________________________________
Name: svn:executable
+ *
More information about the CIG-COMMITS
mailing list