Hi,
I am currently trying to create a heat equation with a non constant diffusion coefficient.
The diffusion coefficient is supposed to look kind of like this: (in 1st to 5th cell its one value for all rest its another one)
Testimg|596x500
I want the coefficient to be based on the cell (I have a list of which cells have what coefficient), but I have no idea how to set the coefficient.
Ive first tried to just set it by the vertex, but this has the obvious issue that the coefficient it calculates with on the border is wrong. Is there a way to find out on which cell its currently working to use that to set the coefficient?
The code ive given here is just a test code, the real code is more complex and in 3d, so I cant just use a cheap workaround.
Code
from dolfin import *
import matplotlib.pyplot as plt
import numpy.linalg as LA
T = 5 # final time
num_steps = 10 # number of time steps
dt = T / num_steps # time step size
mesh = UnitIntervalMesh(10)
V = FunctionSpace(mesh, 'CG', 1)
u = TrialFunction(V)
v = TestFunction(V)
class Source(UserExpression):
def eval(self, values, x):
if(x<=0.5):
values[0] = 1
print(x,values[0],n,"Out1")
else:
values[0] = 0
print(x,values[0],n,"Out5")
u_0 = Expression('10*(x[0]<0.5)',element = V.ufl_element())
u_n = interpolate(u_0, V)
plot(u_n)
plt.show()
f = Source(degree=1)
F = u*v*dx + f*dt*dot(grad(u), grad(v))*dx - (u_n)*v*dx
a, L = lhs(F), rhs(F)
# Time-stepping
u = Function(V)
t = 0
for n in range(num_steps):
print(n)
# Update current time
t += dt
# Compute solution
solve(a == L, u)#, bc
# Plot solution
plot(u)
# Update previous solution
u_n.assign(u)
# Hold plot
plt.show()