How to define a Function in a submesh based on function values present in an adjoining submesh?

1

Hello,

I have two submeshes - 1 and 2 as shown below. They share a common boundary between them.

Domain_Mesh - Copy

I have a function u defined on the function space, Vusb1 in the top submesh. I have obtained this function by solving a PDE in this submesh.

u = Function(Vsub1)  

# ... I have then solved for this u on Vsub1 and obtained a solution

I am now trying to define a function u2 in the bottom submesh that takes the values of u at the shared boundary and zero everywhere else. Something like this:

u2 = Function(Vsub2)
values of u2(at the shared boundary) = values of u at the shared boundary (This u is defined in the top submesh)

I know the the co-ordinates of this shared boundary and I have also marked this boundary as 1 in both the submeshes.

I’m struggling to implement this.

Thanks in advance for the help offered !

2

Hi, Please add a minimal working code example, as described in Read before posting: How do I get my question answered?
Without such an example, it requires more effort from the community to try to help you with your problem, as one has to create example code and hopefully interpret your question correctly. Especially since SubMesh is a relatively new feature, a code example will help others that might face similar problems.

3

Dear @dokken

Thanks a lot for the reply and guidance. I have built a minimal working code that describes the issue for which I have requested help. The code is attached below:

I have a function u defined on the function space, Vsub in the top submesh. I have obtained this function by solving a PDE in this submesh.

from dolfin import *
from fenics import *
from ufl import nabla_div
from ufl import sinh

r = Constant(25.2)
g = -100
C = Constant(19.34)
kappa = 0.04

# create mesh
mesh = RectangleMesh(Point(0,0), Point(100e-6, 200e-6), 100,100)

regions = MeshFunction('size_t', mesh, mesh.topology().dim())
boundaries = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)

class A(SubDomain):
    def inside(self, x, on_boundary):
        return x[1] < 100e-6            #Bottom half of domain

a = A()

class B(SubDomain):
    def inside(self, x, on_boundary):
        return x[1] >= 100e-6           #Top half of domain
b = B()

class TopBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 200e-6)
top = TopBoundary()

class MiddleBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near (x[1], 100e-6)
middle = MiddleBoundary()

regions.set_all(0)                     
b.mark(regions, 1)                      #Top region is marked 1; bottom remains marked as 0

boundaries.set_all(0)
middle.mark(boundaries, 2)            #Middle boundary in the domain
top.mark(boundaries, 1)                 #Top most boudnary in the domain


# Working on the top submesh - To solve a PDE

# define a submesh, composed of region 1 - top region 
new_region = MeshFunction('size_t', mesh, mesh.topology().dim())
new_region.set_all(0)
new_region.array()[regions.array() == 1] = 1
submesh = SubMesh(mesh, new_region, 1)

# define new meshfunctions on this submesh 
submesh_regions = MeshFunction('size_t', submesh, submesh.topology().dim())
submesh_regions.set_all(0)
submesh_boundaries = MeshFunction('size_t', submesh, submesh.topology().dim() - 1)
submesh_boundaries.set_all(0)

middle.mark(submesh_boundaries, 2)
top.mark(submesh_boundaries, 1)

ds = Measure("ds")[submesh_boundaries]

Vsub = FunctionSpace(submesh, 'CG', 1)
u = TrialFunction(Vsub)
v = TestFunction(Vsub)
u = Function(Vsub)
dx = Measure("dx")[submesh_regions]


bcs = []
F = kappa*inner(nabla_grad(u), nabla_grad(v))*dx + g*v*ds(1) + r*sinh(C*u)*v*ds(2)

J = derivative(F, u)
problem = NonlinearVariationalProblem(F, u, bcs, J)
solver  = NonlinearVariationalSolver(problem)
solver.solve()

#Plotting solution u obtained on Vsub - defined on top submesh 
import matplotlib.pyplot as plt
from matplotlib import cm

p = plot(u, cmap = cm.inferno)
plt.colorbar(p)
plt.show()

# Working on the bottom submesh 

# define a submesh, composed of region 0 - bottom region 
new_region2 = MeshFunction('size_t', mesh, mesh.topology().dim())
new_region2.set_all(0)
new_region2.array()[regions.array() == 0] = 1
submesh2 = SubMesh(mesh, new_region2, 1)

# define new meshfunctions on this submesh 
submesh_regions2 = MeshFunction('size_t', submesh2, submesh2.topology().dim())
submesh_regions2.set_all(0)
submesh_boundaries2 = MeshFunction('size_t', submesh2, submesh2.topology().dim() - 1)
submesh_boundaries2.set_all(0)

middle.mark(submesh_boundaries2, 2)

ds2 = Measure("ds")[submesh_boundaries2]

Vsub2 = FunctionSpace(submesh2, 'CG', 1)
u2 = TrialFunction(Vsub2)
v2 = TestFunction(Vsub2)
dx2 = Measure("dx")[submesh_regions2]

u2 = Function(Vsub2)

