James, this is awesome, thanks for posting. It would be a very good
idea to have a DelaunayMesh in FiPy that used Scipy for the
triangulation as it would give another way to make and test
unstructured meshes without the Gmsh dependency. Something like,
mesh = DelaunayMesh(points_in_2D)
In the long run, it would be good to change your code so that it
doesn'y loop in Python. I suspect that your code is quite inefficient
due to these loops. However, you have something that works, which is
great and it is only a one time cost at the beginning of the
simulation. Is efficiency an issue for you with this?
On Thu, Jun 16, 2016 at 4:29 PM, James Pringle <jprin...@unh.edu> wrote:
> 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
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
--
Daniel Wheeler
_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
