Connecting physical subdomains with common interface into a single mesh

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?

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()

1 Like

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.

# 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

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

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

Hope this helps :slight_smile:

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.

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:

# 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.

3 Likes