Hello everyone
I have created such a 2d mesh in GMSH, with a square (label. 999) on the outside and a circle (label. 777) and their interface (label. 888) on the inside. Four borders of the additional square (label: 111, 222, 333, 444).
The codethat generated the mesh is as follows:
import gmsh
import numpy as np
import meshio
gmsh.initialize()
gmsh.model.add("square_with_hole")
lc = 1.0
L = 10.0
R = 2.0
mesh_size01 = 0.5
mesh_size02 = 0.4
# Square vertices
p1 = gmsh.model.occ.addPoint(0.0, 0.0, 0.0, lc)
p2 = gmsh.model.occ.addPoint(L, 0.0, 0.0, lc)
p3 = gmsh.model.occ.addPoint(L, L, 0.0, lc)
p4 = gmsh.model.occ.addPoint(0.0, L, 0.0, lc)
# The center of the circle
pc1 = gmsh.model.occ.addPoint(L/2 + R, L/2, 0.0, lc)
pc = gmsh.model.occ.addPoint(L/2, L/2, 0.0, lc)
pc0 = gmsh.model.occ.addPoint(L/2 - R, L/2, 0.0, lc)
#
l1 = gmsh.model.occ.addLine(p1, p2)
l2 = gmsh.model.occ.addLine(p2, p3)
l3 = gmsh.model.occ.addLine(p3, p4)
l4 = gmsh.model.occ.addLine(p4, p1)
#
c1 = gmsh.model.occ.addCircleArc(pc0,pc ,pc1)
c2 = gmsh.model.occ.addCircleArc(pc0,pc ,pc1)
gmsh.model.occ.rotate(dimTags=[(1, c2)], x=L/2, y=L/2, z=0.0, ax=0, ay=0, az=1,angle=np.pi )
# cur loop
ll_outer = gmsh.model.occ.addCurveLoop([l1, l2, l3, l4])
ll_inner = gmsh.model.occ.addCurveLoop([c1, c2])
# surface
s_outer = gmsh.model.occ.addPlaneSurface([ll_outer, ll_inner])
s_inner = gmsh.model.occ.addPlaneSurface([ll_inner])
# union
gmsh.model.occ.fuse(objectDimTags= [(2, s_outer)], toolDimTags= [(2, s_inner)], removeObject= False, removeTool= False)
gmsh.model.occ.synchronize()
# phycical group
gmsh.model.addPhysicalGroup(2, [s_outer], 999)
gmsh.model.addPhysicalGroup(2, [s_inner], 777)
gmsh.model.addPhysicalGroup(1, [c1, c2], 888)
gmsh.model.addPhysicalGroup(1, [l1], 111)
gmsh.model.addPhysicalGroup(1, [l2], 222)
gmsh.model.addPhysicalGroup(1, [l3], 333)
gmsh.model.addPhysicalGroup(1, [l4], 444)
#
gmsh.model.mesh.generate(2)
# save mesh
gmsh.write("./square_with_hole.msh")
#
gmsh.finalize()
Below is an image from gmsh that illustrates the successful addition of tags
Next, I converted the .msh file to three .xdmf files with meshio:
filename = './square_with_hole.msh'
msh = meshio.read(filename)
#write mesh xdmf file
type02 = 'triangle'
type03 = 'tetra'
type01 = 'line'
meshio.write('./square_with_hole_mesh.xdmf',
meshio.Mesh(points = msh.points,
cells = {type02: msh.cells_dict[type02]}
))
#write boundary subdomains xdmf file
meshio.write('./square_with_hole_boundary_subdomains.xdmf',
meshio.Mesh(points = msh.points,
cells = {type01: msh.cells_dict[type01]},
cell_data = {"bnd_marker": [msh.cell_data_dict["gmsh:physical"][type01]]}
))
#write subdomains xdmf file
meshio.write('./square_with_hole_subdomains.xdmf',
meshio.Mesh(points = msh.points,
cells = {type02 : msh.cells_dict[type02]},
cell_data = {"dom_marker": [msh.get_cell_data("gmsh:physical",type02)]}
))
However, when I try to import the .xdmf file into fenics, I get an error like this:
# fenics load mesh with label
from fenics import *
#create mesh from msh file
mesh01 = Mesh()
#Mesh Omega
# mesh
mesh01 = Mesh()
filename_mesh = './mesh_domain_label/square_with_hole_mesh.xdmf'
with XDMFFile(filename_mesh) as mshfile:
mshfile.read(mesh01)
# Meshfunction over domains
# subdomain marks
filename_subdomain = './mesh_domain_label/square_with_hole_subdomains.xdmf'
mft = MeshFunction('size_t', mesh01, mesh01.topology().dim(), 0)
with XDMFFile(filename_subdomain) as domfile:
domfile.read(mft)
dx = Measure('dx', domain=mesh01, subdomain_data = mft)
# Meshfunction over facets
# load boundary marks
mvc = MeshValueCollection("size_t", mesh01, 1)
filename_boundary_subdomains = './mesh_domain_label/square_with_hole_boundary_subdomains.xdmf'
with XDMFFile(filename_boundary_subdomains) as bndfile:
bndfile.read(mvc, "bnd_marker")
mff = MeshFunction('size_t', mesh01, mvc)
ds = Measure('ds', domain= mesh01, subdomain_id= "everywhere", subdomain_data = mff)
dS = Measure('dS', domain= mesh01, subdomain_id= "everywhere", subdomain_data = mff)
Vh01 = FunctionSpace(mesh01, "P", 1)
print(Vh01.dim())
u01 = Function(Vh01)
u01.interpolate(Constant(1.0))
print(assemble(u01*dx(999)))
print(assemble(u01*dx))
print(assemble(u01*ds(888)))
I’m sure I’m following a standard workflow to complete the above task, but I’m really not sure what went wrong. Hope a kind person answers my question
