Re: scipy's Delaunay output to fipy mesh2D object
Dear all --
A warning about using the above code. Since the choice of triangles
formed by a Delaunay triangulation can be somewhat arbitrary in the case of
a somewhat regular distribution of points, it will choose some
triangulations that lead to severe orthogonality errors in fipy, even if
that was not necessary for the points it was given.
So use with caution and monitor non-orthogonality with the
mesh._nonOrthogonal property. For my mostly regular grid I am resorting to
triangulation-by-hand.
Cheers,
Jamie Pringle
University of New Hampshire
On Fri, Jun 17, 2016 at 10:14 AM, James Pringle wrote:
> 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 <
> daniel.wheel...@gmail.com> 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 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'
Re: scipy's Delaunay output to fipy mesh2D object
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
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 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(ve
Re: scipy's Delaunay output to fipy mesh2D object
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 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 ]
scipy's Delaunay output to fipy mesh2D object
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 ]
