Cannot assign and extract a label to a single line in a mesh

mesh
,
1

Hello,
Given this mesh

Screenshot 2024-06-22 at 23.14.07
Screenshot 2024-06-22 at 23.14.072112×532 22.2 KB

I can successfully label in it the circle and the rectangle, save the mesh to file, and read the circle and rectangle. I prepared a minimal working example to illustrate the issue:

First, generate_mesh.py generates the mesh and saves it into line_mesh.xdmf and triangle_mesh.xdmf

#generate_mesh.py:

import numpy
import meshio
import gmsh
import pygmsh
import argparse

parser = argparse.ArgumentParser()
parser.add_argument("resolution")
args = parser.parse_args()

resolution = (float)(args.resolution)

# Channel parameters
L = 2.2
h = 0.41
r = 0.05
c_r = [0.2, 0.2, 0]

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

circle_r = model.add_circle(c_r, r, mesh_size=resolution/2)
rectangle_Lh = model.add_rectangle(0, L, 0, h, 0, mesh_size=resolution)

points = [model.add_point((0, 0, 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),
          model.add_point((0, 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_loop_1 = model.add_curve_loop(channel_lines_1)

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

model.synchronize()
model.add_physical([plane_surface], "Volume")
model.add_physical(circle_r.curve_loop.curves, "Circle")
model.add_physical(rectangle_Lh.curve_loop.curves, "Rectangle")
# model.add_physical([channel_lines_1[3]], "Inflow")


geometry.generate_mesh(64)
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)

Then mwe.py reads line_mesh.xdmf and triangle_mesh.xdmf and finds the circle and rectangle.

#mwe.py:

import math

from fenics import *
from mshr import *
import numpy as np
import meshio
import ufl as ufl
import argparse

parser = argparse.ArgumentParser()
parser.add_argument("input_directory")
args = parser.parse_args()

L = 2.2
h = 0.41
r = 0.05
c_r = [0.2, 0.2]

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

mesh_coordinates = mesh.coordinates()

#read an object with label subdomain_id from xdmf file and assign to it the ds `ds_inner`
mf = dolfin.cpp.mesh.MeshFunctionSizet(mesh, mvc)
ds_circle = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=2)
ds_rectangle = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=3)
circle_perimeter = assemble(1*ds_circle)
rectangle_perimeter = assemble(1*ds_rectangle)
print("Circle perimeter = ", circle_perimeter)
print("Rectangle perimeter = ", rectangle_perimeter)

The circle and rectangle seem to be read correctly by mwe.py, because the values of their length are correct:

$ python3 generate_mesh.py 0.1
$ python3 mwe.py [folder where generate_mesh.py is stored]
Circle perimeter =  0.3078181289931019
Rectangle perimeter =  5.220000000000003

So far, so good. However, I would like to do the same thing that I did for the circle and the rectangle, but this time for a line. If I comment out # model.add_physical([channel_lines_1[3]], "Inflow") in generate_mesh.py in order to add the left edge of the rectangle, I get an error:

$ python3 generate_mesh.py 0.1
$ python3 mwe.py /home/fenics/shared/mesh/membrane_mesh
Traceback (most recent call last):
  File "mwe.py", line 25, in <module>
    infile.read(mvc, "name_to_read")
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:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------

May you please tell me what I am doing wrong ? Thank you :slight_smile:

2

The error you are getting indicates that gmsh has made some new lines in the process of making your plane surface. This means that the line you have tagged is not part of the mesh.

I would consider the mesh generation example at:
https://jsdokken.com/dolfinx-tutorial/chapter2/ns_code2.html#mesh-generation

3

Thank you, but I am using the legacy version and I don’t know how to read the example that you linked. Do you have an example that uses the legacy version?

I don’t see why the line that I have tagged is not part of the mesh, given that it connects points {0,0} and {0,h}, which are in the mesh. Do you have any ideas?

Thank you

4

If you just look at the Gmsh code, and ignore the dolfinx bit you will see how one should tag the mesh boundaries.

This is due to how Gmsh works. There is nothing I can do about how their backend adds and subtracts lines once you make surfaces.

You can merge duplicate entities with Gmsh functions as shown in how to obtain 1D entity information pre-meshing (post gmsh.model.occ.removeAllDuplicates()) (#2788) · Issues · gmsh / gmsh · GitLab

5

I still find this code unintellegible, but your suggestion

There is nothing I can do about how their backend adds and subtracts lines once you make surfaces

give me a useful (?) hint? I suspect that the model.add_physical(rectangle_Lh.curve_loop.curves, "Rectangle") in generate_mesh.py somewhat deranged the reading of the inflow line, because the rectangle was added as a physical element before the inflow line.

So I made the following change: I no longer include the boundary of the mesh with add_rectangle , but with the channel_lines itself, and then I include as a physical element the 2nd entry of channel_lines. Here are the codes:

#generate_mesh.py:

import numpy
import meshio
import gmsh
import pygmsh
import argparse

parser = argparse.ArgumentParser()
parser.add_argument("resolution")
args = parser.parse_args()

resolution = (float)(args.resolution)

# Channel parameters
L = 2.2
h = 0.41
r = 0.05
c_r = [0.2, 0.2, 0]

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

circle_r = model.add_circle(c_r, r, mesh_size=resolution/2)
rectangle_Lh = model.add_rectangle(0, L, 0, h, 0, mesh_size=resolution)

points = [model.add_point((0, 0, 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),
          model.add_point((0, h, 0), mesh_size=5*resolution)
          ]

channel_lines = [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_loop = model.add_curve_loop(channel_lines)

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

model.synchronize()
model.add_physical([plane_surface], "Volume")

#I will read this tagged element with `ds_circle = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=2)`
model.add_physical(circle_r.curve_loop.curves, "Circle")

#I will read this tagged element with `ds_rectangle = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=3)`
model.add_physical([channel_lines[1]], "Inflow")

geometry.generate_mesh(64)
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)

#mwe.py
import math
from fenics import *
from mshr import *
import numpy as np
import meshio
import ufl as ufl
import argparse

parser = argparse.ArgumentParser()
parser.add_argument("input_directory")
args = parser.parse_args()

L = 2.2
h = 0.41
r = 0.05
c_r = [0.2, 0.2]

#create mesh
mesh=Mesh()
with XDMFFile((args.input_directory) + "/triangle_mesh.xdmf") as infile:
    infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile((args.input_directory) + "/line_mesh.xdmf") as infile:
    infile.read(mvc, "name_to_read")
    
mf = dolfin.cpp.mesh.MeshFunctionSizet(mesh, mvc)

#here I read the tagged element declared as `model.add_physical(circle_r.curve_loop.curves, "Circle")` in generate_mesh.py
ds_circle = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=2)

#here I read the tagged element declared as model.add_physical([channel_lines_1[1]], "Inflow")
ds_inflow = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=3)

#here I integrate \int ds 1 over the circle and store the result of the integral as a double in inner_circumference
circle_length = assemble(1*ds_circle)
inflow_length = assemble(1*ds_inflow)
print("Circle length = ", circle_length)
print("Inflow length = ", inflow_length)

now if I run them, everything seems to work

$ python3 generate_mesh.py 0.1
$ python3 mwe.py /home/fenics/shared/mesh/membrane_mesh
$ python3 mwe.py /home/fenics/shared/mesh/membrane_mesh
Circle length =  0.3078181289931019
Inflow length =  0.41

May you please clarify why this is working now ? Also, now can I safely use ds_inflow in variational problems by coding a boundary part of the functional F = \langle \cdots \rangle_{\partial \Omega_{\rm inflow}} as
F = [...] * ds_inflow?

6

Could you please specify what about that mesh generation code you dont understand?
If we look at the relevant code:

gmsh.initialize()

L = 2.2
H = 0.41
c_x = c_y = 0.2
r = 0.05
gdim = 2
mesh_comm = MPI.COMM_WORLD
model_rank = 0
if mesh_comm.rank == model_rank:
    rectangle = gmsh.model.occ.addRectangle(0, 0, 0, L, H, tag=1)
    obstacle = gmsh.model.occ.addDisk(c_x, c_y, 0, r, r)
    fluid = gmsh.model.occ.cut([(gdim, rectangle)], [(gdim, obstacle)])
    gmsh.model.occ.synchronize()
    volumes = gmsh.model.getEntities(dim=gdim)
    assert (len(volumes) == 1)
    gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], fluid_marker)
    gmsh.model.setPhysicalName(volumes[0][0], fluid_marker, "Fluid")
