How to read a physical object from a xdmf file

General
1

Hello,
I generate a mesh written in a xdmf file with the following script, where I add the “My Volume” and “My Obstacle” physical entities.

geometry = pygmsh.geo.Geometry()
model = geometry.__enter__()


# Add circle
circle_r = model.add_circle(c, r, mesh_size=resolution)

o_in = geometry.add_point([-L/2,0,0])
o_out = geometry.add_point([L/2,0,0])

p1 = geometry.add_point([-L/2,R,0])
p2 = geometry.add_point([-L/2-R,0,0])
p3 = geometry.add_point([-L/2,-R,0])

arc_R_in_up = model.add_circle_arc(p1,o_in,p2)
arc_R_in_down = model.add_circle_arc(p2,o_in,p3)

p4 = geometry.add_point([L/2,-R,0])
p5 = geometry.add_point([L/2+R,0,0])
p6 = geometry.add_point([L/2,R,0])

arc_R_out_down = model.add_circle_arc(p4,o_out,p5)
arc_R_out_up = model.add_circle_arc(p5,o_out,p6)

channel_lines = [arc_R_in_up, arc_R_in_down, model.add_line(p3, p4), arc_R_out_down, arc_R_out_up, model.add_line(p6,p1)]

plane_surface = model.add_plane_surface(channel_loop, holes=[circle_r.curve_loop])

model.synchronize()

model.add_physical([plane_surface], "My Volume")
model.add_physical(circle_r.curve_loop.curves, "My Obstacle")

I then read the xdmf file with

mesh=Mesh()
with XDMFFile("../mesh/triangle_mesh.xdmf") as infile:
    infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("../mesh/line_mesh.xdmf") as infile:
    infile.read(mvc, "name_to_read")



# Define function spaces
V = VectorFunctionSpace(mesh, 'P', 2)
Q = FunctionSpace(mesh, 'P', 1)

# Define boundaries
obstacle = 'on_boundary && (x[0]*x[0] + x[1]*x[1] < (0.5*0.5))'

I would like the last line, to define the obstacle. I would like to define obstacle by reading the pysical object above with tag “My Obstacle” directly. How can I do that?
I looked at similar posts about this, but they are unclear

Thank you

2

First, maybe you can open your .xdmf file in gmsh, choose physical group to check the physcial name is correct or not. below is one of my code

gmsh.model.addPhysicalGroup(3, volumes[0],  31)
gmsh.model.addPhysicalGroup(3, volumes[1],  32)

Here 31 and 32 is physical name, it seems use word like ‘‘My Volume’’ can not distinguish by gmsh.

After that, you can define obstcale by

obstacle= mesh.locate_entities_boundary(
       mesh, fdim, lambda x: x[0]*x[0] + x[1]*x[1]  <=  (0.5*0.5))

Here fdim is dimension of your mesh -1, for your 2d case, it should be 1.

Have a try, good lucky.

3

Thank you, but that does not answer my question. Please read the OP.

4

I would say that @Yanjun does answer your question, as DOLFINx cannot read string tags.

You should use the GMSH Python API directly if you want to have control of the integer given to a physical marker.

Please also note that your mesh generation script is not complete, as it is missing several definitions, including channel_loop, R, L and resolution.

For further questions, please make an effort to make the code reproducible, by including all definitions, all imports and what versions of legacy DOLFIN, pygmsh, how you convert your msh file to xdmf etc is done.

5

Hello,
As far as I can tell, I am not using DolfinX but legacy Dolfin. Here is a revised version of my script:

import numpy
import meshio
import gmsh
import pygmsh
resolution = 0.08
R = 1.0
r = 0.25
c = [0.0, 0.0, 0.0]


geometry = pygmsh.geo.Geometry()
model = geometry.__enter__()


circle_r = model.add_circle([0,-0.2,0], r, mesh_size=resolution)
circle_R = model.add_circle(c, R, mesh_size=resolution)


plane_surface = model.add_plane_surface(     circle_R.curve_loop, holes=[circle_r.curve_loop])

model.synchronize()
model.add_physical([plane_surface], "Volume")
model.add_physical(circle_r.curve_loop.curves, "Obstacle")

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

mesh_from_file = meshio.read("membrane_mesh.msh")
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("line_mesh.xdmf", line_mesh)
triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True)
meshio.write("triangle_mesh.xdmf", triangle_mesh)

The script above runs. It is still unclear to me from the answers to this post how to select and use an object from this script by assigning a given tag to the object. I use the .msh file generated from the script above in the following python code



from __future__ import print_function
from fenics import *
from mshr import *
import numpy as np
import meshio
import ufl as ufl


input_directory = "/home/fenics/shared/mesh/membrane_mesh"

#create mesh with new method
mesh=Mesh()
with XDMFFile(input_directory + "/triangle_mesh.xdmf") as infile:
    infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile(input_directory + "/line_mesh.xdmf") as infile:
    infile.read(mvc, "name_to_read")

# Define function spaces
V = VectorFunctionSpace(mesh, 'P', 2)
Q = FunctionSpace(mesh, 'P', 1)

cylinder = 'on_boundary && (x[0]*x[0] + x[1]*x[1] < (0.5*0.5))'

Can you show me a working example of how I can define cylinder in the second script by assigning a tag to a physical object (circle_r ?) in the first script rather than using the cumbersome notation that the radius must be smaller than 0.5 and that I must be on the boundary?

Thank you

6

Hi Mekong, I think there maybe has better way to solve this problem, but below code works for me:

gdim = 2
domain, cell_markers, facet_markers = gmshio.read_from_msh("membrane_mesh.msh", MPI.COMM_WORLD, 0, gdim)
fdim =gdim -1
cylinder = mesh.locate_entities_boundary( domain, fdim , lambda x:  x[0]**2 + x[1]**2 <= 0.5**2 )

print(cylinder)


Info    : Reading 'membrane_mesh.msh'...
Info    : 17 entities
Info    : 664 nodes
Info    : 1247 elements
Info    : Done reading 'membrane_mesh.msh'
[ 544  545  546  577  579  628  656  709  742  794  867  870  942  963
 1034 1037 1110 1176 1259 1262 1343]

I think above 21 poins forms the curve you want.

Screenshot from 2024-04-07 01-12-19
Screenshot from 2024-04-07 01-12-19743×605 52.5 KB

Above code you can find in Deflection of a membrane .

7

You are using pygmsh, which doesn’t give you an option to choose the integer associated with a physical group. If you open line_mesh.xdmf in Paraview you see that the inner circle is marked with the integer 2.
Thus you can integrate over that part of the boundary as shown below:

import numpy as np
from fenics import *


input_directory = "."

# create mesh with new method
mesh = Mesh()
with XDMFFile(input_directory + "/triangle_mesh.xdmf") as infile:
    infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 1)
with XDMFFile(input_directory + "/line_mesh.xdmf") as infile:
    infile.read(mvc, "name_to_read")


mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
ds_inner = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=2)

inner_circumference = assemble(1*ds_inner)

print(inner_circumference/(2*np.pi))

print((assemble(1*ds(domain=mesh))-inner_circumference)/(2*np.pi))

returning the inner and outer radius as expected

0.2490685406113397
0.9997493049308164
1 Like
8

That works like a charm, thank you !