How can I implement the function below correctly?

Question

from fenics import *
import time

import os
import numpy as np
%matplotlib inline

mesh = RectangleMesh(Point(0.0, 0.0), Point(1.0, 4.0), 10, 40,"crossed")

center = Point(0.5, 0.5)
radius = 0.2

class InitialCondition(Expression):
    def eval(self, values, x):
        values[0] = 0.0
        values[1] = 0.0
        values[2] = 0.0
        values[3] = sqrt((x[0]-center[0])*(x[0]-center[0]) + (x[1]-center[1])*(x[1]-center[1]))-radius
    def value_shape(self):
        return (4,)

V=VectorElement("Lagrange", mesh.ufl_cell(), 2)
P=FiniteElement("Lagrange", mesh.ufl_cell(), 1)
L=FiniteElement("Lagrange", mesh.ufl_cell(), 1)
W=MixedElement([V,P,L])

S0 = FunctionSpace(mesh,V)
S1 = FunctionSpace(mesh,P)
S2 = FunctionSpace(mesh,L)
W0 = FunctionSpace(mesh,W)

ic = InitialCondition(element=S2.ufl_element())
w = Function(W0)
w.assign(interpolate(ic,S2))

---

IndexErrorTraceback (most recent call last)
<ipython-input-31-15dd7fe51e0c> in <module>()
     31 ic = InitialCondition(element=S2.ufl_element())
     32 w = Function(W0)
---> 33 w.assign(interpolate(ic,S2))

/usr/local/lib/python2.7/dist-packages/dolfin/fem/interpolation.py in interpolate(v, V)
     64     # Compute interpolation
     65     Pv = Function(V)
---> 66     Pv.interpolate(v)
     67 
     68     return Pv

<ipython-input-31-15dd7fe51e0c> in eval(self, values, x)
     13     def eval(self, values, x):
     14         values[0] = 0.0
---> 15         values[1] = 0.0
     16         values[2] = 0.0
     17         values[3] = sqrt((x[0]-center[0])*(x[0]-center[0]) + (x[1]-center[1])*(x[1]-center[1]))-radius

IndexError: index 1 is out of bounds for axis 0 with size 1
Thanks for all.

Answer 1

Hi, consider

ic = InitialCondition(element=W0.ufl_element())
w = interpolate(ic, W0)

You would use assign for example in case you would want to have separate initial conditions (expressions associated with S0, S1, S2) that would initialize components of w.