How to extract and use physical labels from GMSH surfaces on FEniCS (2D mesh)

mesh
21

Here is what I have been able to do with the codes you provide in your site:

import pygmsh
resolution = 0.01
# Channel parameters
L = 2
H = 2

# Initialize empty geometry using the build in kernel in GMSH
geometry = pygmsh.geo.Geometry()
# Fetch model we would like to add data to
model = geometry.__enter__()

# Add points with finer resolution on left side
points = [model.add_point((0  , 0, 0), mesh_size=5*resolution),
          model.add_point((L/2, 0, 0), mesh_size=5*resolution),
          model.add_point((L/2, H, 0), mesh_size=5*resolution),
          model.add_point((0  , H, 0), mesh_size=5*resolution),
          model.add_point((L  , 0, 0), mesh_size=5*resolution),
          model.add_point((L  , H, 0), mesh_size=5*resolution)]

channel_lines_1 = [model.add_line(points[0], points[1]),
                   model.add_line(points[1], points[2]),
                   model.add_line(points[2], points[3]),
                   model.add_line(points[3], points[0])]

channel_lines_2 = [model.add_line(points[1], points[4]),
                   model.add_line(points[4], points[5]),
                   model.add_line(points[5], points[2]),
                   model.add_line(points[2], points[1])]

# Create the first line loop and plane surface for meshing
channel_loop_1 = model.add_curve_loop(channel_lines_1)
plane_surface_1 = model.add_plane_surface(channel_loop_1)

# Create the second line loop and plane surface for meshing
channel_loop_2 = model.add_curve_loop(channel_lines_2)
plane_surface_2 = model.add_plane_surface(channel_loop_2)

# Call gmsh kernel before add physical entities
model.synchronize()

model.add_physical([plane_surface_1], "Volume_1")
model.add_physical([plane_surface_2], "Volume_2")
model.add_physical([channel_lines_1[3]], "Inflow")
model.add_physical([channel_lines_2[1]], "Outflow")

geometry.generate_mesh(dim=2)
import gmsh
gmsh.write("mesh.msh")
gmsh.clear()
geometry.__exit__()

import meshio
import numpy as np
msh = meshio.read("mesh.msh")

line_cells = []
for cell in msh.cells:
    if cell.type == "triangle":
        triangle_cells = cell.data
    elif  cell.type == "line":
        if len(line_cells) == 0:
            line_cells = cell.data
        else:
            line_cells = np.vstack([line_cells, cell.data])

line_data = []
for key in msh.cell_data_dict["gmsh:physical"].keys():
    if key == "line":
        if len(line_data) == 0:
            line_data = msh.cell_data_dict["gmsh:physical"][key]
        else:
            line_data = np.vstack([line_data, msh.cell_data_dict["gmsh:physical"][key]])
    elif key == "triangle":
        triangle_data = msh.cell_data_dict["gmsh:physical"][key]

triangle_mesh = meshio.Mesh(points=msh.points[:,:2], cells={"triangle": triangle_cells})
line_mesh =meshio.Mesh(points=msh.points,
                           cells=[("line", line_cells)],
                           cell_data={"name_to_read":[line_data]})
meshio.write("mesh.xdmf", triangle_mesh)

meshio.xdmf.write("mf.xdmf", line_mesh)

from dolfin import * 
mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2) 
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

Now I got this error:

*** -------------------------------------------------------------------------
*** Error: Unable to find entity in map.
*** Reason: Error reading MeshValueCollection.
*** Where: This error was encountered inside HDF5File.cpp.
*** Process: 0

Also, I probably have not been able to mark my subdomains to call later.

22

Could you replace this code by what i refer to in:

23

You are right. Now it works. Thank you very much! Sorry for the inconvenience. At some point it did not work and I removed it. My bad.

But now I don’t know how to mark the areas to integrate and mark the individual lines to define them as boundaries. I know I can call them using:

