For the use of others --
Sometimes when dealing with arbitrary grids, it is useful to use the
output of routines other than Gmsh to create the fipy mesh. In my case, I
need to create a complex grid from ocean bathymetry, and I ended up using
scipy.spatial.Delaunay to do so. With help from Daniel and Jonathan, I
created a code to create a mesh2D object from the output of Delaunay. It is
written to emphasize clarity over speed. I hope others find it useful.
James Pringle
#==================================================
#This code creates a fipy grid from the output of a Delaunay
#triangulation. The code is written to prioratize clarity over speed.
from pylab import *
from numpy import *
import fipy
from scipy.spatial import Delaunay
#make a simple example set of points to triangulate with delaunay
points=array([[ 0., 0.],
[ 0., 1.],
[ 1., 0.],
[ 1., 1.]])
tri=Delaunay(points)
#when this code is run as is, it should produce
# faceVertexIDs=
# [[3, 1, 0, 2, 0
# [1, 0, 3, 3, 2]]
#and
# cellFaceIds=
# [[0, 3],
# [1, 2],
# [2, 4]]
#now create the arrays that fipy.mesh.mesh2D needs. vertexCoords is
#of shape (2,N), where N is the number of points, and contains the
#coordinates of the vertices. It is equivalent to tri.points in
#content, but reshaped.
print 'Making vertexCoords'
vertexCoords=tri.points.T
#faceVertexID is of shape (2,M) where M is the number of faces/edges
#(faces is the fipy terminology, edges is the scipy.spatial.Delaunay
#terminology). faceVertexID contains the points (as indices into
#vertexCoords) which make up each face. This data can be extracted
#from simplices from Delaunay, which contains the faces of each
#triangle. HOWEVER, simplices contains many repeated faces, since each
#triangle will border on others. Only one copy of each face should be
#inserted into faceVertexID.
print 'Making faceVertexIDs'
faceIn={} #this is a dictionary of all faces that have already been
#inserted into the faceVertexID, with the key a tuple of
#indices, sorted
faceVertexIDs_list=[] #structure to build
for ntri in xrange(tri.simplices.shape[0]):
for nface in xrange(3):
if nface==0:
face=tri.simplices[ntri,[0,1]]
elif nface==1:
face=tri.simplices[ntri,[1,2]]
else:
face=tri.simplices[ntri,[2,0]]
faceKey=tuple(sort(face))
if not (faceKey in faceIn):
faceIn[faceKey]=True
faceVertexIDs_list.append(face)
faceVertexIDs=array(faceVertexIDs_list).T
#now create cellFaceIDs of shape (3,P) where P is the number of cells
#in the domain. It contains the faces (as indices into faceVertexIDs)
#that make up each cell.
#first create dictionary with a key made up of a sorted tuple of vertex
#indices which maps from a face to its location in faceVertexIDs
print 'Making cellFaceIDs'
faceMap={}
for nface in xrange(faceVertexIDs.shape[1]):
faceKey=tuple(sort(faceVertexIDs[:,nface]))
faceMap[faceKey]=nface
#now loop over simplices, find each edge/face, and map from points to
#faces
cellFaceIDs=tri.simplices.T*0
for ntri in xrange(tri.simplices.shape[0]):
for nface in xrange(3):
if nface==0:
face=tri.simplices[ntri,[0,1]]
elif nface==1:
face=tri.simplices[ntri,[1,2]]
else:
face=tri.simplices[ntri,[2,0]]
faceKey=tuple(sort(face))
cellFaceIDs[nface,ntri]=faceMap[faceKey]
#ok, now instantiate a mesh2D object with this data
print 'Making mesh'
mesh=fipy.meshes.mesh2D.Mesh2D(vertexCoords,faceVertexIDs,cellFaceIDs)
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
