Reputation: 49
I successfully set up a very simple diffusion problem including a source face (continuum) and sink faces (sidewall, top) on a mesh generated with GMSH.
from fipy import *
mesh = Gmsh2D('''Point(1) = {0, 0, 0, 1.0};
Point(2) = {12, 0, 0, 0.1};
Point(3) = {12, 16, 0, 0.1};
Point(4) = {15, 16, 0, 0.1};
Point(5) = {15, 100, 0, 1.0};
Point(6) = {0, 100, 0, 1.0};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 5};
Line(5) = {5, 6};
Line(6) = {6, 1};
Line Loop(1) = {1, 2, 3, 4, 5, 6};
Plane Surface(1) = {1};
Physical Line("continuum") = {5};
Physical Line("sidewall") = {2};
Physical Line("top") = {3};
Physical Surface("domain") = {1};
''')
c = CellVariable(name='concentration', mesh=mesh, value=0.)
c.faceGrad.constrain([-0.05 * c.harmonicFaceValue], mesh.physicalFaces["sidewall"])
c.faceGrad.constrain([0.05 * c.harmonicFaceValue], mesh.physicalFaces["top"])
c.constrain(1., mesh.physicalFaces["continuum"])
dim = 1.
D = 1.
dt = 50 * dim**2 / (2 * D)
steps = 100
eq = TransientTerm() == DiffusionTerm(coeff=D)
viewer = Viewer(vars=(c, c.grad), datamin=0., datamax=1.)
for step in range(steps):
eq.solve(var=c, dt=dt)
viewer.plot()
TSVViewer(vars=(c, c.grad)).plot(filename="conc.tsv")
raw_input('Press any key...')
The calculation works as intended, however I would like to access the gradient at the physical faces I set up in GMSH. I know that I can get the gradient at faces using c.faceGrad
and the mask representing the physical faces using mesh.physicalFaces['sidewall']
. To get the gradient at the faces contained in the physical face, I would expect to use indexing like c.faceGrad[mesh.physicalFaces['sidewall']]
. However, this does not yield the desired result. Is there a way to access the FaceVariable at locations specified by physicalFaces
?
Upvotes: 0
Views: 275
Reputation: 1068
Remember that c.faceGrad
has shape (2, 22009)
so the mask operates on the second index as the shape of mesh.physicalFaces['sidewall']
is (22009,)
. So, try
c.faceGrad[:, mesh.physicalFaces['sidewall']]
and to access the correct shape use
np.array(c.faceGrad[:, mesh.physicalFaces['sidewall']]).shape
which is (2, 160)
.
Upvotes: 1