Define initial condition with mesh coordinates problem

Question

Hi, I am a novice here so please let me know if someone already asked this question before and had it solved. I am trying to define initial conditions using mesh coordinates. The program runs but does not generate desired result. Please have a look:

from dolfin import *
import numpy as np 

mesh = UnitSquareMesh(16,8)
V = FunctionSpace(mesh,'CG',1)

C1 = Function(V)
C1_array = C1.vector().array()

meshco = mesh.coordinates()

for i in range(mesh.num_vertices()):
    if meshco[i][0]<=0.5+DOLFIN_EPS:
        C1_array[i]=1.0
    else:
        C1_array[i]=0.0

C1.vector()[:]=C1_array

plot(C1, mode='color', elevate=0.0, range_max=1.0, range_min=0.0)
plot(mesh, interactive=True)

I checked the array in C1_array and the mesh coordinates are correct for what I wanted. The issue might arise when I imported C1_array back to C1.vector()[:]. This code results in a semi-random patterns of 0 and 1 across the unit square. I don't know how to embed the plot in here so if you could get a quick run on your computer to see what I meant.

Thank you in advance.

Answer 1

This may do what you want:

from dolfin import *
import numpy as np 

mesh = UnitSquareMesh(16,8)
V = FunctionSpace(mesh,'CG',1)

C1 = Function(V)

class InitialCondition(Expression):
    def eval_cell(self, value, x, ufc_cell):
        if x[0] <= 0.5:
            value[0] = 1.0
        else:
            value[0] = 0.0

C1.interpolate(InitialCondition())
plot(C1, interactive=True)

Comments

Awesome, christianv! It worked.
However, can you help me explain why mine didn't work? Thank you.

Answer 2

Interestingly, following code works for your case (comparing Y coordinate of a node rather than X)

for i,j in enumerate(meshco): if j[1] <= 0.5: C1_array[i]= 1.0 else: C1_array[i]= 0.

The reason might be, that Fenics iterates in X direction in mesh.coordinates subroutine, but assigns the values in Y direction (don't know if it is realy true though!). I am also looking for a similar initialization procedure, mainly because of the Expression class is expensive for larger mesh size. May be someone can throw more light on this.

Best regards,
Omkar.