ds_all = Measure("ds", domain=mesh, subdomain_data=mf)
dx_all = Measure("dx", domain=mesh, subdomain_data=mvc)

But ds_all(0), ds_all(1) etc. give zero. Also, I think cannot call strings.

From other topics I visit I understand that I should use something like:

geometry.add_raw_code("Physical Line(111) = {10};")
geometry.add_raw_code("Physical Line(222) = {12};")
geometry.add_raw_code("Physical Surface(333) = {16};")
geometry.add_raw_code("Physical Surface(444) = {18};")

before geometry.generate_mesh(dim=2) but now it gives another error:
‘Geometry’ object has no attribute ‘add_raw_code’

I also came across to another reply of yours where you mark cells with

cell_data = np.array([33, 55,11],dtype=np.int32)

But again, I have to ask for your help.

24

You cannot use strings, you should use integer markers instead of strings in

as shown in my other tutorials.

25

I know I have bothered you a lot today but can you give me a link? You have so many tutorials I don’t know which one to refer to - and most recent ones are for FEniCSx. Using pygmsh and meshio is also very new for me, I am really struggling. The 3D mesh tutorial in your page is confusing for me right now if that is the one you refer to.

I couldn’t find an example or tutorial on how to mark using integers.

26

See for instance: Test problem 2: Flow past a cylinder (DFG 2D-3 benchmark) — FEniCSx tutorial

Note that by opening the current xdmf files you have generated in paraview, you can visually figure out what tags has been used for each boundary.

27

I am having trouble at two points.

First problem, I checked the tags visually with Paraview and saw that it gave 3 and 4 to inlet and outlet boundaries. So when I integrate for inflow and outflow boundaries with:

print(assemble(Constant(1)*ds_all(3)))

I can see correct results. However, I cannot integrate for subdomain areas. It gives an incompatible function argument error. This is why I am trying to define my subdomains using pygmsh (or gmsh?) in the first place.

Second problem, in the tutorials you are making a kind of boolean operation.

assert(volumes == fluid[0])

But I am trying to create my geometry following a different way. So I couldn’t assert fluid to my domain as you do. Sorry that I couldn’t understand that part and I couldn’t find an example that I can follow.

Also, if we are supposed to use all gmsh for marking boundaries and subdomains, why do we even use pygmsh?

I revised my code a little bit, this way it is a little bit easier to follow now:

import pygmsh

# Initialize empty geometry using the build in kernel in GMSH
geometry = pygmsh.geo.Geometry()
# Fetch model we would like to add data to
geom = geometry.__enter__()
    
lcar = 0.1
p1 = geom.add_point([0, 0], lcar)
p2 = geom.add_point([1, 0], lcar)
p3 = geom.add_point([1, 2], lcar)
p4 = geom.add_point([0, 2], lcar)
p5 = geom.add_point([2, 0], lcar)
p6 = geom.add_point([2, 2], lcar)

l1 = geom.add_line(p1, p2)
l2 = geom.add_line(p2, p3)
l3 = geom.add_line(p3, p4)
l4 = geom.add_line(p4, p1)
l5 = geom.add_line(p2, p5)
l6 = geom.add_line(p5, p6)
l7 = geom.add_line(p6, p3)

q1 = geom.add_curve_loop([l1, l2, l3, l4])
q2 = geom.add_curve_loop([l5, l6, l7,-l2])

s1 = geom.add_plane_surface(q1)
s2 = geom.add_plane_surface(q2)

geom.synchronize()

mesh = geom.generate_mesh()

# Check if the geometry is alright
mesh.write("test.vtk")

geom.add_physical(s1, "Volume_1")
geom.add_physical(s2, "Volume_2")
geom.add_physical(l4, "Inflow")
geom.add_physical(l6, "Outflow")

rest of the code is still the same. I know you don’t approve that final part with strings but I still don’t know how to tag them with integers.

28

