I am trying to debug some boundary conditions in a complex problem. I have
run into a problem with faceGradAverage or my understanding of it. It is
illustrated in a simple problem
\Del^2 Psi =0 on a unit square
normal derivative = 1 on x=0
normal derivative = 0 on y=0
normal derivative = 0 on y=1
Psi =1 on x=1
The analytical solution to this problem is
Psi = x
The code is in the attached file; skip until "if __name__ == "__main__":" ,
the earlier code converts a Delaunay triangulation to a fipy mesh. I am
doing this to match my more complex case.
Anyway, when the code is run, I get the analytical solution, and so
gradient \Del\Psi should be (1,0) everywhere. It is not. Note also that
since \Del\Psi is spatially uniform, how the gradient is averaged in space
should not make a difference.
How am I confused?
So that you don't have to open the attachment, the code that solves the
problem and prints the face gradients is included as text here:
#setup fipy solution for Del^2 Psi = 1 with BC of 0 everywhere
#first, make grid variable
psi=fipy.CellVariable(name='psi',mesh=mesh,value=0.0)
#apply boundary of Del\Psi \dot \nhat=-1 on x=0 boundary (equivalent to
d Psi/dx=1)
# Del\Psi \dot \nhat=0 on y=0,1 boundary
#and Psi=1 on x=1
xFace,yFace=mesh.faceCenters
faceNorm=mesh.faceNormals
#y=0 boundary
yeq0bound=yFace==0.0
psi.faceGrad.constrain(0.0*faceNorm,where=yeq0bound)
#y=1.0 boundary
yeq1bound=yFace==1.0
psi.faceGrad.constrain(0.0*faceNorm,where=yeq1bound)
#x=0 boundary
xeq0bound=xFace==0.0
psi.faceGrad.constrain(-1.0*faceNorm,where=xeq0bound)
#x=1 boundary
xeq1bound=xFace==1.0
psi.constrain(0.0,where=xeq1bound)
#make equatiion and solve
eq=(fipy.DiffusionTerm(var=psi,coeff=1.0))
eq.solve(var=psi)
#now plot solution
subplot(2,1,2)
tripcolor(psi.mesh.vertexCoords[0,:],psi.mesh.vertexCoords[1,:],
tri.simplices,psi.numericValue)
colorbar()
title('solution')
#now print faceGradients
print 'Face Gradients on x=0, should be
1,0',psi.faceGradAverage[:,xeq0bound]
print ' '
print 'Face Gradients on y=0, should be
1,0',psi.faceGradAverage[:,yeq0bound]
Thanks,
Jamie
# These modules takes the output of scipy.spatial.Delaunay's
# triangulation of data points, and returns both the FiPy mesh for the
# points, and a list of the boundary edges so that boundary conditions
# can be set for the points.
#
# This code handles the case where triangles/simplices have been
# removed from output of the Delaunay triangulation after the
# triangulation. This triangle removal can be useful to create
# interior holes in the domain, or to handle concave boundaries of the
# domain.
#
# The documentation for the routines can be found in there docstrings,
# and examples of there use can be found in __main__
from numpy import *
import fipy
from collections import defaultdict
import copy as Copy
def ScipyDelaunay2mesh(points,simplices,debug=False):
'''ScipyDelaunay2mesh(points,simplices):
Creates a FiPy 2Dmesh object from the output of a
scipy.spatial.Delaunay triangulation,
Also returns a list of border faces for the application of
boundary conditions. The list of border conditions is computed
within this routine; it will handle the case where
triangles/simplices have been removed from the output of the
Delaunay() triangulation.
Parameters
----------
points : array
an (N,2) array of the location of N points that made the
triangulation. If Delaunay is called with
tri=scipy.spatial.Delaunay(), points is tri.points
simplices : array
an (M,3) array of the location of M triangles returned by
Delaunay. If Delaunay is called with
tri=scipy.spatial.Delaunay(), simplices is tri.simplices
debug=False : Boolean
if this is set to true, the mesh will be drawn in blue,
boundaries will be drawn in red, boundary point numbers will be
written at vertices, and numbers for triangles will be drawn in
the triangle. This output is useful for finding where triangles
have been removed incorrectly.
Returns
-------
mesh: a FiPy mesh2D object
borderList: a list of lists
Each list in borderList is a list of the indices of the faces
as indexed in mesh.cellFaceIDs. Each list is for a single
contiguous border of the domain. These will be used to set the
boundary conditions on the mesh.
'''
#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 points in
#content, but reshaped.
print 'Entering ScipyDelaunay2Mesh; Making vertexCoords'
vertexCoords=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(simplices.shape[0]):
for nface in xrange(3):
if nface==0:
face=simplices[ntri,[0,1]]
elif nface==1:
face=simplices[ntri,[1,2]]
else:
face=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=simplices.T*0
for ntri in xrange(simplices.shape[0]):
for nface in xrange(3):
if nface==0:
face=simplices[ntri,[0,1]]
elif nface==1:
face=simplices[ntri,[1,2]]
else:
face=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)
#now make borderList. We cannot assume the neighbors property of
#the Delaunay triangulation is valid because some triangles may
#have been removed.
#first make vertexPathList, a list of all points that are on
#borders
#first, redo neighbors. Create a defaultdict that when called
#with a new key returns an empty list
whichTriangles=defaultdict(list);
#the entry into the dictionary are the edges of the triangle, sorted
#so smaller vertex is first
print ' Making new neighbor data'
for n in range(simplices.shape[0]):
edge1=tuple(sort(simplices[n,[0,1]]))
edge2=tuple(sort(simplices[n,[1,2]]))
edge3=tuple(sort(simplices[n,[2,0]]))
whichTriangles[edge1].append(n)
whichTriangles[edge2].append(n)
whichTriangles[edge3].append(n)
if debug:
triplot(points[:,0],points[:,1],triangles=simplices)
#to debug, plot over all triangles and plot there number in the centers
for n in range(simplices.shape[0]):
xave=0.0
yave=0.0
for ne in range(3):
xave=xave+points[simplices[n,ne],0]/3.0
yave=yave+points[simplices[n,ne],1]/3.0
text(xave,yave,str(n),fontsize=6)
#to debug, iterate over all edges in whichTriangle, and plot number of neighbors
for edge in whichTriangles.keys():
xedge=0.5*(points[edge[0],0]+points[edge[1],0])
yedge=0.5*(points[edge[0],1]+points[edge[1],1])
if len(whichTriangles[edge])==1:
plot(xedge,yedge,'r*',mec='none')
text(points[edge[0],0],points[edge[0],1],str(edge[0]),color='r',fontsize=9)
text(points[edge[1],0],points[edge[1],1],str(edge[1]),color='r',fontsize=9)
elif len(whichTriangles[edge])==2:
plot(xedge,yedge,'g*',mec='none')
else:
assert False, 'wrong number of neighbors'
#now, for each edge/face, see if it has any neighbors
neighbors=simplices*0
for n in range(simplices.shape[0]):
edge1=tuple(sort(simplices[n,[0,1]]))
edge2=tuple(sort(simplices[n,[1,2]]))
edge3=tuple(sort(simplices[n,[2,0]]))
edges=[edge1,edge2,edge3]
for ne in range(3):
edge=edges[ne]
if len(whichTriangles[edge])==1:
neighbors[n,ne]=-1
elif len(whichTriangles[edge])==2:
if whichTriangles[edge][0]==n:
neighbors[n,ne]=whichTriangles[edge][1]
else:
neighbors[n,ne]=whichTriangles[edge][0]
assert whichTriangles[edge][1]==n,'Umm, how can this be true?!?!'
else:
assert False, 'I am confused, simplex has other than 0 or 1 neighbor...'
#identify border edges/face; they are the ones where
#neighbor=-1. Collect edgePoints, a list of tuples that define
#points that connect along the edges of the domain
edgePoints=[]
for n in range(simplices.shape[0]):
edge1=sort(simplices[n,[0,1]])
edge2=sort(simplices[n,[1,2]])
edge3=sort(simplices[n,[2,0]])
if neighbors[n,0]==-1:
edgePoints.append((edge1[0],edge1[1]))
if debug:
plot(points[[edge1[0],edge1[1]],0],
points[[edge1[0],edge1[1]],1],'r-',linewidth=4,alpha=0.2)
if neighbors[n,1]==-1:
edgePoints.append((edge2[0],edge2[1]))
if debug:
plot(points[[edge2[0],edge2[1]],0],
points[[edge2[0],edge2[1]],1],'r-',linewidth=4,alpha=0.2)
if neighbors[n,2]==-1:
edgePoints.append((edge3[0],edge3[1]))
if debug:
plot(points[[edge3[0],edge3[1]],0],
points[[edge3[0],edge3[1]],1],'r-',linewidth=4,alpha=0.2)
#ok, now make vertexPathList. Make a dictionary of paths so
#paths(npoint) returns list of points the path is connected to.
print ' Now make vertexPathList'
paths=defaultdict(list)
for np in xrange(len(edgePoints)):
edge=edgePoints[np]
paths[edge[0]].append(edge[1])
paths[edge[1]].append(edge[0])
#now follow each path to find a cycle; store each cycle in vertexPathList,
#which is a list of lists. Be sure to remove each path after taking it.
vertexPathList=[]
keys=paths.keys()
for key in keys:
edges=paths[key]
if len(edges)==0:
assert True,'nothing to do here'
else:
thispath=[key] #start of path
nhere=paths[key].pop()
nfrom=key+0
cols='rgbcmyk'
while (nhere != key):
if False:
plot([points[nhere,0],points[nfrom,0]]
,[points[nhere,1],points[nfrom,1]],cols[remainder(key,7)]+'*-',
mec='none',markersize=12,linewidth=2)
#remove path whence we came
paths[nhere].remove(nfrom)
#add this point to path
thispath.append(nhere)
#now move to next point, also remove that path:
nfrom=nhere+0
nhere=paths[nhere].pop()
vertexPathList.append(thispath)
#vertexPathList contains a list of lists of vertex points that
#surround all islands and (at item nOcean) the model domain
#boundary. All of paths are closed. I need to convert this into a
#list of faces, as defined in faceMap and faceVertexIDs. I also
#need to make sure to include the face that connects the first and
#last element of the list.
#remember faceMap is a dictionary with a key made up of a sorted tuple of vertex
#indices which maps from a face to its location in faceVertexIDs
#first, make a copy of vertexPathList, since all lists will be the same size
borderList=Copy.deepcopy(vertexPathList)
#now fill fipy BoundaryList with good data
for ncycle in xrange(len(vertexPathList)):
for npoint in xrange(len(vertexPathList[ncycle])):
nOtherPoint=remainder(npoint+1,len(vertexPathList[ncycle]))
faceKey=tuple(sort((vertexPathList[ncycle][npoint],
vertexPathList[ncycle][nOtherPoint])))
borderList[ncycle][npoint]=faceMap[faceKey]
return (mesh,borderList)
if __name__ == "__main__":
from pylab import *
from scipy.spatial import Delaunay
print ' '
print 'Test case -- square domain with Del^2 Psi=0 and various BC'
#set up figure
figure(1)
clf()
#make points in square domain of size 1
points=[]
xvec=linspace(0.0,1.0,30)
xmat,ymat=meshgrid(xvec,xvec)
points=zip(xmat.flatten(),ymat.flatten())
#to debug shuffle points
shuffle(points)
#now make triangulation, and convert to mesh
tri=Delaunay(points)
mesh,borderList=ScipyDelaunay2mesh(tri.points,tri.simplices)
#now plot mesh
figure(1)
clf()
subplot(2,1,1)
triplot(mesh.vertexCoords[0,:],mesh.vertexCoords[1,:],triangles=tri.simplices)
title('simple square mesh')
draw()
show()
#setup fipy solution for Del^2 Psi = 1 with BC of 0 everywhere
#first, make grid variable
psi=fipy.CellVariable(name='psi',mesh=mesh,value=0.0)
#apply boundary of Del\Psi \dot \nhat=1 on x=0 boundary
# Del\Psi \dot \nhat=0 on y=0,1 boundary
#and Psi=0 on x=1
xFace,yFace=mesh.faceCenters
faceNorm=mesh.faceNormals
#y=0 boundary
yeq0bound=yFace==0.0
psi.faceGrad.constrain(0.0*faceNorm,where=yeq0bound)
#y=1.0 boundary
yeq1bound=yFace==1.0
psi.faceGrad.constrain(0.0*faceNorm,where=yeq1bound)
#x=0 boundary
xeq0bound=xFace==0.0
psi.faceGrad.constrain(-1.0*faceNorm,where=xeq0bound)
#x=1 boundary
xeq1bound=xFace==1.0
psi.constrain(0.0,where=xeq1bound)
#make equatiion and solve
eq=(fipy.DiffusionTerm(var=psi,coeff=1.0))
eq.solve(var=psi)
#now plot solution
subplot(2,1,2)
tripcolor(psi.mesh.vertexCoords[0,:],psi.mesh.vertexCoords[1,:],
tri.simplices,psi.numericValue)
colorbar()
title('solution')
#now print faceGradients
print 'Face Gradients on x=0, should be 1,0',psi.faceGradAverage[:,xeq0bound]
print ' '
print 'Face Gradients on y=0, should be 1,0',psi.faceGradAverage[:,yeq0bound]
_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