Following this, I need to define a function u2 in the bottom submesh that takes the values of u at the shared boundary and zero everywhere else. Something like this:

values of u2(at the shared boundary) = values of u at the shared boundary (This u is defined in the top submesh)

Since, both the boundaries are marked and share the same co-ordinates, I thought of extracting the dofs and modifying the function u2 as follows:

from numpy import where
top_dofs = where(submesh_boundaries.array()==2)
bottom_dofs = where(submesh_boundaries2.array()==2)

u2.vector()[bottom_dofs[0]] = u.vector()[top_dofs[0]]

It shows the following error.

u2.vector()[bottom_dofs[0]] = u.vector()[top_dofs[0]]
IndexError: Vector index out of range

I’m doing some fundamental mistake in assigning the function values based on the dofs, and not able to update the function u2 properly.

Please help me resolving this issue. Thanks in advance !

4

Dear @Sabarish_Gopi,

As this is not my field of expertise, my solution to this problem is highly specialized to CG-1 spaces.
Maybe @cdaversin has a better solution than this.

Please note that your example is far from a minimal example,as you have included solving a non-linear PDE. Let the following code be a guide on how to produce a Minimal code example in the future.

from  dolfin  import *
import numpy as np
n = 10

mesh = UnitSquareMesh(n, n)

# Divide a mesh into two parts, x[0]<0.7 and x[0]>=0.7
marker = MeshFunction("size_t",mesh ,mesh.topology().dim(),0)
for c in  cells(mesh):
    marker[c] = 0.7 < c.midpoint().x()
marker_array = marker.array()
submesh0 = MeshView.create(marker ,0)
submesh1 = MeshView.create(marker ,1)

# Initialize  function  spaces  and  basis  functions
V0 = FunctionSpace(submesh0 , "CG", 1)
V1 = FunctionSpace(submesh1 , "CG", 1)
W = MixedFunctionSpace(V0,V1)

# Create empty function, and interpolate u=x[1] onto part 0
u = Function(W)
u.sub(0).interpolate(Expression("x[1]", degree=1))

# Create mapping from submesh to parent mesh
mv0 = submesh0.topology().mapping()[mesh.id()]
vertex_map0 = np.array(mv0.vertex_map(),dtype=np.int32)
mv1 = submesh1.topology().mapping()[mesh.id()]
vertex_map1 = np.array(mv1.vertex_map(), dtype=np.int32)
# Find which vertices are in common for the two submeshes
common_vertices = np.flatnonzero(np.isin(vertex_map0, vertex_map1))

dofs0 = vertex_to_dof_map(V0)
dofs1 = vertex_to_dof_map(V1)
# Map dofs from submesh 0 to submesh 1 through their common vertex.
# To do this with other function-spaces than CG-1, mappings
# between facets has to be added
for vertex in common_vertices:
    vertex_0 = vertex_map0[vertex]
    dof0 = dofs0[vertex]
    vertex_1 = np.flatnonzero(vertex_map1 == vertex_0)[0]
    dof1 = dofs1[vertex_1]
    u.sub(1).vector().vec().setValueLocal(dof1,  u.sub(0).vector()[dof0])

File("u0.pvd") << u.sub(0)
File("u1.pvd") << u.sub(1)
1 Like
5

Dear @dokken,

Thanks a lot for spending time to help me. I now understand how to provide a minimal working code.

In the code you have shared, it shows the following error for MeshView and MixedFunctionSpace.

NameError: name 'MeshView' is not defined

Any leads on what I’m missing ? Thanks again !

6

To use the meshview functionality was added as part of the mixed dimensional FE implementation in FEniCS (see arxiv paper).

This is added to the development version of dolfin, and can for instance be found with docker:

docker run -ti -v $(pwd):/home/fenics/shared -w /home/fenics/shared/ --rm quay.io/fenicsproject/dev:latest
1 Like
7

Dear @dokken

Is there a way around the usage of MeshView and MixedFunctionSpace provided just by mapping the dofs from submesh0 to submesh1 through the common vertices.
I have split the domain into two submeshes, also I have a function u defined on one of the submeshes.

If not, how to get started with the docker ? I’m very new to it.

Thank you very much.

8

As i have not used submeshes I do not know. There is probably possible to do something similar to what I’ve done above. If you can make your problem similar to the code I posted I might Give it a go. Your current problem is too big for me to take time in reducing it.
There is a long tutorial on docker here: https://fenics.readthedocs.io/projects/containers/en/latest/introduction.html#installing-docker

1 Like
9

Hi! I have a question about modeling: if my model consists of two regions with a common boundary and I want to use FEniCS to calculate heat transfer or mechanical problems, when drawing the mesh using Gmsh, should the common boundary be drawn as a single edge or as two overlapping edges (such as the interface between a solid and a liquid)? What are the differences between these two methods? If it should be drawn as two overlapping edges, then how can I match these two edges in FEniCS? I hope to get an answer, thank you!