This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Interpolating over time

+1 vote

Hi,
Is there any way in Fenics, that I can interpolate the solution over the time, for example I am solving dt=0.3 but I am interested in knowing the solution at time t=0.4, can I do that with out making my time step smaller?

    from dolfin import *
    import numpy

    # Create mesh and define function space
    nx = ny = 2
    mesh = UnitSquare(nx, ny)
    V = FunctionSpace(mesh, 'Lagrange', 1)

    # Define boundary conditions
    alpha = 3; beta = 1.2
    u0 = Expression('1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t',
                    alpha=alpha, beta=beta, t=0)

    class Boundary(SubDomain):  # define the Dirichlet boundary
        def inside(self, x, on_boundary):
            return on_boundary

    boundary = Boundary()
    bc = DirichletBC(V, u0, boundary)

    # Initial condition
    u_1 = interpolate(u0, V)
    #u_1 = project(u0, V)  # will not result in exact solution!

    dt = 0.3      # time step

    # Define variational problem
    u = TrialFunction(V)
    v = TestFunction(V)
    f = Constant(beta - 2 - 2*alpha)
    a = u*v*dx + dt*inner(nabla_grad(u), nabla_grad(v))*dx
    L = (u_1 + dt*f)*v*dx

    A = assemble(a)   # assemble only once, before the time stepping
    b = None          # necessary for memory saving assemeble call

    # Compute solution
    u = Function(V)   # the unknown at a new time level
    T = 1.9           # total simulation time
    t = dt
    while t <= T:
        print 'time =', t
        b = assemble(L, tensor=b)
        u0.t = t
        bc.apply(A, b)
        solve(A, u.vector(), b)

Thanks.

asked Mar 22, 2014 by Bahram FEniCS Novice (400 points)

1 Answer

+3 votes

Hi, you can perform linear interpolation between u and u0 as follows

from dolfin import *

mesh = UnitSquareMesh(25, 25)
V = FunctionSpace(mesh, 'CG', 1)

foo = Expression('sin(2*pi*x[0])*sin(2*pi*x[1])*pow(t, 2)', t=0.3)
u0 = interpolate(foo, V)

foo.t=0.6
u = interpolate(foo, V)

v = Function(V)
alpha = 1/3.
v.vector()[:] = (1 - alpha)*u0.vector().array() + alpha*u.vector().array() 
answered Mar 22, 2014 by MiroK FEniCS Expert (80,920 points)
...