Setting different expressions in different subdomains from python code

Question

Say I want two sub domains and a different function to be interpolated in each of these. I want the functions to be solved in the same manner as one of my previous questions "interpolating vector function from python code to fenics". Basically, I want to extract the dofs of each sub domain and compute the expressions in python code, then I want to port these back into the global domain. have attached a code that is as short as I can make it where I try to solve this using submesh. The code is giving the wrong result when compared to hand calculations, and this seems to be something with the dof structure in each submesh that I'm doing wrong. Any method around this or simpler approach is greatly appreciated.

Regards
Daniel

from dolfin import *
import scitools.std as sc

# Constants
L = 1
H = 2
l = 2
m = 4
ys = 1

# Domains
mesh = RectangleMesh(0,0,L,H,l,m)
class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return (between(x[1], (0,ys)))

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return (between(x[1], (ys,H)))

bottommesh = Bottom()
topmesh = Top()
cf = CellFunction("size_t", mesh)
cf.set_all(0)
bottommesh.mark(cf,0)
topmesh.mark(cf,1)
bottom = SubMesh(mesh, bottommesh)
top = SubMesh(mesh, topmesh)

# Functionspaces
V = VectorFunctionSpace(mesh, "CG", 1)
Vs = VectorFunctionSpace(bottom, "CG", 1)
Vf = VectorFunctionSpace(top, "CG", 1)

def extract_nodes(mesh, V):
    """
    Extract x and y coordinates of mesh and
    align with dof structure
    """
    dim = V.dim()
    N = mesh.geometry().dim()
    coor = V.dofmap().tabulate_all_coordinates(mesh).reshape(dim,N)
    xdof = V.sub(0).dofmap().dofs()
    ydof = V.sub(1).dofmap().dofs()
    x = coor[:,0]
    y = coor[:,1]

    # x and y values of the vector components 
    xx, xy = x[xdof], y[xdof]
    yx, yy = x[ydof], y[ydof]

    return xx, xy, yx, yy, xdof, ydof

def func1(mesh, V):
    """
    Function to be solved in bottom layer
    """
    xx, xy, yx, yy,xdof, ydof = extract_nodes(mesh, V)
    fx = 0*xx
    fy = sc.cos(yy)

    return fx, fy, xdof, ydof

def func2(mesh, V):
    """
    Expression to be solved in top layer
    """
    xx, xy, yx, yy, xdof, ydof = extract_nodes(mesh, V)
    fx = 0*xx
    fy = sc.sin(yy)

    return fx, fy, xdof, ydof

u = Function(V)
us = Function(V)

# Get function values and dofs from each sub domain
fx1, fy1, xdof1, ydof1 = func1(bottom, Vs)
fx2, fy2, xdof2, ydof2 = func2(top, Vf)

# Set function values into vector
u.vector()[xdof1] = fx1
u.vector()[ydof1] = fy1
u.vector()[xdof2] = fx2
u.vector()[ydof2] = fy2

us.assign(interpolate(u,V))
plot(us, interactive=True)

Answer 1

Hi, you should use the map of submesh cell indices back to the parent mesh indices to extract the dofs correctly. Add something like this to you extact function

from dolfin import *
import numpy as np

def submesh_dofs(submesh, V):
    'Return dofs of V that belong to cells in submesh.'
    mesh = V.mesh()
    tdim = mesh.topology().dim()
    dofmap = V.dofmap()

    submesh_dofs = set()
    parent_cell_indices = submesh.data().array('parent_cell_indices', tdim)
    for index in range(submesh.num_cells()):
        cell = parent_cell_indices[index]
        [submesh_dofs.add(dof) for dof in dofmap.cell_dofs(cell)]

    return np.array(list(submesh_dofs))

mesh = RectangleMesh(-1, -1, 1, 1, 10, 10)

Left = AutoSubDomain(lambda x: x[0] < DOLFIN_EPS) 
Right = AutoSubDomain(lambda x: x[0] > -DOLFIN_EPS)

cell_f = CellFunction('size_t', mesh, 0)
Left.mark(cell_f, 1)
Right.mark(cell_f, 2)

left_mesh = SubMesh(mesh, cell_f, 1)
right_mesh = SubMesh(mesh, cell_f, 2)

V = FunctionSpace(mesh, 'CG', 1)
left_dofs = submesh_dofs(left_mesh, V)
right_dofs = submesh_dofs(right_mesh, V)

