Extracting node points from subdomain

Question

I'm trying to compute a time dependent expression object in a domain that is divided into two subdomains. After running the program, I would like to extract the values of the solution at the node points in the first subdomain. I have a code under that tries to implement this, but I end up with an array of size 3 no matter how fine my mesh is.

The code is given below. Any help or advice is greatly appreciated!



from dolfin import *

import math as mt

# size of Mesh
L = 1
h = 1
l = 39
m = 39

# Mesh, functionspace, subdomains
mesh = RectangleMesh(0,0,L,h,l,m)
cf = CellFunction("size_t", mesh, 0)
V = FunctionSpace(mesh, "Lagrange", 1)
o2 = AutoSubDomain(lambda x: x[0] > 0.5)
o2.mark(cf,1)

# Constants
omega = 10
c1 = 1
lk = mt.pi/(2*h)
kk = mt.sqrt(omega**2/c1**2 - lk**2)

# Time parameters
dt = 0.1
T = 10
t = 0

# Exact solution
u_e = Expression("sin(%s*t - %s*x[0])*cos(%s*x[1])"% (omega, kk, lk), t=t)

# Main loop and plot
u_p = Function(V)
while t <= 10:
    begin("Computing at time t=%g" % t)
    u_e.t = t
    u_p.assign(interpolate(u_e, V))
    plot(u_p, rescale=False)
    end()

    t += dt

# Extract solution array from domain 1
for cell in cells(mesh):
    i = cell.index()
    c = cf[i]
    if c == 1:
        dofs = V.dofmap().cell_dofs(i)
        vals = u_p.vector()[dofs]

print vals.array()

MiroK · Answer 1 · Jan 14, 2014

Hi, the way you wrote the code, you will end up with vals pointing to the values at cell_dofs of the last cell in the cells(mesh) iterator. The length of vals is 3 since you are using cG1 elements on triangles. You are not adding vals to some longer array, but always pointing it to something else! As a suggestion, in the code below the values of function in the domain of interest are used to build a dictionary and in the rest of the domain the values are set to meaningless values so you can clearly distinguish the two in the plot.

u_pv = u_p.vector()                                                              
vals = {}                                                                        
for cell in cells(mesh):                                                         
    i = cell.index()                                                             
    c = cf[i]                                                                    
    dofs = V.dofmap().cell_dofs(i)                                               
    if c == 0: # set domain0 values to something unphysical                      
      u_pv[dofs] = -1.5                                                          
    else:                                                                        
      vals[i] = u_pv.array()[dofs]                                               

print vals                                                                       
plot(u_p, interactive=True)

Extracting node points from subdomain

1 Answer