Workaround for submesh or how to extract a function on a subdomain?

mesh
1

Hello !
I am struggling to sort my problem out. Maybe someone here could help!
I created a toy problem to explain my problem. Say I have a domain \Omega and a subdomain \Omega_1 as drawn here:

domain
domain600×600 5.5 KB

I solve a first equation E_1(u_\text{ref}) = 0 which defines a function u_\text{ref} on the whole domain \Omega.
Then, I have a second equation E_2(u) = u_\text{ref} but only in the subdomain \Omega_1, and defined thanks to u_\text{ref}.

I have defined a code, which works well, thanks to the Submesh function:

from dolfin import *
from mshr import Rectangle, generate_mesh
import matplotlib.pyplot as plt

# Create classes for defining parts of the boundaries and the interior
# of the domain
class Small_domain(SubDomain):
    def inside(self, x, on_boundary):
        return between(x[0], (0.0, 4.0))

class Small_boundary(SubDomain):
    def inside(self, x, on_boundary):
        return (on_boundary and between(x[0], (0.0, 4.0)))

# Initialize sub-domain instances
small_domain = Small_domain()
small_boundary = Small_boundary()

geometry = Rectangle(Point(0,0), Point(6,4))
mesh = generate_mesh(geometry, 50)
small_mesh = SubMesh(mesh, small_domain)

# Initialize mesh function for boundary domains
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
boundaries.set_all(0)
small_boundary.mark(boundaries, 1)

# Define function space and basis functions
V = FunctionSpace(mesh, "CG", 2)
u_ = TrialFunction(V)
v = TestFunction(V)
uRef = Function(V)

# Define Dirichlet cdts
bcs_big = [DirichletBC(V, 0.0, boundaries, 0), DirichletBC(V, 5.0, boundaries, 1)]

# Define variational form
F_big = inner(grad(u_), grad(v))*dx + inner(u_, v)*dx
a_big, L_big = lhs(F_big), rhs(F_big)
solve(a_big == L_big, uRef, bcs_big)

fig = plt.figure()
cntr = plot(uRef, title="uRef")
plt.colorbar(cntr)
plt.show()


#Extract function on subdomain
V_small = FunctionSpace(small_mesh, "CG", 2)
small_uRef = Function(V_small)
small_uRef = project(uRef, V_small)

#define problem on subdomain
boundaries = MeshFunction("size_t", small_mesh, mesh.topology().dim()-1, 0)
boundaries.set_all(0)
small_boundary.mark(boundaries, 1)

u_ = TrialFunction(V_small)
v = TestFunction(V_small)
u = Function(V_small)

bcs_small = [DirichletBC(V_small, 5.0, boundaries, 1)]
F_small = inner(grad(u_), grad(v))*dx + inner(small_uRef, v)*dx

# Separate left and right hand sides of equation
a_small, L_small = lhs(F_small), rhs(F_small)

# Solve problem
solve(a_small == L_small, u, bcs_small)

# Plot solution
fig = plt.figure()
cntr = plot(u, title="u")
plt.colorbar(cntr)
plt.show()

(The implemented equations are not necessarily meaningful, it’s not the point.)
The thing is: this works fine when it’s launch in serial, but my real problem should be solved using MPI (because it’s really big). And SubMesh does not work with MPI.

Do you know any way I could do this?
Thank you !

2

Have you tried to use MeshView?
There are several demos showing the applicability of meshview in dolfin: Bitbucket
As well as @cdaversin’s paper: https://arxiv.org/pdf/1911.01166.pdf

1 Like
3

It seems like it’s working, once I activate the extrapolation option. If I have a problem with my real problem I’ll come back here.
Thanks for the references!

I’m posting the code I ran, which works with MPI:

from dolfin import *
from mshr import Rectangle, generate_mesh
import matplotlib.pyplot as plt

# Create classes for defining parts of the boundaries and the interior
# of the domain
class Small_domain(SubDomain):
    def inside(self, x, on_boundary):
        return between(x[0], (0.0, 4.0))

class Small_boundary(SubDomain):
    def inside(self, x, on_boundary):
        return (on_boundary and between(x[0], (0.0, 4.0)))

# Initialize sub-domain instances
small_domain = Small_domain()
small_boundary = Small_boundary()

geometry = Rectangle(Point(0,0), Point(6,4))
mesh = generate_mesh(geometry, 50)

small_domain_function = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
small_domain.mark(small_domain_function, 1)
small_mesh = MeshView.create(small_domain_function, 1)


# Initialize mesh function for boundary domains
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
boundaries.set_all(0)
small_boundary.mark(boundaries, 1)

# Define function space and basis functions
V = FunctionSpace(mesh, "CG", 2)
u_ = TrialFunction(V)
v = TestFunction(V)
uRef = Function(V)

# Define Dirichlet cdts
bcs_big = [DirichletBC(V, 0.0, boundaries, 0), DirichletBC(V, 5.0, boundaries, 1)]

# Define variational form
F_big = inner(grad(u_), grad(v))*dx + inner(u_, v)*dx
a_big, L_big = lhs(F_big), rhs(F_big)
solve(a_big == L_big, uRef, bcs_big)

#Extract function on subdomain
uRef.set_allow_extrapolation(True)
V_small = FunctionSpace(small_mesh, "CG", 2)
small_uRef = Function(V_small)
small_uRef = project(uRef, V_small)

#define problem on subdomain
boundaries = MeshFunction("size_t", small_mesh, mesh.topology().dim()-1, 0)
boundaries.set_all(0)
small_boundary.mark(boundaries, 1)

u_ = TrialFunction(V_small)
v = TestFunction(V_small)
u = Function(V_small)

bcs_small = [DirichletBC(V_small, 5.0, boundaries, 1)]
F_small = inner(grad(u_), grad(v))*dx + inner(small_uRef, v)*dx

# Separate left and right hand sides of equation
a_small, L_small = lhs(F_small), rhs(F_small)

# Solve problem
solve(a_small == L_small, u, bcs_small)

# Plot solutions
out_big = XDMFFile(MPI.comm_world, 'big_domain_output.xdmf')
out_big.write(uRef)
out_small = XDMFFile(MPI.comm_world, 'small_domain_output.xdmf')
out_small.write(u)
4

Actually… After further tests, I noticed it is only working with a limited number of processors. Namely, It stops working once I run it using MPI on 6 or 8 processors.
I think the MeshView function leaves some processors without any node, which makes the rest of the code bug. Is there a way to remap this ?

5

It seems that MeshView needs the parent mesh distributed in a way that each process have a part of the child and the parent\child meshes. When this is not the case, the child mesh created by MeshView is empty in some processes, which breaks some dolfin functionalities. I had this issue before and tracked the error (presumably) to dolfin/mesh/DistributedMeshTools.cpp (but I don’t know how to fix it).

1 Like