myprojects-hg-reborn

Fork

#!/usr/bin/env python


from fipy import *


cellSize = 0.05
radius = 1.

transient=0

mesh = GmshImporter2D('''
                       cellSize = %(cellSize)g;
                       radius = %(radius)g;
                       Point(1) = {0, 0, 0, cellSize};
                       Point(2) = {-radius, 0, 0, cellSize};
                       Point(3) = {0, radius, 0, cellSize};
                       Point(4) = {radius, 0, 0, cellSize};
                       Point(5) = {0, -radius, 0, cellSize};
                       Circle(6) = {2, 1, 3};
                       Circle(7) = {3, 1, 4};
                       Circle(8) = {4, 1, 5};
                       Circle(9) = {5, 1, 2};
                       Line Loop(10) = {6, 7, 8, 9};
                       Plane Surface(11) = {10};
                       ''' % locals())



phi = CellVariable(name = "solution variable",mesh = mesh,value = 0.)

viewer = None
if __name__ == '__main__':
    try:
        viewer = Viewer(vars=phi, datamin=-1, datamax=1.)
        viewer.plotMesh()
        raw_input("Irregular circular mesh. Press <return> to proceed...")
    except:
        print "Unable to create a viewer for an irregular mesh (try Gist2DViewer, Matplotlib2DViewer, or MayaviViewer)"

D = 1.

X, Y = mesh.getFaceCenters()

BCs = (FixedValue(faces=mesh.getExteriorFaces(), value=X),)


if (transient==1):

    eq = TransientTerm() == DiffusionTerm(coeff=D)


    timeStepDuration = 10 * 0.9 * cellSize**2 / (2 * D)
    steps = 10
    for step in range(steps):
        eq.solve(var=phi,boundaryConditions=BCs,dt=timeStepDuration)
        if viewer is not None:
            viewer.plot()


    TSVViewer(vars=(phi, phi.getGrad())).plot(filename="myTSV.tsv")

# SS solution


DiffusionTerm(coeff=D).solve(var=phi,boundaryConditions=BCs)


print "phi.allclose(X, atol = 0.02), ", phi.allclose(X, atol = 0.02)

if viewer is not None:
    viewer.plot()
    raw_input("Steady-state diffusion. Press <return> to proceed...")




print "So far so good"