Trouble layering 2D plots using FiPy viewer functions
Hi FiPy community,
I'm solving a system of 2-dimensional PDEs in time, trying to plot the
changing contour lines of individual solutions together with a map of their
sum values.
In the attached example, I defined four uncoupled diffusion equations with
separate centers. I solved them simultaneously, plotted each solution using
the contour function in matplotlib, and plotted their sum using a FiPy
viewer function.
However, I noticed that the plots are layered only when my FiPy viewer
calls *Matplotlib2DGridContourViewer *(line 55), but not when I try to use
*Matplotlib2DGridViewer *or *Matplotlib2DViewer*, both of which failed to
render any image when the other contour functions are called.
My problem requires me to generate a heat map of the solutions' sum, so I'm
wondering if there's any way to work around this pesky issue.
Thanks,
Yun
--
Yun Tao
Postdoc
Center for Infectious Disease Dynamics
Pennsylvania State University
State College, PA 16803
#!/usr/bin/env python
# encoding: utf-8
"""
Created by Yun Tao on 2012-04-01.
Copyright (c) 2012 __MyCompanyName__. All rights reserved.
"""
from fipy import *
from fipy import numerix
from matplotlib.mlab import bivariate_normal
import matplotlib.pyplot as plt
import numpy as np
fig = plt.figure()
txt = fig.suptitle('time zero', fontsize=14)
plt.show()
ax1 = plt.subplot()
# mesh config
l = 6.
nx = 201.
ny = nx
dx = l / nx
dy = dx
mesh = Grid2D(nx=nx, ny=ny, dx=dx, dy=dy) + [[-l/2.]]
# initial condition config
x, y = mesh.getCellCenters()
X = x.reshape(201, -1)
Y = y.reshape(201, -1)
z1 = bivariate_normal(x, y, .2, .2, .3, .3)
z2 = bivariate_normal(x, y, .2, .2, -.3, -.3)
z3 = bivariate_normal(x, y, .2, .2, -.3, .3)
z4 = bivariate_normal(x, y, .2, .2, .3, -.3)
# variable config
phi1 = CellVariable(name='$\phi1$', mesh=mesh, value=z1)
phi2 = CellVariable(name='$\phi2$', mesh=mesh, value=z2)
phi3 = CellVariable(name='$\phi3$', mesh=mesh, value=z3)
phi4 = CellVariable(name='$\phi4$', mesh=mesh, value=z4)
# Boundary Conditions
facesTopRight = ((mesh.getFacesRight())
| (mesh.getFacesTop()))
facesBottomLeft = ((mesh.getFacesLeft())
| (mesh.getFacesBottom()))
BCs = (FixedFlux(faces=mesh.getFacesRight(), value=0.),
FixedFlux(faces=mesh.getFacesLeft(), value=0.))
# Viewer config
if __name__ == '__main__':
viewer = Matplotlib2DGridContourViewer(vars=phi1+phi2+phi3+phi4,
limits={'xmin': -l/4, 'xmax': l/4, 'ymin': -l/4, 'ymax': l/4},
datamin=0.,
axes=ax1,
legend=None)
# equation config
D = 0.1
# define eqs
eq1 = (TransientTerm(var=phi1) == DiffusionTerm(coeff = D, var=phi1))
eq2 = (TransientTerm(var=phi2) == DiffusionTerm(coeff = D, var=phi2))
eq3 = (TransientTerm(var=phi3) == DiffusionTerm(coeff = D, var=phi3))
eq4 = (TransientTerm(var=phi4) == DiffusionTerm(coeff = D, var=phi4))
eq = (eq1 & eq2 & eq3 & eq4)
# temporal config
dt = 0.05
steps = 15
txt.set_text(r'Time Step 0')
for step in range(steps):
print 'steps:', step+1
eq.solve(boundaryConditions=BCs, dt=dt)
txt.set_text(r'Time Step %s'%(step+1))
ax1.cla()
levels=np.arange(0.,.8,.1)
cs1=ax1.contour(X, Y, phi1.value.reshape(201,-1), levels, colors='k', alpha=.3)
cs2=ax1.contour(X, Y, phi2.value.reshape(201,-1), levels, colors='k', alpha=.3)
cs3=ax1.contour(X, Y, phi3.value.reshape(201,-1), levels, colors='k', alpha=.3)
cs4=ax1.contour(X, Y, phi4.value.reshape(201,-1), levels, colors='k', alpha=.3)
ax1.clabel(cs1, fontsize=9, inline=1)
ax1.clabel(cs2, fontsize=9, inline=1)
ax1.clabel(cs3, fontsize=9, inline=1)
ax1.clabel(cs4, fontsize=9, inline=1)
viewer.plot()___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: how to create large grid with holes and odd geometry
Jonathan and Daniel -- Thanks. In a few days or so, I shall post what I come up with for others to use. Jamie On Wed, Jun 15, 2016 at 3:58 PM, Guyer, Jonathan E. Dr. (Fed) < jonathan.gu...@nist.gov> wrote: > You can generate a simple mesh in order to probe what's stored by doing > something like: > > mesh1 = fp.Grid2D(nx=2, ny=2) > mesh2 = fp.Grid2D(nx=2, ny=2) + [[2], [1]] > mesh = mesh1 + mesh2 > > > > On Jun 15, 2016, at 2:06 PM, James Pringle wrote: > > > > Jonathan -- > > > > Thank you, this is exactly what I need. Two more questions. > > • Is the order of the vertices in faceVertexIDs important? > > • Is the order of faces in cellFaceIDs important? Must they wind > clockwise or counterclockwise? > > Thanks, > > Jamie > > > > On Wed, Jun 15, 2016 at 1:56 PM, Guyer, Jonathan E. Dr. (Fed) < > jonathan.gu...@nist.gov> wrote: > > The first three arguments to the Mesh2D constructor are (all that is) > required: > > > > class Mesh2D(Mesh): > > def __init__(self, vertexCoords, faceVertexIDs, cellFaceIDs, ...): > > > > All other arguments have default values assigned. > > > > For your case: > > > > vertexCoords is of shape (2, N) where N is the number of vertices > > faceVertexIDs is of shape (2, M) where M is the number of faces > > cellFaceIDs is of shape (3, P) where P is the number of cells > > > > faceVertexIDs and cellFaceIDs are 0-based, as they are indices into the > preceding array > > > > > On Jun 15, 2016, at 1:29 PM, James Pringle wrote: > > > > > > Well, I am motivated to give it a go, since I only have the summer to > make progress on project and it is blocking my research progress right now. > Can you give me a pointer to where the appropriate quantities are defined? > I can certainly write code to make the transformations, but it is a bit > hard without understanding precisely what is defined in the Mesh2D object. > I have made a simple Mesh2D object, but I am not sure which of the > attributes, etc, are absolutely required, or how they are precisely defined. > > > > > > Perhaps most useful for me would be > > > • a definition of which parts of the Mesh2D object must exist > > > • and the format of those parts, in particular the a face to > vertex array and a > > > cell to face array > > > Or you could point me to the appropriate part of the Gmsh code so I > can go from there. I presume poking around in fipy.Gmsh2D would be a good > place to start? Is there a better place to start? > > > > > > I would love any documentation on the details of the Mesh2D object. > > > > > > Thanks, > > > Jamie Pringle > > > University of New Hampshire > > > > > > On Wed, Jun 15, 2016 at 12:21 PM, Daniel Wheeler < > daniel.wheel...@gmail.com> wrote: > > > Hi Jamie, > > > > > > There is no automated way to make a FiPy mesh from Scipy's Delaunay > > > object. The problem is that FiPy needs a face to vertex array and a > > > cell to face array while the Delaunay object just has the cell to > > > vertex array. The functionality to extract the face to vertex array > > > and cell to face array is in FiPy because it must be done for Gmsh, > > > however, it is not abstracted in so that it can be reused. > > > > > > It is certainly possible to make the correct arrays with Numpy and > > > pass them to the Mesh2D object, but it's a bit of work to write the > > > code. If I find some time I might give it a go, but I don't know when > > > I will get to it. > > > > > > Cheers, > > > > > > Daniel > > > > > > On Tue, Jun 14, 2016 at 4:15 PM, James Pringle > wrote: > > > > Daniel et al. -- > > > > > > > > As referenced earlier in this thread, I have a complex domain > with holes > > > > in it; I now have it broken up into triangles, based on the Delaunay > package > > > > in SciPy. I have > > > > > > > > Locations of all vertex coordinates, > > > > list of vertices that make up faces > > > > list of faces that make cells > > > > list of faces that make up all the internal and external boundaries. > > > > > > > > Can you point me to any code or documentation that would help me > understand > > > > how to make this into a Mesh2D object? I am having a devil of a time > > > > figuring out from the manual online. The best would be something > that used > > > > the output of either the triangles or scipy.spatial.Delaunay() > packaged. My > > > > equation is of the form > > > > > > > > 0=J(Psi,A(x,y)) + \Del(B(x,y)*\Del Psi) > > > > > > > > > > > > and I can get the coefficients A(x,y) and B(x,y) on either the faces > or in > > > > the cell centers are needed. > > > > > > > > Thank you. > > > > Jamie Pringle > > > > University of New Hampshire > > > > > > -- > > > Daniel Wheeler > > > ___ > > > fipy mailing list > > > fipy@nist.gov > > > > 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=Q6M6XDjzGDlxMDgrY-riZ7Q70Abl5Aq7SRrTwZNYsws&s=yiO9QA
Re: how to create large grid with holes and odd geometry
You can generate a simple mesh in order to probe what's stored by doing something like: mesh1 = fp.Grid2D(nx=2, ny=2) mesh2 = fp.Grid2D(nx=2, ny=2) + [[2], [1]] mesh = mesh1 + mesh2 > On Jun 15, 2016, at 2:06 PM, James Pringle wrote: > > Jonathan -- > > Thank you, this is exactly what I need. Two more questions. > • Is the order of the vertices in faceVertexIDs important? > • Is the order of faces in cellFaceIDs important? Must they wind > clockwise or counterclockwise? > Thanks, > Jamie > > On Wed, Jun 15, 2016 at 1:56 PM, Guyer, Jonathan E. Dr. (Fed) > wrote: > The first three arguments to the Mesh2D constructor are (all that is) > required: > > class Mesh2D(Mesh): > def __init__(self, vertexCoords, faceVertexIDs, cellFaceIDs, ...): > > All other arguments have default values assigned. > > For your case: > > vertexCoords is of shape (2, N) where N is the number of vertices > faceVertexIDs is of shape (2, M) where M is the number of faces > cellFaceIDs is of shape (3, P) where P is the number of cells > > faceVertexIDs and cellFaceIDs are 0-based, as they are indices into the > preceding array > > > On Jun 15, 2016, at 1:29 PM, James Pringle wrote: > > > > Well, I am motivated to give it a go, since I only have the summer to make > > progress on project and it is blocking my research progress right now. Can > > you give me a pointer to where the appropriate quantities are defined? I > > can certainly write code to make the transformations, but it is a bit hard > > without understanding precisely what is defined in the Mesh2D object. I > > have made a simple Mesh2D object, but I am not sure which of the > > attributes, etc, are absolutely required, or how they are precisely defined. > > > > Perhaps most useful for me would be > > • a definition of which parts of the Mesh2D object must exist > > • and the format of those parts, in particular the a face to vertex > > array and a > > cell to face array > > Or you could point me to the appropriate part of the Gmsh code so I can go > > from there. I presume poking around in fipy.Gmsh2D would be a good place to > > start? Is there a better place to start? > > > > I would love any documentation on the details of the Mesh2D object. > > > > Thanks, > > Jamie Pringle > > University of New Hampshire > > > > On Wed, Jun 15, 2016 at 12:21 PM, Daniel Wheeler > > wrote: > > Hi Jamie, > > > > There is no automated way to make a FiPy mesh from Scipy's Delaunay > > object. The problem is that FiPy needs a face to vertex array and a > > cell to face array while the Delaunay object just has the cell to > > vertex array. The functionality to extract the face to vertex array > > and cell to face array is in FiPy because it must be done for Gmsh, > > however, it is not abstracted in so that it can be reused. > > > > It is certainly possible to make the correct arrays with Numpy and > > pass them to the Mesh2D object, but it's a bit of work to write the > > code. If I find some time I might give it a go, but I don't know when > > I will get to it. > > > > Cheers, > > > > Daniel > > > > On Tue, Jun 14, 2016 at 4:15 PM, James Pringle wrote: > > > Daniel et al. -- > > > > > > As referenced earlier in this thread, I have a complex domain with > > > holes > > > in it; I now have it broken up into triangles, based on the Delaunay > > > package > > > in SciPy. I have > > > > > > Locations of all vertex coordinates, > > > list of vertices that make up faces > > > list of faces that make cells > > > list of faces that make up all the internal and external boundaries. > > > > > > Can you point me to any code or documentation that would help me > > > understand > > > how to make this into a Mesh2D object? I am having a devil of a time > > > figuring out from the manual online. The best would be something that used > > > the output of either the triangles or scipy.spatial.Delaunay() packaged. > > > My > > > equation is of the form > > > > > > 0=J(Psi,A(x,y)) + \Del(B(x,y)*\Del Psi) > > > > > > > > > and I can get the coefficients A(x,y) and B(x,y) on either the faces or in > > > the cell centers are needed. > > > > > > Thank you. > > > Jamie Pringle > > > University of New Hampshire > > > > -- > > Daniel Wheeler > > ___ > > fipy mailing list > > fipy@nist.gov > > 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=Q6M6XDjzGDlxMDgrY-riZ7Q70Abl5Aq7SRrTwZNYsws&s=yiO9QAq_3EUy9lVY5F76Tc9mnaB5Ubclilpum_QCUOE&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=Q6M6XDjzGDlxMDgrY-riZ7Q70Abl5Aq7SRrTwZNYsws&s=yIJJw-iQknimZfYCRrZmk7V2eHSJrCtcAehgVbGeiDk&e= > > ] > > > > ___ > >
Re: how to create large grid with holes and odd geometry
The answer to both is "I don't think so" > On Jun 15, 2016, at 2:06 PM, James Pringle wrote: > > Jonathan -- > > Thank you, this is exactly what I need. Two more questions. > • Is the order of the vertices in faceVertexIDs important? > • Is the order of faces in cellFaceIDs important? Must they wind > clockwise or counterclockwise? > Thanks, > Jamie > > On Wed, Jun 15, 2016 at 1:56 PM, Guyer, Jonathan E. Dr. (Fed) > wrote: > The first three arguments to the Mesh2D constructor are (all that is) > required: > > class Mesh2D(Mesh): > def __init__(self, vertexCoords, faceVertexIDs, cellFaceIDs, ...): > > All other arguments have default values assigned. > > For your case: > > vertexCoords is of shape (2, N) where N is the number of vertices > faceVertexIDs is of shape (2, M) where M is the number of faces > cellFaceIDs is of shape (3, P) where P is the number of cells > > faceVertexIDs and cellFaceIDs are 0-based, as they are indices into the > preceding array > > > On Jun 15, 2016, at 1:29 PM, James Pringle wrote: > > > > Well, I am motivated to give it a go, since I only have the summer to make > > progress on project and it is blocking my research progress right now. Can > > you give me a pointer to where the appropriate quantities are defined? I > > can certainly write code to make the transformations, but it is a bit hard > > without understanding precisely what is defined in the Mesh2D object. I > > have made a simple Mesh2D object, but I am not sure which of the > > attributes, etc, are absolutely required, or how they are precisely defined. > > > > Perhaps most useful for me would be > > • a definition of which parts of the Mesh2D object must exist > > • and the format of those parts, in particular the a face to vertex > > array and a > > cell to face array > > Or you could point me to the appropriate part of the Gmsh code so I can go > > from there. I presume poking around in fipy.Gmsh2D would be a good place to > > start? Is there a better place to start? > > > > I would love any documentation on the details of the Mesh2D object. > > > > Thanks, > > Jamie Pringle > > University of New Hampshire > > > > On Wed, Jun 15, 2016 at 12:21 PM, Daniel Wheeler > > wrote: > > Hi Jamie, > > > > There is no automated way to make a FiPy mesh from Scipy's Delaunay > > object. The problem is that FiPy needs a face to vertex array and a > > cell to face array while the Delaunay object just has the cell to > > vertex array. The functionality to extract the face to vertex array > > and cell to face array is in FiPy because it must be done for Gmsh, > > however, it is not abstracted in so that it can be reused. > > > > It is certainly possible to make the correct arrays with Numpy and > > pass them to the Mesh2D object, but it's a bit of work to write the > > code. If I find some time I might give it a go, but I don't know when > > I will get to it. > > > > Cheers, > > > > Daniel > > > > On Tue, Jun 14, 2016 at 4:15 PM, James Pringle wrote: > > > Daniel et al. -- > > > > > > As referenced earlier in this thread, I have a complex domain with > > > holes > > > in it; I now have it broken up into triangles, based on the Delaunay > > > package > > > in SciPy. I have > > > > > > Locations of all vertex coordinates, > > > list of vertices that make up faces > > > list of faces that make cells > > > list of faces that make up all the internal and external boundaries. > > > > > > Can you point me to any code or documentation that would help me > > > understand > > > how to make this into a Mesh2D object? I am having a devil of a time > > > figuring out from the manual online. The best would be something that used > > > the output of either the triangles or scipy.spatial.Delaunay() packaged. > > > My > > > equation is of the form > > > > > > 0=J(Psi,A(x,y)) + \Del(B(x,y)*\Del Psi) > > > > > > > > > and I can get the coefficients A(x,y) and B(x,y) on either the faces or in > > > the cell centers are needed. > > > > > > Thank you. > > > Jamie Pringle > > > University of New Hampshire > > > > -- > > Daniel Wheeler > > ___ > > fipy mailing list > > fipy@nist.gov > > 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=Q6M6XDjzGDlxMDgrY-riZ7Q70Abl5Aq7SRrTwZNYsws&s=yiO9QAq_3EUy9lVY5F76Tc9mnaB5Ubclilpum_QCUOE&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=Q6M6XDjzGDlxMDgrY-riZ7Q70Abl5Aq7SRrTwZNYsws&s=yIJJw-iQknimZfYCRrZmk7V2eHSJrCtcAehgVbGeiDk&e= > > ] > > > > ___ > > fipy mailing list > > fipy@nist.gov > > https://urldefense.proofpoint.com/v2/url?u=http-3A__www.ctcms.nist.gov_fipy&d=CwIGaQ&c=c6MrceVCY5m5A_KAUkrdoA&
Re: how to create large grid with holes and odd geometry
I think that you need to instatiate one of these https://github.com/usnistgov/fipy/blob/develop/fipy/meshes/mesh2D.py#L70 with the correct arrays. For examples, >>> points = [[0,0], [0, 1], [1, 0], [1, 1]] >>> from scipy.spatial import Delaunay >>> tri = Delaunay(points) >>> tri.points array([[ 0., 0.], [ 0., 1.], [ 1., 0.], [ 1., 1.]]) >>> tri.simplices array([[3, 1, 0], [2, 3, 0]], dtype=int32) So the vertices are just the points but reshaped to (2, N) from (N, 2). The face to vertex array would be [[3, 1, 0, 2, 0 [1, 0, 3, 3, 2]] Its all the faces or edges in 2D, and the cell to face ID would be [[0, 3], [1, 2], [2, 4]] Then you can instatiate the Mesh2D object with those three arrrays. Look at https://github.com/usnistgov/fipy/blob/develop/fipy/meshes/mesh2D.py#L269 to see how to do it as well. Notice that the test example has both quads and triangles (the -1 values are masked). On Wed, Jun 15, 2016 at 1:29 PM, James Pringle wrote: > Well, I am motivated to give it a go, since I only have the summer to make > progress on project and it is blocking my research progress right now. Can > you give me a pointer to where the appropriate quantities are defined? I can > certainly write code to make the transformations, but it is a bit hard > without understanding precisely what is defined in the Mesh2D object. I have > made a simple Mesh2D object, but I am not sure which of the attributes, etc, > are absolutely required, or how they are precisely defined. > > Perhaps most useful for me would be > > a definition of which parts of the Mesh2D object must exist > and the format of those parts, in particular the a face to vertex array and > a > cell to face array > > Or you could point me to the appropriate part of the Gmsh code so I can go > from there. I presume poking around in fipy.Gmsh2D would be a good place to > start? Is there a better place to start? > > I would love any documentation on the details of the Mesh2D object. > > Thanks, > Jamie Pringle > University of New Hampshire > > On Wed, Jun 15, 2016 at 12:21 PM, Daniel Wheeler > wrote: >> >> Hi Jamie, >> >> There is no automated way to make a FiPy mesh from Scipy's Delaunay >> object. The problem is that FiPy needs a face to vertex array and a >> cell to face array while the Delaunay object just has the cell to >> vertex array. The functionality to extract the face to vertex array >> and cell to face array is in FiPy because it must be done for Gmsh, >> however, it is not abstracted in so that it can be reused. >> >> It is certainly possible to make the correct arrays with Numpy and >> pass them to the Mesh2D object, but it's a bit of work to write the >> code. If I find some time I might give it a go, but I don't know when >> I will get to it. >> >> Cheers, >> >> Daniel >> >> On Tue, Jun 14, 2016 at 4:15 PM, James Pringle wrote: >> > Daniel et al. -- >> > >> > As referenced earlier in this thread, I have a complex domain with >> > holes >> > in it; I now have it broken up into triangles, based on the Delaunay >> > package >> > in SciPy. I have >> > >> > Locations of all vertex coordinates, >> > list of vertices that make up faces >> > list of faces that make cells >> > list of faces that make up all the internal and external boundaries. >> > >> > Can you point me to any code or documentation that would help me >> > understand >> > how to make this into a Mesh2D object? I am having a devil of a time >> > figuring out from the manual online. The best would be something that >> > used >> > the output of either the triangles or scipy.spatial.Delaunay() packaged. >> > My >> > equation is of the form >> > >> > 0=J(Psi,A(x,y)) + \Del(B(x,y)*\Del Psi) >> > >> > >> > and I can get the coefficients A(x,y) and B(x,y) on either the faces or >> > in >> > the cell centers are needed. >> > >> > Thank you. >> > Jamie Pringle >> > University of New Hampshire >> >> -- >> Daniel Wheeler >> ___ >> fipy mailing list >> fipy@nist.gov >> >> 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=Q6M6XDjzGDlxMDgrY-riZ7Q70Abl5Aq7SRrTwZNYsws&s=yiO9QAq_3EUy9lVY5F76Tc9mnaB5Ubclilpum_QCUOE&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=Q6M6XDjzGDlxMDgrY-riZ7Q70Abl5Aq7SRrTwZNYsws&s=yIJJw-iQknimZfYCRrZmk7V2eHSJrCtcAehgVbGeiDk&e= >> ] > > > > ___ > 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
Re: how to create large grid with holes and odd geometry
Jonathan -- Thank you, this is exactly what I need. Two more questions. 1. Is the order of the vertices in faceVertexIDs important? 2. Is the order of faces in cellFaceIDs important? Must they wind clockwise or counterclockwise? Thanks, Jamie On Wed, Jun 15, 2016 at 1:56 PM, Guyer, Jonathan E. Dr. (Fed) < jonathan.gu...@nist.gov> wrote: > The first three arguments to the Mesh2D constructor are (all that is) > required: > > class Mesh2D(Mesh): > def __init__(self, vertexCoords, faceVertexIDs, cellFaceIDs, ...): > > All other arguments have default values assigned. > > For your case: > > vertexCoords is of shape (2, N) where N is the number of vertices > faceVertexIDs is of shape (2, M) where M is the number of faces > cellFaceIDs is of shape (3, P) where P is the number of cells > > faceVertexIDs and cellFaceIDs are 0-based, as they are indices into the > preceding array > > > On Jun 15, 2016, at 1:29 PM, James Pringle wrote: > > > > Well, I am motivated to give it a go, since I only have the summer to > make progress on project and it is blocking my research progress right now. > Can you give me a pointer to where the appropriate quantities are defined? > I can certainly write code to make the transformations, but it is a bit > hard without understanding precisely what is defined in the Mesh2D object. > I have made a simple Mesh2D object, but I am not sure which of the > attributes, etc, are absolutely required, or how they are precisely defined. > > > > Perhaps most useful for me would be > > • a definition of which parts of the Mesh2D object must exist > > • and the format of those parts, in particular the a face to > vertex array and a > > cell to face array > > Or you could point me to the appropriate part of the Gmsh code so I can > go from there. I presume poking around in fipy.Gmsh2D would be a good place > to start? Is there a better place to start? > > > > I would love any documentation on the details of the Mesh2D object. > > > > Thanks, > > Jamie Pringle > > University of New Hampshire > > > > On Wed, Jun 15, 2016 at 12:21 PM, Daniel Wheeler < > daniel.wheel...@gmail.com> wrote: > > Hi Jamie, > > > > There is no automated way to make a FiPy mesh from Scipy's Delaunay > > object. The problem is that FiPy needs a face to vertex array and a > > cell to face array while the Delaunay object just has the cell to > > vertex array. The functionality to extract the face to vertex array > > and cell to face array is in FiPy because it must be done for Gmsh, > > however, it is not abstracted in so that it can be reused. > > > > It is certainly possible to make the correct arrays with Numpy and > > pass them to the Mesh2D object, but it's a bit of work to write the > > code. If I find some time I might give it a go, but I don't know when > > I will get to it. > > > > Cheers, > > > > Daniel > > > > On Tue, Jun 14, 2016 at 4:15 PM, James Pringle wrote: > > > Daniel et al. -- > > > > > > As referenced earlier in this thread, I have a complex domain with > holes > > > in it; I now have it broken up into triangles, based on the Delaunay > package > > > in SciPy. I have > > > > > > Locations of all vertex coordinates, > > > list of vertices that make up faces > > > list of faces that make cells > > > list of faces that make up all the internal and external boundaries. > > > > > > Can you point me to any code or documentation that would help me > understand > > > how to make this into a Mesh2D object? I am having a devil of a time > > > figuring out from the manual online. The best would be something that > used > > > the output of either the triangles or scipy.spatial.Delaunay() > packaged. My > > > equation is of the form > > > > > > 0=J(Psi,A(x,y)) + \Del(B(x,y)*\Del Psi) > > > > > > > > > and I can get the coefficients A(x,y) and B(x,y) on either the faces > or in > > > the cell centers are needed. > > > > > > Thank you. > > > Jamie Pringle > > > University of New Hampshire > > > > -- > > Daniel Wheeler > > ___ > > fipy mailing list > > fipy@nist.gov > > > 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=Q6M6XDjzGDlxMDgrY-riZ7Q70Abl5Aq7SRrTwZNYsws&s=yiO9QAq_3EUy9lVY5F76Tc9mnaB5Ubclilpum_QCUOE&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=Q6M6XDjzGDlxMDgrY-riZ7Q70Abl5Aq7SRrTwZNYsws&s=yIJJw-iQknimZfYCRrZmk7V2eHSJrCtcAehgVbGeiDk&e= > ] > > > > ___ > > fipy mailing list > > fipy@nist.gov > > > https://urldefense.proofpoint.com/v2/url?u=http-3A__www.ctcms.nist.gov_fipy&d=CwIGaQ&c=c6MrceVCY5m5A_KAUkrdoA&r=oMf1-WHpUHeD7kN3S7612CDdF2TDPqDF9R-n-71Ks1Y&m=iGk4QC0NEBSkRipnXJcMmV0qOdCYUWc6oQMIUpAmqyI&s=1tQAVl8UCyd9BiwQBhF_rRh-tByp33v
Re: how to create large grid with holes and odd geometry
The first three arguments to the Mesh2D constructor are (all that is) required: class Mesh2D(Mesh): def __init__(self, vertexCoords, faceVertexIDs, cellFaceIDs, ...): All other arguments have default values assigned. For your case: vertexCoords is of shape (2, N) where N is the number of vertices faceVertexIDs is of shape (2, M) where M is the number of faces cellFaceIDs is of shape (3, P) where P is the number of cells faceVertexIDs and cellFaceIDs are 0-based, as they are indices into the preceding array > On Jun 15, 2016, at 1:29 PM, James Pringle wrote: > > Well, I am motivated to give it a go, since I only have the summer to make > progress on project and it is blocking my research progress right now. Can > you give me a pointer to where the appropriate quantities are defined? I can > certainly write code to make the transformations, but it is a bit hard > without understanding precisely what is defined in the Mesh2D object. I have > made a simple Mesh2D object, but I am not sure which of the attributes, etc, > are absolutely required, or how they are precisely defined. > > Perhaps most useful for me would be > • a definition of which parts of the Mesh2D object must exist > • and the format of those parts, in particular the a face to vertex > array and a > cell to face array > Or you could point me to the appropriate part of the Gmsh code so I can go > from there. I presume poking around in fipy.Gmsh2D would be a good place to > start? Is there a better place to start? > > I would love any documentation on the details of the Mesh2D object. > > Thanks, > Jamie Pringle > University of New Hampshire > > On Wed, Jun 15, 2016 at 12:21 PM, Daniel Wheeler > wrote: > Hi Jamie, > > There is no automated way to make a FiPy mesh from Scipy's Delaunay > object. The problem is that FiPy needs a face to vertex array and a > cell to face array while the Delaunay object just has the cell to > vertex array. The functionality to extract the face to vertex array > and cell to face array is in FiPy because it must be done for Gmsh, > however, it is not abstracted in so that it can be reused. > > It is certainly possible to make the correct arrays with Numpy and > pass them to the Mesh2D object, but it's a bit of work to write the > code. If I find some time I might give it a go, but I don't know when > I will get to it. > > Cheers, > > Daniel > > On Tue, Jun 14, 2016 at 4:15 PM, James Pringle wrote: > > Daniel et al. -- > > > > As referenced earlier in this thread, I have a complex domain with holes > > in it; I now have it broken up into triangles, based on the Delaunay package > > in SciPy. I have > > > > Locations of all vertex coordinates, > > list of vertices that make up faces > > list of faces that make cells > > list of faces that make up all the internal and external boundaries. > > > > Can you point me to any code or documentation that would help me understand > > how to make this into a Mesh2D object? I am having a devil of a time > > figuring out from the manual online. The best would be something that used > > the output of either the triangles or scipy.spatial.Delaunay() packaged. My > > equation is of the form > > > > 0=J(Psi,A(x,y)) + \Del(B(x,y)*\Del Psi) > > > > > > and I can get the coefficients A(x,y) and B(x,y) on either the faces or in > > the cell centers are needed. > > > > Thank you. > > Jamie Pringle > > University of New Hampshire > > -- > Daniel Wheeler > ___ > fipy mailing list > fipy@nist.gov > 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=Q6M6XDjzGDlxMDgrY-riZ7Q70Abl5Aq7SRrTwZNYsws&s=yiO9QAq_3EUy9lVY5F76Tc9mnaB5Ubclilpum_QCUOE&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=Q6M6XDjzGDlxMDgrY-riZ7Q70Abl5Aq7SRrTwZNYsws&s=yIJJw-iQknimZfYCRrZmk7V2eHSJrCtcAehgVbGeiDk&e= > ] > > ___ > fipy mailing list > fipy@nist.gov > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: how to create large grid with holes and odd geometry
Well, I am motivated to give it a go, since I only have the summer to make progress on project and it is blocking my research progress right now. Can you give me a pointer to where the appropriate quantities are defined? I can certainly write code to make the transformations, but it is a bit hard without understanding precisely what is defined in the Mesh2D object. I have made a simple Mesh2D object, but I am not sure which of the attributes, etc, are absolutely required, or how they are precisely defined. Perhaps most useful for me would be 1. a definition of which parts of the Mesh2D object must exist 2. and the format of those parts, in particular the a face to vertex array and a cell to face array Or you could point me to the appropriate part of the Gmsh code so I can go from there. I presume poking around in fipy.Gmsh2D would be a good place to start? Is there a better place to start? I would love any documentation on the details of the Mesh2D object. Thanks, Jamie Pringle University of New Hampshire On Wed, Jun 15, 2016 at 12:21 PM, Daniel Wheeler wrote: > Hi Jamie, > > There is no automated way to make a FiPy mesh from Scipy's Delaunay > object. The problem is that FiPy needs a face to vertex array and a > cell to face array while the Delaunay object just has the cell to > vertex array. The functionality to extract the face to vertex array > and cell to face array is in FiPy because it must be done for Gmsh, > however, it is not abstracted in so that it can be reused. > > It is certainly possible to make the correct arrays with Numpy and > pass them to the Mesh2D object, but it's a bit of work to write the > code. If I find some time I might give it a go, but I don't know when > I will get to it. > > Cheers, > > Daniel > > On Tue, Jun 14, 2016 at 4:15 PM, James Pringle wrote: > > Daniel et al. -- > > > > As referenced earlier in this thread, I have a complex domain with > holes > > in it; I now have it broken up into triangles, based on the Delaunay > package > > in SciPy. I have > > > > Locations of all vertex coordinates, > > list of vertices that make up faces > > list of faces that make cells > > list of faces that make up all the internal and external boundaries. > > > > Can you point me to any code or documentation that would help me > understand > > how to make this into a Mesh2D object? I am having a devil of a time > > figuring out from the manual online. The best would be something that > used > > the output of either the triangles or scipy.spatial.Delaunay() > packaged. My > > equation is of the form > > > > 0=J(Psi,A(x,y)) + \Del(B(x,y)*\Del Psi) > > > > > > and I can get the coefficients A(x,y) and B(x,y) on either the faces or > in > > the cell centers are needed. > > > > Thank you. > > Jamie Pringle > > University of New Hampshire > > -- > Daniel Wheeler > ___ > fipy mailing list > fipy@nist.gov > > 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=Q6M6XDjzGDlxMDgrY-riZ7Q70Abl5Aq7SRrTwZNYsws&s=yiO9QAq_3EUy9lVY5F76Tc9mnaB5Ubclilpum_QCUOE&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=Q6M6XDjzGDlxMDgrY-riZ7Q70Abl5Aq7SRrTwZNYsws&s=yIJJw-iQknimZfYCRrZmk7V2eHSJrCtcAehgVbGeiDk&e= > ] > ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: casting implicit Boundary Conditions in FiPy
On Wed, Jun 15, 2016 at 7:27 AM, Gopalakrishnan, Krishnakumar wrote: > > Dan. I was able validate that your code correctly implements the Implicit > Neumann Boundary Condition (in 1-D). > > Here is a link to the plot of the solution after a transient simulation for > 0.2 seconds. The plot on the left is generated in Mathematica, and that to > the right is that generated by the FiPy’s matplotlib viewer class. The two > are identical. > > https://imperialcollegelondon.box.com/s/lb8v1iqxb3rqhxsphn9cmhwcmh5jzpiz Awesome that it worked. I think that you're right though regarding first versus second order boundary condition. I think it is first order in space. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: how to create large grid with holes and odd geometry
Hi Jamie, There is no automated way to make a FiPy mesh from Scipy's Delaunay object. The problem is that FiPy needs a face to vertex array and a cell to face array while the Delaunay object just has the cell to vertex array. The functionality to extract the face to vertex array and cell to face array is in FiPy because it must be done for Gmsh, however, it is not abstracted in so that it can be reused. It is certainly possible to make the correct arrays with Numpy and pass them to the Mesh2D object, but it's a bit of work to write the code. If I find some time I might give it a go, but I don't know when I will get to it. Cheers, Daniel On Tue, Jun 14, 2016 at 4:15 PM, James Pringle wrote: > Daniel et al. -- > > As referenced earlier in this thread, I have a complex domain with holes > in it; I now have it broken up into triangles, based on the Delaunay package > in SciPy. I have > > Locations of all vertex coordinates, > list of vertices that make up faces > list of faces that make cells > list of faces that make up all the internal and external boundaries. > > Can you point me to any code or documentation that would help me understand > how to make this into a Mesh2D object? I am having a devil of a time > figuring out from the manual online. The best would be something that used > the output of either the triangles or scipy.spatial.Delaunay() packaged. My > equation is of the form > > 0=J(Psi,A(x,y)) + \Del(B(x,y)*\Del Psi) > > > and I can get the coefficients A(x,y) and B(x,y) on either the faces or in > the cell centers are needed. > > Thank you. > Jamie Pringle > University of New Hampshire -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: casting implicit Boundary Conditions in FiPy
> On Jun 15, 2016, at 7:27 AM, Gopalakrishnan, Krishnakumar > wrote: > > Can we add this as a feature request in the project’s github page, perhaps > for a future release? By all means, please file issues on GitHub for any features you'd like or bugs you find. ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
RE: casting implicit Boundary Conditions in FiPy
Hi, Dan. I was able validate that your code correctly implements the Implicit Neumann Boundary Condition (in 1-D). Here is a link to the plot of the solution after a transient simulation for 0.2 seconds. The plot on the left is generated in Mathematica, and that to the right is that generated by the FiPy’s matplotlib viewer class. The two are identical. https://imperialcollegelondon.box.com/s/lb8v1iqxb3rqhxsphn9cmhwcmh5jzpiz Can we add this as a feature request in the project’s github page, perhaps for a future release? Does Implicit Neumann BCs occur often in other disciplines ? Once again, thanks for your help in providing this solution approach. Krishna From: fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] On Behalf Of Daniel Wheeler Sent: 10 June 2016 18:18 To: Multiple recipients of list Subject: Re: casting implicit Boundary Conditions in FiPy Hi Krishina and Ray, Thanks for the interesting discussion. I'm not 100% sure about everything that Krishina is asking for in the latter part of the discussion so I'm just going to address the code that Ray has developed below (my code is below). I think there is a way to handle the right hand side boundary condition both implicitly and with second order accuracy using an ImplicitSourceTerm boundary condition trick. Sort of similar to what is here, http://www.ctcms.nist.gov/fipy/documentation/USAGE.html#applying-fixed-flux-boundary-conditions. Like I said, I think the code below is both second order accurate (for fixed dx) and implicit. Extending this to 2D might raise a few issues. The grid spacing needs to be in the coefficient so, obviously, dx needs to be fixed in both x and y directions. Also, I'm not sure if I'm missing a factor of dx in the source term in 1D as I'm using both a divergence and an ImplicitSourceTerm so there is some question about volumes and face areas in 2D as well. I'm also confused about the signs, I had to flip the sign in front of the source to make it work. It seems to look right though in 1D and you can just take one massive time step to get to the answer. This will only work if k is negative otherwise it's unstable, right? The way I came up with the source was the following n.grad(phi) = k * (phi_P + n.grad(phi) * dx / 2) and then solve for n.grad(phi) which gives n.grad(phi) = k * phi_P / (1 - dx * k / 2) where phi_P is the value of phi at the cell next to the boundary with the implicit boundary condition. Then we fake the outbound flux to be the expression on the right of the equation. Just as a general note it would be great in FiPy if we could come up with a nice way to write boundary conditions in scripts that did all these tricks implicitly without having to know all these background details about FV and how FiPy works. import fipy as fp nx = 50 dx = 1. mesh = fp.Grid1D(nx=nx, dx=dx) phi = fp.CellVariable(name="field variable", mesh=mesh, value=1.0) D = 1. k = -1. diffCoeff = fp.FaceVariable(mesh=mesh, value=D) diffCoeff.constrain(0., mesh.facesRight) valueLeft = 0.0 phi.constrain(valueLeft, mesh.facesLeft) #phi.faceGrad.constrain([phi], mesh.facesRight) # This is the problematic BC #phi.faceGrad.constrain(phi.harmonicFaceValue, mesh.facesRight) # This is the problematic BC #phi.faceGrad.constrain([k * phi.harmonicFaceValue], mesh.facesRight) # This is the problematic BC implicitCoeff = -D * k / (1. - k * dx / 2.) eq = (fp.TransientTerm() == fp.DiffusionTerm(diffCoeff) - \ fp.ImplicitSourceTerm((mesh.faceNormals * implicitCoeff * mesh.facesRight).divergence)) timeStep = 0.9 * dx**2 / (2 * D) timeStep = 10.0 steps = 800 viewer = fp.Viewer(vars=phi, datamax=1., datamin=0.) for step in range(steps): eq.solve(var=phi, dt=timeStep) viewer.plot() On Thu, Jun 9, 2016 at 12:02 PM, Raymond Smith mailto:smit...@mit.edu>> wrote: Hi, Krishna. Perhaps I'm misunderstanding something, but I'm still not convinced the second version you suggested -- c.faceGrad.constrain([-(j_at_c_star + partial_j_at_op_point*(c.faceValue - c_star))], mesh.facesTop) -- isn't working like you want. Could you look at the example I suggested to see if that behaves differently than you expect? Here's the code I used. To me it looks very similar in form to c constraint above and at first glance it seems to behave exactly like we want -- that is, throughout the time stepping, n*grad(phi) is constrained to the value, -phi at the surface. Correct me if I'm wrong, but my impression is that this is the behavior you desire. from fipy import * nx = 50 dx = 1. mesh = Grid1D(nx=nx, dx=dx) phi = CellVariable(name="field variable", mesh=mesh, value=1.0) D = 1. valueLeft = 0.0 phi.constrain(valueLeft, mesh.facesLeft) #phi.faceGrad.constrain([phi], mesh.facesRight) # This is the problematic BC #phi.faceGrad.constrain(phi.harmonicFaceValue, mesh.facesRight) # This is the problematic BC phi.faceGrad.constrain([-phi.harmonicFaceValue], mesh.facesRight) # This is the pr
