[cig-commits] r4448 - in cs: . benchmark benchmark/trunk
benchmark/trunk/examples
luis at geodynamics.org
luis at geodynamics.org
Tue Aug 29 12:09:57 PDT 2006
Author: luis
Date: 2006-08-29 12:09:57 -0700 (Tue, 29 Aug 2006)
New Revision: 4448
Added:
cs/benchmark/
cs/benchmark/branches/
cs/benchmark/tags/
cs/benchmark/trunk/
cs/benchmark/trunk/bm.py
cs/benchmark/trunk/examples/
cs/benchmark/trunk/examples/eric-import.py
cs/benchmark/trunk/examples/greg-import.py
Log:
New import
Added: cs/benchmark/trunk/bm.py
===================================================================
--- cs/benchmark/trunk/bm.py 2006-08-29 16:37:40 UTC (rev 4447)
+++ cs/benchmark/trunk/bm.py 2006-08-29 19:09:57 UTC (rev 4448)
@@ -0,0 +1,388 @@
+#!/usr/bin/env python
+
+import numpy
+import tables
+
+from numpy import array, zeros
+
+class TetrahedralMesh(object):
+
+
+ def __init__(self, **kwargs):
+ self.nno = kwargs.get('nodes') # number of nodes
+ self.nel = kwargs.get('elements') # number of elements
+ self.deg = kwargs.get('degree', 1) # degree of element
+ self.qd = kwargs.get('quadrature', 1) # degree of quadrature rule
+ self.nqd = self.qd ** 3 # number of quadrature pts
+
+
+ def open(self, filename, mode='r'):
+ self.filename = filename
+ if mode in ('r', 'a'):
+ self.h5 = tables.openFile(filename, mode)
+ self.nno = self.h5.root._v_attrs.nno
+ self.nel = self.h5.root._v_attrs.nel
+ self.qd = self.h5.root._v_attrs.qd
+ self.nqd = self.qd ** 3
+ elif mode == 'w':
+ self.h5 = tables.openFile(filename, 'w')
+ self.h5.createGroup("/", "solution")
+ self.h5.createGroup("/", "geometry")
+ self.h5.createGroup("/", "topology")
+ self.h5.createGroup("/", "tmp")
+ self.h5.createGroup("/geometry", "map")
+ self.h5.createGroup("/geometry", "quadrature")
+ self.h5.root._v_attrs.nno = self.nno
+ self.h5.root._v_attrs.nel = self.nel
+ self.h5.root._v_attrs.qd = self.qd
+ return
+
+ def load_connect(self, file, skip=0):
+ """
+ Load node connectivity from file.
+ """
+ for n in xrange(skip):
+ line = file.readline()
+ print "Loading /topology/connect"
+ connect = zeros((self.nel, 4), dtype=int)
+ for e in xrange(self.nel):
+ line = file.readline()
+ e,etp,mat,inf,n1,n2,n3,n4 = map(int, line.split())
+ connect[e-1] = (n1, n2, n3, n4)
+ self.h5.createArray("/topology", "connect", connect)
+ file.close()
+
+ def load_coord(self, file, skip=0):
+ """
+ Load node coordinates from file.
+ """
+ for n in xrange(skip):
+ line = file.readline()
+ print "Loading /geometry/coord"
+ coord = zeros((self.nno, 3), dtype=float)
+ for n in xrange(self.nno):
+ line = file.readline()
+ cols = line.split()
+ x,y,z = map(float, (cols[1], cols[2], cols[3]))
+ coord[n] = (1000*x,1000*y,1000*z)
+ self.h5.createArray("/geometry", "coord", coord)
+ file.close()
+
+ def load_solution(self, file, skip=0, **kw):
+ """
+ Load solution from file.
+ """
+ for n in xrange(skip):
+ line = file.readline()
+ print "Loading /solution/{displacement,velocity}"
+ nc = kw.get('ncol', 0)
+ xc = kw.get('xcol', 1)
+ vxc = kw.get('vxcol', 4)
+ yc,zc = xc+1, xc+2
+ vyc, vzc = vxc+1, vxc+2
+ displacement = zeros((self.nno, 3), dtype=float)
+ velocity = zeros((self.nno, 3), dtype=float)
+ for n in xrange(self.nno):
+ line = file.readline()
+ cols = line.split()
+ n = int(cols[nc]) - 1
+ x,y,z = map(float, (cols[xc], cols[yc], cols[zc]))
+ vx,vy,vz = map(float, (cols[vxc], cols[vyc], cols[vzc]))
+ displacement[n] = (x,y,z)
+ velocity[n] = (vx,vy,vz)
+ self.h5.createArray("/solution", "displacement", displacement)
+ self.h5.createArray("/solution", "velocity", velocity)
+
+ def load_map(self):
+ print "Loading /geometry/map/{jacobian,to_ref,from_ref,origin}"
+ from numpy.linalg import inv, det
+ coord = self.h5.root.geometry.coord[:]
+ connect = self.h5.root.topology.connect[:]
+ jacobian = zeros((self.nel,), dtype=float)
+ map_to_ref = zeros((self.nel, 3, 3), dtype=float)
+ map_from_ref = zeros((self.nel, 3, 3), dtype=float)
+ map_origin = zeros((self.nel, 3), dtype=float)
+ for e in xrange(self.nel):
+ # print "processing element", e
+ a,b,c,d = connect[e]
+ x0,y0,z0 = coord[a-1] # a = (x0,y0,z0)
+ x1,y1,z1 = coord[b-1] # b = (x1,y1,z1)
+ x2,y2,z2 = coord[c-1] # c = (x2,y2,z2)
+ x3,y3,z3 = coord[d-1] # d = (x3,y3,z3)
+ j = array([[(x1 - x0)/2, (x2 - x0)/2, (x3 - x0)/2],
+ [(y1 - y0)/2, (y2 - y0)/2, (y3 - y0)/2],
+ [(z1 - z0)/2, (z2 - z0)/2, (z3 - z0)/2]], dtype=float)
+ jacobian[e] = det(j)
+ map_from_ref[e] = j
+ map_to_ref[e] = inv(j)
+ map_origin[e] = ((x1+x2+x3-x0)/2, (y1+y2+y3-y0)/2, (z1+z2+z3-z0)/2)
+ self.h5.createArray("/geometry/map", "jacobian", jacobian)
+ self.h5.createArray("/geometry/map", "to_ref", map_to_ref)
+ self.h5.createArray("/geometry/map", "from_ref", map_from_ref)
+ self.h5.createArray("/geometry/map", "origin", map_origin)
+
+ def load_normals(self):
+ print "Loading /geometry/normals"
+ from numpy import cross
+ connect = self.h5.root.topology.connect[:]
+ coord = self.h5.root.geometry.coord[:]
+ normals = zeros((self.nel, 4, 3), dtype=float)
+ for e in xrange(self.nel):
+ a,b,c,d = connect[e]
+ ra = coord[a-1]
+ rb = coord[b-1]
+ rc = coord[c-1]
+ rd = coord[d-1]
+ normals[e,0] = cross(rc-ra, rb-ra)
+ normals[e,1] = cross(ra-rd, rb-rd)
+ normals[e,2] = cross(rb-rd, rc-rd)
+ normals[e,3] = cross(rc-rd, ra-rd)
+ self.h5.createArray("/geometry", "normals", normals)
+
+ def load_quadrature(self):
+ print "Loading /geometry/quadrature/{points,weights}"
+ from numpy import dot
+ from FIAT.shapes import TETRAHEDRON
+ from FIAT.quadrature import make_quadrature
+ Q = make_quadrature(TETRAHEDRON, self.qd)
+ x = array(Q.x)
+ w = array(Q.w)
+ self.h5.createArray("/geometry/quadrature", "points", x)
+ self.h5.createArray("/geometry/quadrature", "weights", w)
+
+ def load_quadrature2(self, file, skip=0):
+ """
+ Instead of using FIAT, load quadrature points and weights
+ from an external file.
+ """
+ for n in xrange(skip):
+ line = file.readline()
+ print "Loading /geometry/quadrature/{points,weights}"
+ nel, nqd = self.nel, self.nqd
+ x = zeros((nqd, 3), dtype=float)
+ w = zeros((nqd,), dtype=float)
+ for q in xrange(nqd):
+ line = file.readline()
+ wq,xq,yq,zq = map(float, line.split())
+ x[q] = (xq,yq,zq)
+ w[q] = wq
+ file.close()
+ self.h5.createArray("/geometry/quadrature", "points", x)
+ self.h5.createArray("/geometry/quadrature", "weights", w)
+
+ def load_qpoints(self):
+ print "Loading /tmp/qpoints"
+ from numpy import dot
+ nel, nqd = self.nel, self.nqd
+ map_origin = self.h5.root.geometry.map.origin[:]
+ map_from_ref = self.h5.root.geometry.map.from_ref[:]
+ qpoints = zeros((nel, nqd, 3), dtype=float)
+ x = self.h5.root.geometry.quadrature.points[:]
+ for e in xrange(nel):
+ j = map_from_ref[e]
+ origin = map_origin[e]
+ for q in xrange(nqd): # TODO: vectorize this loop
+ qpoints[e,q] = origin + dot(j,x[q])
+ self.h5.createArray("/tmp", "qpoints", qpoints)
+
+ def write_qpoints(self, file):
+ qpoints = self.h5.root.tmp.qpoints[:]
+ nel = qpoints.shape[0]
+ nqd = qpoints.shape[1]
+ for e in xrange(nel):
+ for q in xrange(nqd):
+ point = qpoints[e,q]
+ print >> file, point[0], point[1], point[2]
+ file.close()
+
+ def load_projection(self):
+ """
+ Project solutions to own /tmp/qpoints, using FIAT to evaluate
+ the shape functions at the specified locations.
+ """
+ from numpy import dot
+ from FIAT.Lagrange import Lagrange
+ from FIAT.shapes import TETRAHEDRON
+
+ nel, nqd = self.nel, self.nqd
+
+ U = Lagrange(TETRAHEDRON, self.deg)
+ Ufs = U.function_space()
+
+ points = self.h5.root.geometry.quadrature.points[:]
+ qx = [map(float, tuple(points[q])) for q in xrange(nqd)]
+ N = array(Ufs.tabulate(qx)).transpose()
+
+ connect = self.h5.root.topology.connect[:]
+ x = self.h5.root.solution.displacement[:]
+ vx = self.h5.root.solution.velocity[:]
+
+ displacement = zeros((nel, nqd, 3), dtype=float)
+ velocity = zeros((nel, nqd, 3), dtype=float)
+
+ vert = zeros((4,3), dtype=float)
+
+ print "Loading /tmp/displacement"
+ for e in xrange(nel):
+ a,b,c,d = connect[e]
+ vert[:] = (x[a-1], x[b-1], x[c-1], x[d-1])
+ displacement[e] = dot(N, vert)
+ self.h5.createArray("/tmp", "displacement", displacement)
+
+ print "Loading /tmp/velocity"
+ for e in xrange(nel):
+ a,b,c,d = connect[e]
+ vert[:] = (vx[a-1], vx[b-1], vx[c-1], vx[d-1])
+ velocity[e] = dot(N, vert)
+ self.h5.createArray("/tmp", "velocity", velocity)
+
+ def load_projection2(self, file):
+ """
+ Load an external projection (from a file).
+ """
+ print "Loading /tmp/{displacement2,velocity2}"
+ nel, nqd = self.nel, self.nqd
+ displacement = zeros((nel, nqd, 3), dtype=float)
+ velocity = zeros((nel, nqd, 3), dtype=float)
+ for e in xrange(nel):
+ for q in xrange(nqd):
+ line = file.readline()
+ x,y,z,vx,vy,vz = map(float, line.split())
+ displacement[e,q] = (x,y,z)
+ velocity[e,q] = (vx,vy,vz)
+ file.close()
+ self.h5.createArray("/tmp", "displacement2", displacement)
+ self.h5.createArray("/tmp", "velocity2", velocity)
+
+ def load_displacement2(self, file):
+ """
+ Load an external projection (from a file).
+ """
+ print "Loading /tmp/displacement2"
+ nel, nqd = self.nel, self.nqd
+ displacement2 = zeros((nel, nqd, 3), dtype=float)
+ for e in xrange(nel):
+ for q in xrange(nqd):
+ line = file.readline()
+ x,y,z,dx,dy,dz = map(float, line.split())
+ displacement2[e,q] = (dx,dy,dz)
+ file.close()
+ self.h5.createArray("/tmp", "displacement2", displacement2)
+
+
+ def project(A, B, field=None):
+ """
+ Project own mesh A into destination mesh B. That is, evaluate
+ the specified field in A by using B's quadrature points. Assumes
+ B has been opened in append mode. The result is stored on B
+ as /tmp/field2.
+ """
+ from numpy import dot
+ from FIAT.Lagrange import Lagrange
+ from FIAT.shapes import TETRAHEDRON
+
+ U = Lagrange(TETRAHEDRON, A.deg)
+ Ufs = U.function_space()
+
+ connect = A.h5.root.topology.connect[:]
+ coord = A.h5.root.geometry.coord[:]
+ normals = A.h5.root.geometry.normals[:]
+ map_to_ref = A.h5.root.geometry.map.to_ref[:]
+ origin = A.h5.root.geometry.map.origin[:]
+ qpoints = B.h5.root.tmp.qpoints[:]
+
+ f = getattr(A.h5.root.solution, field)[:]
+ f2 = zeros(qpoints.shape, dtype=float)
+
+ vert = zeros((4,3), dtype=float)
+
+ def contained(e,x):
+ a = connect[e,0]
+ d = connect[e,3]
+ xa = x - coord[a-1]
+ xd = x - coord[d-1]
+ if dot(normals[e,0], xa) > 0:
+ return False
+ if dot(normals[e,1], xd) > 0:
+ return False
+ if dot(normals[e,2], xd) > 0:
+ return False
+ if dot(normals[e,3], xd) > 0:
+ return False
+ return True
+
+ def brute_force_find(x):
+ for e in xrange(A.nel):
+ if contained(e,x): return e
+ return None
+
+ find_element = brute_force_find
+
+ print "Loading /tmp/%s2 on %s" % (field, B.filename)
+ for eb in xrange(B.nel):
+ for qb in xrange(B.nqd): # iterate over B's quadrature pts
+
+ # coordinates of quadrature point
+ qx = qpoints[eb,qb]
+
+ # find corresponding element in A
+ ea = find_element(qx)
+
+ # map point to reference domain
+ xi = dot(map_to_ref[ea], qx - origin[ea])
+
+ # evaluate A's shape functions at xi
+ N = Ufs.tabulate([tuple(xi)]).transpose()
+ a,b,c,d = connect[ea]
+ vert[:] = (f[a-1], f[b-1], f[c-1], f[d-1])
+ f2[eb,qb] = dot(N,vert)
+
+ B.h5.createArray("/tmp", field+'2', f2)
+
+ def project_displacement(A, B):
+ A.project(B, field='displacement')
+
+ def project_velocity(A,B):
+ A.project(B, field='velocity')
+
+
+
+ def error(self, field=None):
+ """
+ Calculates L2 error between /tmp/field and /tmp/field2
+ """
+ from numpy import dot, sqrt
+ nno, nel, nqd = self.nno, self.nel, self.nqd
+ x1 = getattr(self.h5.root.tmp, field)[:]
+ x2 = getattr(self.h5.root.tmp, field+'2')[:]
+ w = self.h5.root.geometry.quadrature.weights[:].reshape(1, nqd)
+ jacobian = self.h5.root.geometry.map.jacobian[:]
+ global_err = 0.0
+ for e in xrange(nel):
+ dx = x2[e] - x1[e]
+ global_err += jacobian[e] * sum(sum(dot(w,dx**2)))
+ return sqrt(global_err)
+
+ def error_displacement(self):
+ """
+ Calculates L2 error between /tmp/displacement and /tmp/displacement2
+ """
+ return self.error(field='displacement')
+
+ def error_velocity(self):
+ """
+ Calculates L2 error between /tmp/velocity and /tmp/velocity
+ """
+ return self.error(field='velocity')
+
+
+
+ def close(self):
+ self.h5.close()
+
+
+def main():
+ pass
+
+if __name__ == '__main__':
+ main()
Added: cs/benchmark/trunk/examples/eric-import.py
===================================================================
--- cs/benchmark/trunk/examples/eric-import.py 2006-08-29 16:37:40 UTC (rev 4447)
+++ cs/benchmark/trunk/examples/eric-import.py 2006-08-29 19:09:57 UTC (rev 4448)
@@ -0,0 +1,38 @@
+#!/usr/bin/env python
+
+from bm import TetrahedralMesh
+
+
+mesh = TetrahedralMesh(nodes=14874, elements=75653, quadrature=2)
+
+# Open h5 file
+mesh.open('eric-hetland.h5', 'w')
+
+# Load coord file
+fn = 'Old_BCs/bmrsnog_1000m_1.0.coord'
+f = open(fn, 'r')
+mesh.load_coord(f, skip=6)
+
+# Load connect file
+fn = 'Old_BCs/bmrsnog_1000m_1.0.connect'
+f = open(fn, 'r')
+mesh.load_connect(f, skip=3)
+
+# Load displacement and velocities file
+fn = 'Old_BCs/bmrsnog_1000m_1.0.mesh.time.00000.inp'
+f = open(fn, 'r')
+mesh.load_solution(f, skip=7)
+
+# Calculate some things (order is important)
+mesh.load_map()
+mesh.load_normals()
+mesh.load_quadrature()
+mesh.load_qpoints()
+mesh.load_projection()
+
+# Write down list of quadrature points
+f = open('eric-qpoints.in', 'w')
+mesh.write_qpoints(f)
+
+# done
+mesh.close()
Property changes on: cs/benchmark/trunk/examples/eric-import.py
___________________________________________________________________
Name: svn:executable
+
Added: cs/benchmark/trunk/examples/greg-import.py
===================================================================
--- cs/benchmark/trunk/examples/greg-import.py 2006-08-29 16:37:40 UTC (rev 4447)
+++ cs/benchmark/trunk/examples/greg-import.py 2006-08-29 19:09:57 UTC (rev 4448)
@@ -0,0 +1,38 @@
+#!/usr/bin/env python
+
+import bm
+from bm import TetrahedralMesh
+
+mesh = TetrahedralMesh(nodes=14874, elements=75653, quadrature=2)
+
+# Open h5 file
+mesh.open('greg-lyzenga.h5', 'w')
+
+# Load coord file
+fn = 'bmrsnog_1000m.coord'
+f = open(fn, 'r')
+mesh.load_coord(f, skip=1)
+
+# Load connect file
+fn = 'bmrsnog_1000m.connect'
+f = open(fn, 'r')
+mesh.load_connect(f, skip=0)
+
+# Load displacement and velocities file
+fn = 'bmrsnog_1kmA.out'
+f = open(fn, 'r')
+mesh.load_solution(f, skip=9, ncol=1, xcol=5, vxcol=8)
+
+# Calculate some things
+mesh.load_map()
+mesh.load_normals()
+mesh.load_quadrature()
+mesh.load_qpoints()
+mesh.load_projection()
+
+# Write down list of quadrature points
+f = open('greg-qpoints.in', 'w')
+mesh.write_qpoints(f)
+
+# done
+mesh.close()
Property changes on: cs/benchmark/trunk/examples/greg-import.py
___________________________________________________________________
Name: svn:executable
+
More information about the cig-commits
mailing list