Load inner boundary marks error

xiaohe · March 14, 2024, 2:59pm

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

Ekrem_Ekici · March 14, 2024, 3:09pm

Please share the full error message in the encapsulation ```

xiaohe · March 15, 2024, 2:13am

Thanks for the reply Below is my full error message

File "/home/xiaohe/fenics_code/complex_mesh_in_gmsh/mesh_labels01.py", line 28, in <module>
    bndfile.read(mvc, "bnd_marker")
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to find entity in map.
*** Reason:  Error reading MeshValueCollection.
*** Where:   This error was encountered inside HDF5File.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  2e001bd1aae8e14d758264f77382245e6eed04b0

dokken · March 15, 2024, 10:39am

First of all, you have hanging nodes in your domain. This means nodes not connected to a cell.
You can see this by calling:

print( len(np.unique(msh.cells_dict["triangle"].reshape(-1))))
print(msh.points.shape)

which returns

162
(170, 3)

which means that 8 nodes in your mesh is not in any cell.

You can find these in some of your facet markers by calling:

print(np.setdiff1d(np.unique(msh.cells_dict["line"].reshape(-1)), np.unique(msh.cells_dict["triangle"].reshape(-1))))

giving

array([ 4,  5, 50, 51, 52, 53, 54, 55])

The issue in particular is:

gmsh.model.addPhysicalGroup(1, [c1, c2], 888)

as they are being changed when calling

xiaohe · March 15, 2024, 11:22am

Thanks for the reply, I now understand my error. When I changed the code in the file that generated the mesh(remove fuse), the above error persisted.

import gmsh
import numpy as np
import meshio
gmsh.initialize()
gmsh.model.add("square_with_hole")

lc = 1.0  # mesh_size
L = 10.0   # 
R = 2.0    # 
mesh_size01 = 0.5
mesh_size02 = 0.4


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

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

ll_outer = gmsh.model.occ.addCurveLoop([l1, l2, l3, l4])
ll_inner = gmsh.model.occ.addCurveLoop([c1, c2])


s_outer = gmsh.model.occ.addPlaneSurface([ll_outer, ll_inner])
s_inner = gmsh.model.occ.addPlaneSurface([ll_inner])
gmsh.model.occ.synchronize()

# physical 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)
# generate mesh
gmsh.model.mesh.generate(2)
# save mesh
gmsh.write("./square_with_hole.msh")

# 
gmsh.finalize()

When I look from gmsh, I don’t find “hanging nodes”, so I ask first if there is a workaround for this kind of marker inner boundary problem?

dokken · March 15, 2024, 2:06pm

You need to make your physical markers on the final segments, not those defined before making any operators. i.e. once you have called

you should fetch all entities of dim 1, and tag them based on some criteria (usually I use center of mass of the entity to determine which one it should be).

xiaohe · March 15, 2024, 3:38pm

Thanks for the reply, I modified the code as suggested, but the same error occurs

import gmsh
import numpy as np
path_mesh = "D:/work files/Gmsh vs code/mesh_domain_label/square_circle01.msh"
path_of_mesh = "D:/work files/Gmsh vs code/mesh_domain_label/3d_mesh01_msh.xdmf"
path_of_boundary = "D:/work files/Gmsh vs code/mesh_domain_label/3d_mesh01_boundary.xdmf"
gmsh.initialize()
lc = 0.5
L = 1.0
R = 0.2
# 
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)

# 
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 )
# 
ll_outer = gmsh.model.occ.addCurveLoop([l1, l2, l3, l4])
ll_inner = gmsh.model.occ.addCurveLoop([c1, c2])

# 
s_outer = gmsh.model.occ.addPlaneSurface([ll_outer, ll_inner])
s_inner = gmsh.model.occ.addPlaneSurface([ll_inner])
gmsh.model.occ.synchronize()

inner01 = []
i = 0
for line in gmsh.model.getEntities(dim = 1):
    # print(line)
    com = gmsh.model.occ.getCenterOfMass(line[0], line[1])
    print(com)
    bool01 = abs(com[0] - 0.5) < 1e-6 and abs(com[1] - 1.0) > 1e-6 and abs(com[1]) > 1e-6
    if (bool01): inner01.append(line[1])

print(f"Nb of left surfaces: {len(inner01)}")
gmsh.model.addPhysicalGroup(2, [s_outer], 999)
gmsh.model.addPhysicalGroup(2, [s_inner], 777)
gmsh.model.addPhysicalGroup(1, inner01, 888)
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.2)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.1)
gmsh.model.mesh.generate(2)
gmsh.write(path_mesh)
gmsh.finalize()

xiaohe · March 15, 2024, 3:40pm

This image illustrates that my tag is right in .msh

xiaohe · March 19, 2024, 3:21am

Is there a soluation for this problem?

Topic		Replies	Views
Gmsh import error mesh	5	1411	January 30, 2021
Error in `extract_topology_markers` mesh	10	511	August 26, 2022
Gmsh, boundaries, Output,	5	1847	September 2, 2021
Marking the boundaries of a 2D mesh for a mixed function-space problem	4	647	July 17, 2022
Solving mixed poisson in circunference mesh	7	112	February 29, 2024

Load inner boundary marks error

Related Topics