# Connecting physical subdomains with common interface into a single mesh

**URL:** <https://fenicsproject.discourse.group/t/connecting-physical-subdomains-with-common-interface-into-a-single-mesh/6606>\
**Category:** mesh\
**Created:** [September 17, 2021, 8:00am UTC](https://fenicsproject.discourse.group/t/connecting-physical-subdomains-with-common-interface-into-a-single-mesh/6606 "2021-09-17T08:00:39Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Kostas\_Bardis](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/kostas_bardis/32/2968_2.png) [@Kostas\_Bardis](https://fenicsproject.discourse.group/u/Kostas_Bardis)\
**Post date:** [September 17, 2021, 8:00am UTC](https://fenicsproject.discourse.group/t/connecting-physical-subdomains-with-common-interface-into-a-single-mesh/6606/1 "2021-09-17T08:00:39Z")

</div>

Hello,

I have generated a mesh using gmsh from python API. For the shake of simplicity, in my minimal example below I created two rectangles with a common interface. It seems that fenics solves the PDE in each of the domain individually while I would like to treat the two subdomains as one and connected. How could I do that?

 ![hi](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/f/fba902481c5d93202bb141230b0bde995b939eaa.png)

```auto
from matplotlib import pyplot
import numpy as np
import math as mt
import gmsh
import meshio
from fenics import *
from dolfin import *
from sympy.polys.polytools import degree

# initialise gmsh engine
gmsh.initialize()

# assign name to geomtry
gmsh.model.add("simple_2d")
units = 1

reactangle_1 = gmsh.model.occ.addRectangle(0.0, 0.0, 0.0, units, 0.5*units, tag=1)
reactangle_2 = gmsh.model.occ.addRectangle(0.0, 0.5*units, 0.0, units, 0.5*units, tag=2)

# sycrhonise geometry with gmsh
gmsh.model.occ.synchronize()

# gather all surfaces and curves
surfaces = gmsh.model.occ.getEntities(dim=2)
curves = gmsh.model.occ.getEntities(dim=1)

gmsh.model.addPhysicalGroup(2, [1], 1)
gmsh.model.setPhysicalName(2, 1, "reactangle_1")

gmsh.model.addPhysicalGroup(2, [2], 2)
gmsh.model.setPhysicalName(2, 2, "reactangle_2")

# sycrhonise geometry with gmsh
gmsh.model.occ.synchronize()

# gnerate 2D mesh, write mesh and convert to xdmf
gmsh.model.mesh.generate(1)
gmsh.model.mesh.generate(2)
gmsh.write("mesh_2d_minimal.msh")

msh = meshio.read("mesh_2d_minimal.msh")
msh.prune_z_0()

area_mesh = meshio.Mesh(points=msh.points, cells={"triangle": msh.get_cells_type("triangle")},
 cell_data={"name_to_read": [msh.get_cell_data("gmsh:physical", "triangle")]})

meshio.write("surface.xdmf", area_mesh)

# Reading mesh data stored in .xdmf files.
mesh = Mesh()
with XDMFFile("surface.xdmf") as infile:
    infile.read(mesh)
mvc_2d = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("surface.xdmf") as infile:
    infile.read(mvc_2d, "name_to_read")
mf_2d = cpp.mesh.MeshFunctionSizet(mesh, mvc_2d)

dx = Measure("dx", domain=mesh, subdomain_data=mf_2d)

V = FunctionSpace(mesh, "Lagrange", 1)

# Define Dirichlet boundary 
def boundary(x, on_boundary):
    return on_boundary

# Define boundary condition
u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Expression("100", degree = 2)
a = inner(grad(u), grad(v))*dx
L = f*v*dx

# Compute solution
u = Function(V)
solve(a == L, u, bc)

plot_u = plot(u, title='Temperature')
plot_u.set_cmap("viridis")
pyplot.title ("Cahn - Hilliard solution ")
pyplot.colorbar(plot_u)
pyplot.savefig("hi.png", dpi = 500)
pyplot.close()

```

---

<div class="post-metadata">

**Author:** ![Emanuel\_Gebauer](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/emanuel_gebauer/32/3004_2.png) [@Emanuel\_Gebauer](https://fenicsproject.discourse.group/u/Emanuel_Gebauer)\
**Post date:** [September 17, 2021, 9:03am UTC](https://fenicsproject.discourse.group/t/connecting-physical-subdomains-with-common-interface-into-a-single-mesh/6606/2 "2021-09-17T09:03:09Z")

</div>

I think the problem is due to your boundary condition (I assume, that the parameter on\_boundary) is also true for boudnaries between subdomains). If I create 4 different boundary conditions like this, then it works.

```auto
# Define Dirichlet boundary 
def boundary_Left(x, on_boundary):
    return near(x[0],0) and on_boundary
# Define Dirichlet boundary 
def boundary_Right(x, on_boundary):
    return near(x[0],1)and on_boundary
# Define Dirichlet boundary 
def boundary_Bottom(x, on_boundary):
    return near(x[1],0)and on_boundary
# Define Dirichlet boundary 
def boundary_Top(x, on_boundary):
    return near(x[1],1)and on_boundary

# Define boundary condition
u0 = Constant(0.0)
bc0 = DirichletBC(V, u0, boundary_Left)
bc1 = DirichletBC(V, u0, boundary_Right)
bc2 = DirichletBC(V, u0, boundary_Bottom)
bc3 = DirichletBC(V, u0, boundary_Top)
bc = [bc0, bc1, bc2, bc3]

```

![CahnHillard](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/f/fa2af95ee319b51eefc07a5fd2d3f8cc2ccb45a9.png)

However I assume my quick solution is not very practical to use for you as it’s kind of dirty (at least in my opinion).  
I would define all boundaries in gmsh (name lines in your mesh) and create a 1d-MeshFunction in Fenics representing the boundaries. Applying a boundary condition is then easily done via

```auto
    bc = DirichletBC(V, u0, mf_1d, Your_Boundary_Index) 

```

Hope this helps 🙂

---

<div class="post-metadata">

**Author:** ![Kostas\_Bardis](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/kostas_bardis/32/2968_2.png) [@Kostas\_Bardis](https://fenicsproject.discourse.group/u/Kostas_Bardis)\
**Post date:** [September 17, 2021, 9:18am UTC](https://fenicsproject.discourse.group/t/connecting-physical-subdomains-with-common-interface-into-a-single-mesh/6606/3 "2021-09-17T09:18:29Z")

</div>

Dear Emanuel,

Many thanks for your reply. Unfortunately, my actual problem is not so simple and this solution would not really apply to me actual case. The issues are:

1. That I have multiple complex domains that cannot be analytically defined.
2. I would actually like to have a continuous solution of the PDE over the whole domain. The individual subdomains will be helpful in specifying different materials properties.

My actual problem is how to treat a domain that includes multiple meshes as one for the solution of the PDE.

---

<div class="post-metadata">

**Author:** ![Jesus\_Amador](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/jesus_amador/32/2894_2.png) [@Jesus\_Amador](https://fenicsproject.discourse.group/u/Jesus_Amador)\
**Post date:** [September 17, 2021, 9:18am UTC](https://fenicsproject.discourse.group/t/connecting-physical-subdomains-with-common-interface-into-a-single-mesh/6606/4 "2021-09-17T09:18:34Z")

</div>

Another possible solution in case you don’t need the line linking the two domains is to use the gmsh command:

> gmsh.model.occ.fuse

This links your two domains. The gmsh part of the code could look like this:

```auto
# initialise gmsh engine
gmsh.initialize()

# assign name to geomtry
gmsh.model.add("simple_2d")

units = 1
rectangle_1 = gmsh.model.occ.addRectangle(0.0, 0.0, 0.0, units, 0.5*units, tag=1)
rectangle_2 = gmsh.model.occ.addRectangle(0.0, 0.5*units, 0.0, units, 0.5*units, tag=2)
# link both domains and remove old rectangles
rectangles = gmsh.model.occ.fuse([(2,1)], [(2,2)], tag=3, removeObject=True, removeTool=True)

# sycrhonise geometry with gmsh
gmsh.model.occ.synchronize()

# create groups
gmsh.model.addPhysicalGroup(2, [3], 1)
gmsh.model.setPhysicalName(2, 1, "rectangles")

# gnerate 2D mesh, write mesh and convert to xdmf
gmsh.model.mesh.generate(1)
gmsh.model.mesh.generate(2)
gmsh.write("mesh_2d_minimal.msh")

# finalize gmsh engine
gmsh.finalize()

```

In case you want to keep the two domains separate you should handle the boundary conditions as indicated in the previous tip.