inlet_marker, outlet_marker, wall_marker, obstacle_marker = 2, 3, 4, 5
inflow, outflow, walls, obstacle = [], [], [], []
if mesh_comm.rank == model_rank:
    boundaries = gmsh.model.getBoundary(volumes, oriented=False)
    for boundary in boundaries:
        center_of_mass = gmsh.model.occ.getCenterOfMass(boundary[0], boundary[1])
        if np.allclose(center_of_mass, [0, H / 2, 0]):
            inflow.append(boundary[1])
        elif np.allclose(center_of_mass, [L, H / 2, 0]):
            outflow.append(boundary[1])
        elif np.allclose(center_of_mass, [L / 2, H, 0]) or np.allclose(center_of_mass, [L / 2, 0, 0]):
            walls.append(boundary[1])
        else:
            obstacle.append(boundary[1])
    gmsh.model.addPhysicalGroup(1, walls, wall_marker)
    gmsh.model.setPhysicalName(1, wall_marker, "Walls")
    gmsh.model.addPhysicalGroup(1, inflow, inlet_marker)
    gmsh.model.setPhysicalName(1, inlet_marker, "Inlet")
    gmsh.model.addPhysicalGroup(1, outflow, outlet_marker)
    gmsh.model.setPhysicalName(1, outlet_marker, "Outlet")
    gmsh.model.addPhysicalGroup(1, obstacle, obstacle_marker)
    gmsh.model.setPhysicalName(1, obstacle_marker, "Obstacle")
    gmsh.model.mesh.generate(gdim)

This code creates a channel mesh with the relevant markers. This has nothing to do with dolfinx, and is pure Gmsh.

7

Sure, for example, here

gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], fluid_marker)
gmsh.model.setPhysicalName(volumes[0][0], fluid_marker, "Fluid")
  • I don’t know what the addPhysicalGroup does
  • I don’t know what volumes is
  • I don’t know what the [0][0] entry (?) of volumes is
  • … and plenty of other issues, both in these lines and in the other parts of the cited code.
8

There is only one volume in the mesh. You can verify this by printing volumes when running the code.
What I do is that I add a physical group for the volume/surface that we should mesh, to ensure that it is meshed properly by gmsh.
Each entry in volume is a tuple (dim_of_entity, tag), where one in 2D has (2, tag) for a cell.
Next we want to create a physical group for that tag, which we give the value fluid marker.

is optional to give it a string in the GMSH GUI.

Please note that these are pure gmsh functions, which can easily be looked up in the gmsh documentation.