Extending a function defined on a submesh to the parent mesh

1

Hi,

I have a mesh that looks like this:

1
1905×655 33 KB

As you can see the inner circle is labeled with 3 and outside of it with 4.
Now suppose I have a function f(x,y) that is exclusively defined on the circle. I need to extend it to a function on the whole domain such that the values of the extended function is f(x,y) inside of the circle and zero outside. Inspired by this link, I came up with the following code:

from fenics import *
import numpy as np
#parent mesh 
mesh = Mesh('Mesh.xml')  
Volume = MeshFunction('size_t' , mesh , 'Mesh_physical_region.xml' ) 
bnd_mesh = MeshFunction('size_t', mesh , 'Mesh_facet_region.xml')  
###############################################################

#Submesh
Sub_Mesh=SubMesh(mesh,Volume, 3)
surface_marker = MeshFunction("size_t", Sub_Mesh, Sub_Mesh.topology().dim() - 1, 0)
volume_marker = MeshFunction("size_t", Sub_Mesh, Sub_Mesh.topology().dim())
###############################################################

#Function spaces
S = FunctionSpace(mesh,'P',1)
S_sub = FunctionSpace(Sub_Mesh,'P',1)
###############################################################

#dof and vertex maps
dofb_vb = np.array(dof_to_vertex_map(S_sub), dtype=int)
v_dof = np.array(vertex_to_dof_map(S), dtype=int)
###############################################################

#Function on the parent mesh and array for value assignments
ff = Function(S)
array = ff.vector().array()
###############################################################


#The function on the submesh
f = Expression('x[0]*x[0]',degree=0)
f_sub = project(f,S_sub)
###############################################################

#See what the function on the submesh looks like
File('sub_domain.pvd')<<f_sub

#Extend
array[v_dof[dofb_vb]] = f_sub.vector().array()
ff.vector()[:] = array
###############################################################

#See what the extension looks like
File('domain.pvd')<<ff

But the result comes out wrong. Here is the function on the submesh:

2
2905×655 34.8 KB

and this is what my method gives as the extension.

3
3905×655 61 KB

Could you please help me understand what the issue is?

Best,

2

How about using a UserExpression instead?
Suggested code showed below:

from fenics import *
import numpy as np
mesh = UnitSquareMesh(10,10)

class left(SubDomain):
    def inside(self, x, on_boundary):
        return x[0]<0.5


class expr(UserExpression):
    def __init__(self, mesh, marker, **kwargs):
        super().__init__(kwargs)
        self.mesh = mesh
        self.marker = marker

    def eval_cell(self, value, x, ufc_cell):
        cell = Cell(self.mesh, ufc_cell.index)
        if self.marker[cell] == 3:
            value[0] = x[1]
        else:
            value[0] = 0

    def value_shape(self):
        return ()


volume_marker = MeshFunction("size_t", mesh, mesh.topology().dim())
left().mark(volume_marker, 3)
V = FunctionSpace(mesh,'P',1)
f_expr = expr(mesh, volume_marker)
V = FunctionSpace(mesh, "DG", 1)
u = project(f_expr, V)

XDMFFile("u_dg.xdmf").write_checkpoint(u, "u", 0)

For another time, please post a code example using a built in mesh, as it is then easy to reproduce your issue.

1 Like
3

Thank you. It helped a lot.