f = Function(V)
F = f.vector()
gdim = mesh.geometry().dim()
X = V.dofmap().tabulate_all_coordinates(mesh).reshape((-1, gdim))
X_l, X_r = X[left_dofs], X[right_dofs]

x_l, y_l = X_l[:, 0], X_l[:, 1]
x_r, y_r = X_r[:, 0], X_r[:, 1]

F[left_dofs] = x_l + 0
F[right_dofs] = x_r**2

plot(f, interactive=True) 

Comments

Hi Miro. Thank you for your response.
I'm currently trying to implement your code with V as a vector function space. I need to access each component of the velocity vectors in each sub domain separately. So in a way, after your code snippet:

X = V.dofmap().tabulate_all_coordinates(mesh).reshape((-1, gdim))
X_l, X_r = X[left_dofs], X[right_dofs]

I need the program to extract the values in a way that can look like:

xxl, xyl = X[left_dofs[xdofs]], X[left_dofs[xdofs]] # x component in sub domain 1
yxl, yyl = X[left_dofs[ydofs]], X[left_dofs[ydofs]] # y component in sub domain 1
xxr, xyr = X[right_dofs[xdofs]], X[right_dofs[xdofs]] # x component in sub domain 2
yxr, yyr = X[right_dofs[ydofs]], X[right_dofs[ydofs]] # y component in sub domain 2

And then I can define the functions in a way

fxl = xxl*xyl # x component of vector function in layer 1
fyl = yxl*yyl # y component of vector function in layer 1
fxr = xxr*xyr # x component of vector function in layer 2
fyr = yxr*yyr # y component of vector function in layer 2

And in the end set these into the function:

u.vector()[left_dofs[xdofs]] = fxl
u.vector()[left_dofs[ydofs]] = fyl
u.vector()[right_dofs[xdofs]] = fxr
u.vector()[right_dofs[ydofs]] = fyr

And interpolate this into the domain.

---

I figured it out! By using your example I have constructed a code that gets the components in each domain. The code is found below. Thank you again for excellent help Miro!

from dolfin import *
import scitools.std as sc

def submesh_dofs(mesh, submesh, functionspace):
    tdim = mesh.topology().dim()
    dofmap = V.dofmap()
    xdof = V.sub(0).dofmap()
    ydof = V.sub(1).dofmap()

    #submesh_dofs = set()
    submesh_dofx = set()
    submesh_dofy = set()
    parent_cell_indices = submesh.data().array('parent_cell_indices', tdim)
    for index in range(submesh.num_cells()):
        cell = parent_cell_indices[index]
        [submesh_dofx.add(dof) for dof in xdof.cell_dofs(cell)]
        [submesh_dofy.add(dof) for dof in ydof.cell_dofs(cell)]

    return sc.array(list(submesh_dofx)), sc.array(list(submesh_dofy))

mesh = RectangleMesh(0,0,1,1,2,2)
left = AutoSubDomain(lambda x: x[0] < 0.5 + DOLFIN_EPS)
righ = AutoSubDomain(lambda x: x[0] > 0.5 - DOLFIN_EPS)

cf = CellFunction("size_t", mesh, 0)
left.mark(cf, 0)
righ.mark(cf, 1)

lmesh = SubMesh(mesh, cf, 0)
rmesh = SubMesh(mesh, cf, 1)

V = VectorFunctionSpace(mesh, "CG", 1)

ldofx, ldofy = submesh_dofs(mesh, lmesh, V)
rdofx, rdofy = submesh_dofs(mesh, rmesh, V)

gdim = mesh.geometry().dim()
X = V.dofmap().tabulate_all_coordinates(mesh).reshape((-1, gdim)) 
x = X[:,0]
y = X[:,1]

# Vector coordinates in left layer
xxl, xyl = x[ldofx], y[ldofx]
yxl, yyl = x[ldofy], y[ldofy]

# Vector coordinates in right layer
xxr, xyr = x[rdofx], y[rdofx]
yxr, yyr = x[rdofy], y[rdofy]

u = Function(V)
u.vector()[ldofx] = xxl
u.vector()[ldofy] = 0*yyl
u.vector()[rdofx] = xxr
u.vector()[rdofy] = 0*yyr

print u.vector().array()

plot(u, interactive=True, range_max=1.0, range_min=-1.0)