Connecting physical subdomains with common interface into a single mesh

mesh
1

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
hi3200×2400 215 KB 

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
2

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:

3

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.

4

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