Hi! I'm new to Ubuntu and Fenics and got problems with the display of my plot. Is there any special Program I have to install to get the plot to show up with command plot()? It's mostly empty or sometimes filled with glibberish. I already tried ParaViewer but that didn't really worked out either. I want to plot u[0] and u[1] for different times.
Or is my code not working at all?
# Create mesh and define function space
nx = 120
xmax = 1
mesh = IntervalMesh(nx, 0, xmax)
V = FunctionSpace(mesh, 'Lagrange', degree=1)
W=V*V
# Define const
D= 0.0005
Dc= 0.45
lamb= 30.0
alpha= 0.1
h= 10.0
k= 0.01
cM= 40.0
# Define boundary condition n=c=1 on the right side for x=1
class BoundaryRight1(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and \
(near(x[0], 1.0))
boundaryRight = BoundaryRight()
bsc = DirichletBC(V, 1.0, boundaryRight)
# Initial condition for t=0 x<1
n0=0
c0=0
u0_function= Expression(("n0", "c0"), n0=n0, c0=c0, element = W.ufl_element())
u0= Function(W)
u0.interpolate(u0_function)
#Def functions
def g (n):
return (n*(1.0+alpha**2.0))/(n**2.0+alpha**2.0)
B =(1.0+cM**2.0-2.0*h*cM)/((1.0-cM)**2.0)
def s (c):
return ((2.0*cM*(h-B)*c)/(cM**2.0+c**2.0))+B
# Vectors for saving
uk = Function(W)
uk.vector()[:]=u0.vector()
ut = Function(W)
ut.vector()[:]=ut.vector()
#saving time dependant solutions
f = File("solution/AK.pvd")
# timesteps
dt=0.3
t=0
T=50.0
while t <= T:
t = t + dt
v = TestFunction(W)
u = Function(W)
du= TrialFunction(W)
F1 = (1.0/dt * (u[0] - ut[0])-s(u[1])*(2-u[0])*u[0]+u[0])*v[0]*dx \
+ D*inner(nabla_grad (u[0]), nabla_grad (v[0]))*dx
F2 = (1.0/dt * (u[1] - ut[1]) -lamb*g(u[0]) +lamb*u[1])*v[1]*dx \
+inner(Dc*nabla_grad (u[1]), nabla_grad (v[1]))*dx
F=F1+F2
u.vector()[:]=uk.vector()
J=derivative(F,u,du)
# Solver
problem = NonlinearVariationalProblem(F, u, bsc, J)
solver = NonlinearVariationalSolver(problem)
solver.solve()
f << u
plot(u[0])
interactive()