Hello Everyone,
I have the following equation:
Is this the correct way of implementing it in fenics?
ui = Expression('-40.*sin(pi*x[0])*sin(pi*x[1])*exp(-1000.*(pow(x[0]-0.5,2)+pow(x[1]-0.5,2)-c*t))'), c=c, t=0.0)
dt = 0.005
u0 = interpolate(ui, V)
t = dt
ui.t = t
u1 = interpolate(ui, V)
a = 1/dt**2*inner(u, v)*dx \
+ c**2*beta*inner(nabla_grad(v), nabla_grad(u))*dx
L = 2.0/dt**2*inner(v,u1)*dx \
- 1/dt**2*inner(v, u0)*dx \
- (1.0 - 2.0*beta)*c**2*inner(nabla_grad(v), nabla_grad(u1))*dx \
- beta*c**2*inner(nabla_grad(v), nabla_grad(u0))*dx
S = SystemAssembler(a, L, bc) # Assemble linear system
A = PETScMatrix()
b = PETScVector()
S.assemble(A)
u = Function(V) # The unknown u at the new time level
solver = LUSolver()
solver.set_operator(A)
while t <= T:
S.assemble(b)
solver.solve(u.vector(), b)
u0.assign(u1)
u1.assign(u)
file << (u1,t)
t += dt