To add my probably incorrect recollection of events:

  • The GMSH python API was very verbose if you wanted to create simple meshes
  • pygmsh was developed as a way of writing gmsh meshes (by writing geo files)
  • Over time, pygmsh became a lightweight wrapper around the gmsh Python API

Me personally, do no longer use pygmsh, as you can see in all the DOLFINx tutorials.

As you have not supplied the whole error message, and has changed how you read in meshes, it is not a trivial job for anybody to reproduce your error. Please post your revised version of all scripts as a single post, and add the full error message.

29

I reproduced my geometry with only using gmsh now. I believe tagged the PhysicalGroups (PhysicalLines and PhysicalSurfaces) and generated .vtk and .msh files. I can see the mesh checking with gmsh and tagged boundaries with Paraview (though I am still not sure about the PhysicalSurfaces).

This time, I am having trouble generating .xdmf files.

Here is the revised version of the whole script:

import gmsh

gmsh.initialize()

geom  =  gmsh.model.geo   # Line to abbreviate geometry commands

reso = 0.1
geom.addPoint(0, 0, 0, reso, tag=1)
geom.addPoint(1, 0, 0, reso, tag=2)
geom.addPoint(1, 2, 0, reso, tag=3)
geom.addPoint(0, 2, 0, reso, tag=4)
geom.addPoint(2, 0, 0, reso, tag=5)
geom.addPoint(2, 2, 0, reso, tag=6)

geom.addLine(1, 2, tag=7)
geom.addLine(2, 3, tag=8)
geom.addLine(3, 4, tag=10)
geom.addLine(4, 1, tag=11)
geom.addLine(2, 5, tag=12)
geom.addLine(5, 6, tag=13)
geom.addLine(6, 3, tag=14)

geom.addCurveLoop([7, 8, 10, 11], tag=15)
geom.addCurveLoop([12, 13, 14,-8], tag=16)

geom.addPlaneSurface([15], tag=17)
geom.addPlaneSurface([16], tag=18)

geom.synchronize()

geom.addPhysicalGroup(2, [17], tag=19) # first region
geom.addPhysicalGroup(2, [18], tag=20) # second region
geom.addPhysicalGroup(1, [11], tag=21) # inlet
geom.addPhysicalGroup(1, [13], tag=22) # outlet

gmsh.model.geo.synchronize()
gmsh.model.mesh.generate() 
gmsh.model.mesh.set_order(2)

# check if the markings are alright
gmsh.write("unitsquare.vtk")

# generate mesh
gmsh.write("mesh.msh")

# close gmsh
gmsh.finalize()

import meshio
mesh_from_file = meshio.read("mesh.msh")

import numpy
def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    points = mesh.points[:,:2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    return out_mesh

line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)
meshio.write("facet_mesh.xdmf", line_mesh)

triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True)
meshio.write("mesh.xdmf", triangle_mesh)

ValueError Traceback (most recent call last)
in
58 return out_mesh
59
—> 60 line_mesh = create_mesh(mesh_from_file, “line”, prune_z=True)
61 meshio.write(“facet_mesh.xdmf”, line_mesh)
62

in create_mesh(mesh, cell_type, prune_z)
53 def create_mesh(mesh, cell_type, prune_z=False):
54 cells = mesh.get_cells_type(cell_type)
—> 55 cell_data = mesh.get_cell_data(“gmsh:physical”, cell_type)
56 points = mesh.points[:,:2] if prune_z else mesh.points
57 out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={“name_to_read”:[cell_data]})

~/.local/lib/python3.6/site-packages/meshio/_mesh.py in get_cell_data(self, name, cell_type)
226 def get_cell_data(self, name: str, cell_type: str):
227 return np.concatenate(
→ 228 [d for c, d in zip(self.cells, self.cell_data[name]) if c.type == cell_type]
229 )
230

ValueError: need at least one array to concatenate

30

This is due to the fact that you are generating a second order mesh

and therefore need to change

to

and

to
triangle_mesh = create_mesh(mesh_from_file, "triangle6", prune_z=True)