gistfile1.txt

"""
Solve a metamorphosis problem.
"""
from firedrake import *
from numpy import array

Ncells = 50
dz = 1.0/Ncells

m = UnitIntervalMesh(Ncells)
mesh = ExtrudedMesh(m, layers=Ncells, layer_height=dz)

##Setting up the finite element spaces
V = FunctionSpace(mesh,"CG",1)

template = Expression('(tanh(x[0]/100 - 0.00875)/2 + 1.0/2.0)*(tanh(x[0]/100 - 0.0075)/2 + 1.0/2.0) + (tanh(x[0]/100 - 0.0025)/2 + 1.0/2.0)*(tanh(x[0]/100 - 0.00125)/2 + 1.0/2.0)')

target = Expression('(tanh(x[0]/100 - 0.0045)/2 + 1.0/2.0)*(tanh(x[0]/100 - 0.004)/2 + 1.0/2.0) + (tanh(x[0]/100 - 0.003)/2 + 1.0/2.0)*(tanh(x[0]/100 - 0.001)/2 + 1.0/2.0)')

#spaces are u, Helmholtz u = v, I
W = MixedFunctionSpace((V,V,V))

du,dv,dI = TestFunctions(W)
c = 1.0e-2

w = Function(W)
u,v,I = split(w)

a = (
    du*u + du.dx(0)*u.dx(0) - du*v +
    dv*v + dv.dx(0)*v.dx(0) + c*dv*u*I.dx(0)*(I.dx(1) + u*I.dx(0)) +
    (dI.dx(1) + u*dI.dx(0))*(I.dx(1) + u*I.dx(0))
    )*dx

#Boundary conditions
bcs = [DirichletBC(W[2], template, "bottom"),
       DirichletBC(W[2], target, "top")]

w0 = Function(W)

solve(a == 0, w0, bcs=bcs)