Thanks -- Once I have pounded on this a bit more, I will write the function. It should be straightforward. I concentrated on making it simple and easy to understand; since it is run only to make the grid, efficiency was not a great issue for me.
Jamie On Fri, Jun 17, 2016 at 10:03 AM, Daniel Wheeler <[email protected]> wrote: > 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 <[email protected]> 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 > > [email protected] > > > https://urldefense.proofpoint.com/v2/url?u=http-3A__www.ctcms.nist.gov_fipy&d=CwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=oMf1-WHpUHeD7kN3S7612CDdF2TDPqDF9R-n-71Ks1Y&m=Dq-AXpFxiZiFH2Xhbvb8v4KBbo7s5Yvv49ugj35wl4Q&s=9DekPXt2ye1Peb1uvJttFUQWcZmdREgN8LWfezsaUR4&e= > > [ NIST internal ONLY: > https://urldefense.proofpoint.com/v2/url?u=https-3A__email.nist.gov_mailman_listinfo_fipy&d=CwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=oMf1-WHpUHeD7kN3S7612CDdF2TDPqDF9R-n-71Ks1Y&m=Dq-AXpFxiZiFH2Xhbvb8v4KBbo7s5Yvv49ugj35wl4Q&s=NtcGSFkPViP2zPUUsmgcOmtCf0ljwqmrVVm8zSJnKnw&e= > ] > > > > > > -- > Daniel Wheeler > _______________________________________________ > fipy mailing list > [email protected] > > https://urldefense.proofpoint.com/v2/url?u=http-3A__www.ctcms.nist.gov_fipy&d=CwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=oMf1-WHpUHeD7kN3S7612CDdF2TDPqDF9R-n-71Ks1Y&m=Dq-AXpFxiZiFH2Xhbvb8v4KBbo7s5Yvv49ugj35wl4Q&s=9DekPXt2ye1Peb1uvJttFUQWcZmdREgN8LWfezsaUR4&e= > [ NIST internal ONLY: > https://urldefense.proofpoint.com/v2/url?u=https-3A__email.nist.gov_mailman_listinfo_fipy&d=CwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=oMf1-WHpUHeD7kN3S7612CDdF2TDPqDF9R-n-71Ks1Y&m=Dq-AXpFxiZiFH2Xhbvb8v4KBbo7s5Yvv49ugj35wl4Q&s=NtcGSFkPViP2zPUUsmgcOmtCf0ljwqmrVVm8zSJnKnw&e= > ] >
_______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
