Tranfinite gmsh mesh and read_meshtags lead to unexpected crash

flo1 · February 22, 2024, 3:25pm

Good afternoon everyone,

I decided to generate a transfinite mesh with gmsh and use the standard routine that you can find in this tutorial: Defining subdomains for different materials — FEniCSx tutorial to get tagged physical group. The weirdest thing is if you put a break point inside the last python with statement you will get the tags. But then, I’ve got either:

PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range, …
Python malloc error with missing objects,
…
And once it worked, but I got leaked semaphore at the end…

And here is a long MWE, sorry I’m not using open Cascade so the gmsh part is quite long…

import pathlib

import gmsh
from dolfinx import io
from dolfinx.mesh import MeshTagsMetaClass
from meshio import Mesh as MeshIO
from meshio import read, write
from mpi4py.MPI import COMM_SELF, COMM_WORLD


def create_mesh(meshio_mesh: MeshIO, cell_type: str, prune_z=False) -> MeshIO:
    cells = meshio_mesh.get_cells_type(cell_type)
    cell_data = meshio_mesh.get_cell_data("gmsh:physical", cell_type)
    points = meshio_mesh.points[:, : 2] if prune_z else meshio_mesh.points
    out_mesh = MeshIO(points=points, cells={cell_type: cells}, cell_data={"name_to_read": [cell_data]})
    return out_mesh


L = 160.0
H = 40.0
mesh_coef = 1
vertical_cells = int(H * mesh_coef) if int(H * mesh_coef) % 2 == 0 else int(H * mesh_coef) + 1
horizontal_cells = int(round(vertical_cells * L / H, 0))

gmsh_file = pathlib.Path.cwd().joinpath('mesh').with_suffix('.msh')
gmsh.initialize()

gmsh.model.add('poutre_1_2D')
# Add points for the rectangular
p_sw = gmsh.model.geo.addPoint(0, 0, 0, 1.0)
p_nw = gmsh.model.geo.addPoint(0, H, 0, 1.0)
p_ne = gmsh.model.geo.addPoint(L, H, 0, 1.0)
p_se = gmsh.model.geo.addPoint(L, 0, 0, 1.0)
# Add a point in the middle of the right face
p_middle = gmsh.model.geo.addPoint(L, H / 2, 0, 1.0)
# Add lines
l_left = gmsh.model.geo.addLine(p_sw, p_nw)
l_upper = gmsh.model.geo.addLine(p_nw, p_ne)
l_right = gmsh.model.geo.addLine(p_se, p_ne)
l_bottom = gmsh.model.geo.addLine(p_se, p_sw)
# Add curve loop
cl = gmsh.model.geo.addCurveLoop([-l_left, -l_upper, l_right, -l_bottom])
# Add plane surface
ps = gmsh.model.geo.addPlaneSurface([cl])
# Synchronization
gmsh.model.geo.synchronize()
# Adding physical groups
# Tag the middle point
middle_point = gmsh.model.geo.addPhysicalGroup(0, [p_middle], name="point_force")
# Tag the left face
left_face = gmsh.model.geo.addPhysicalGroup(1, [l_left], name="clamped")
# Tag all the surface ps
all_surfaces = gmsh.model.geo.addPhysicalGroup(2, [ps], name="design_domain")
# Synchronization
gmsh.model.geo.synchronize()
# Transfinite mesh
# Constraints explicitly the location of the nodes on the curve.
gmsh.model.geo.mesh.setTransfiniteCurve(l_left, vertical_cells + 1)
gmsh.model.geo.mesh.setTransfiniteCurve(l_right, vertical_cells + 1)
gmsh.model.geo.mesh.setTransfiniteCurve(l_upper, horizontal_cells + 1)
gmsh.model.geo.mesh.setTransfiniteCurve(l_bottom, horizontal_cells + 1)
# Uses a transfinite interpolation algorithm in the parametric plane of the surface to connect the nodes on the
# boundary using a structured grid.
gmsh.model.geo.mesh.setTransfiniteSurface(ps)
# To create quadrangles instead of triangles, one can use the 'setRecombine' constraint:
gmsh.model.geo.mesh.setRecombine(2, ps)
# Before meshing the model, we need to use the synchronize command
gmsh.model.geo.synchronize()
gmsh.model.mesh.generate(2)
# Save the mesh
gmsh.write(str(gmsh_file))
# Generate the mesh
gmsh.finalize()

meshio_mesh: MeshIO = read(str(gmsh_file))
# First line
meshio_to_save = create_mesh(meshio_mesh, 'line', True)
file_dir_line = pathlib.Path.cwd().joinpath('line_mesh')
write(str(file_dir_line.with_suffix('.xdmf')), meshio_to_save)
# Then quad
meshio_to_save = create_mesh(meshio_mesh, 'quad', True)
file_dir_quad = pathlib.Path.cwd().joinpath('quad_mesh')
write(str(file_dir_quad.with_suffix('.xdmf')), meshio_to_save)
# Load mesh: ok
with io.XDMFFile(COMM_SELF, file_dir_quad.with_suffix('.xdmf'), 'r', io.XDMFFile.Encoding.HDF5) as xdmf:
    # name="Grid" for gmsh and name="mesh" for fenicsx
    domain = xdmf.read_mesh(name='Grid')
print('Domain was load')
# Get line tags: bad
with io.XDMFFile(COMM_SELF, str(file_dir_line.with_suffix('.xdmf')), "r", io.XDMFFile.Encoding.HDF5) as xdmf:
    tags: MeshTagsMetaClass = xdmf.read_meshtags(domain, name="Grid")
print('Tags were read')

PSs:

I’m running (still) on fenicsx 0.6. But on 0.7 it gives me the same error. I’m working with a Mac, but I got the same error on a Linux platform (with fenicsx 0.6).
And yes the code works without this transfinite stuff, and it gives me the same errors with 3D transfinite mesh.

dokken · February 22, 2024, 4:46pm

Remove the physical marker on the center point of the facet. It creates an extra point in your mesh that shouldn’t be there (im not sure why gmsh creates an extra point though).
Also added create connectivity:

import pathlib

import gmsh
from dolfinx import io
from meshio import Mesh as MeshIO
from meshio import read, write
from mpi4py.MPI import COMM_SELF, COMM_WORLD


def create_mesh(meshio_mesh: MeshIO, cell_type: str, prune_z=False) -> MeshIO:
    cells = meshio_mesh.get_cells_type(cell_type)
    cell_data = meshio_mesh.get_cell_data("gmsh:physical", cell_type)
    points = meshio_mesh.points[:, : 2] if prune_z else meshio_mesh.points
    out_mesh = MeshIO(points=points, cells={cell_type: cells}, cell_data={"name_to_read": [cell_data]})
    return out_mesh


L = 160.0
H = 40.0
mesh_coef = 1
vertical_cells = int(H * mesh_coef) if int(H * mesh_coef) % 2 == 0 else int(H * mesh_coef) + 1
horizontal_cells = int(round(vertical_cells * L / H, 0))

gmsh_file = pathlib.Path.cwd().joinpath('mesh').with_suffix('.msh')
gmsh.initialize()

gmsh.model.add('poutre_1_2D')
# Add points for the rectangular
p_sw = gmsh.model.geo.addPoint(0, 0, 0, 1.0)
p_nw = gmsh.model.geo.addPoint(0, H, 0, 1.0)
p_ne = gmsh.model.geo.addPoint(L, H, 0, 1.0)
p_se = gmsh.model.geo.addPoint(L, 0, 0, 1.0)
# Add a point in the middle of the right face
p_middle = gmsh.model.geo.addPoint(L, H / 2, 0, 1.0)
# Add lines
l_left = gmsh.model.geo.addLine(p_sw, p_nw)
l_upper = gmsh.model.geo.addLine(p_nw, p_ne)
l_right = gmsh.model.geo.addLine(p_se, p_ne)
l_bottom = gmsh.model.geo.addLine(p_se, p_sw)
# Add curve loop
cl = gmsh.model.geo.addCurveLoop([-l_left, -l_upper, l_right, -l_bottom])
# Add plane surface
ps = gmsh.model.geo.addPlaneSurface([cl])
# Synchronization
gmsh.model.geo.synchronize()
# Adding physical groups
# Tag the middle point
#middle_point = gmsh.model.geo.addPhysicalGroup(0, [p_middle], name="point_force")
# Tag the left face
left_face = gmsh.model.geo.addPhysicalGroup(1, [l_left], name="clamped")
# Tag all the surface ps
all_surfaces = gmsh.model.geo.addPhysicalGroup(2, [ps], name="design_domain")
# Synchronization
gmsh.model.geo.synchronize()
# Transfinite mesh
# Constraints explicitly the location of the nodes on the curve.
gmsh.model.geo.mesh.setTransfiniteCurve(l_left, vertical_cells + 1)
gmsh.model.geo.mesh.setTransfiniteCurve(l_right, vertical_cells + 1)
gmsh.model.geo.mesh.setTransfiniteCurve(l_upper, horizontal_cells + 1)
gmsh.model.geo.mesh.setTransfiniteCurve(l_bottom, horizontal_cells + 1)
# Uses a transfinite interpolation algorithm in the parametric plane of the surface to connect the nodes on the
# boundary using a structured grid.
gmsh.model.geo.mesh.setTransfiniteSurface(ps)
# To create quadrangles instead of triangles, one can use the 'setRecombine' constraint:
gmsh.model.geo.mesh.setRecombine(2, ps)
# Before meshing the model, we need to use the synchronize command
gmsh.model.geo.synchronize()
gmsh.model.mesh.generate(2)
# Save the mesh
gmsh.write(str(gmsh_file))
# Generate the mesh
gmsh.finalize()

meshio_mesh: MeshIO = read(str(gmsh_file))
# First line
meshio_to_save = create_mesh(meshio_mesh, 'line', True)
file_dir_line = pathlib.Path.cwd().joinpath('line_mesh')
write(str(file_dir_line.with_suffix('.xdmf')), meshio_to_save)
# Then quad
meshio_to_save = create_mesh(meshio_mesh, 'quad', True)
file_dir_quad = pathlib.Path.cwd().joinpath('quad_mesh')
write(str(file_dir_quad.with_suffix('.xdmf')), meshio_to_save)
# Load mesh: ok
with io.XDMFFile(COMM_SELF, file_dir_quad.with_suffix('.xdmf'), 'r', io.XDMFFile.Encoding.HDF5) as xdmf:
    # name="Grid" for gmsh and name="mesh" for fenicsx
    domain = xdmf.read_mesh(name='Grid')
print('Domain was load')
# Get line tags: bad
domain.topology.create_connectivity(domain.topology.dim - 1, domain.topology.dim)
with io.XDMFFile(COMM_SELF, str(file_dir_line.with_suffix('.xdmf')), "r", io.XDMFFile.Encoding.HDF5) as xdmf:
    tags = xdmf.read_meshtags(domain, name="Grid")

flo1 · February 26, 2024, 8:45pm

Many thanks @dokken for your answer! Everything works fine now. I will try to use Open Cascade and divide my geometry in multiple part, so this point must be here !

Have a good evening!

dokken · February 26, 2024, 8:54pm

You can embed the point: How to set 'boundary condition' inside the boundary? - #11 by dokken

LGerard · February 26, 2024, 10:30pm

Hello

I exported a mesh from gmsh, but I have problems to imposse a boundary condicition

Note, the file .msh it was created from gmsh software

import gmsh
import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.fem import (Constant, Function, FunctionSpace,
                         assemble_scalar, dirichletbc, form, locate_dofs_topological, set_bc, locate_dofs_geometrical)
from dolfinx.io import (VTXWriter, distribute_entity_data, gmshio)
from ufl import (FacetNormal, FiniteElement, Identity, Measure, TestFunction, TrialFunction, VectorElement,
                 as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym)

mesh, cell_tags, facet_tags = gmshio.read_from_msh(
    "mesh_v1.msh", MPI.COMM_WORLD, 0, gdim=2)

gmsh.initialize()
if MPI.COMM_WORLD.rank == 0:
    gmsh.model.add("Mesh from file")
    gmsh.merge("mesh_v1.msh")
output = gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=2)
gmsh.finalize()

v_cg2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
s_cg1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, v_cg2)
Q = FunctionSpace(mesh, s_cg1)

u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)

u_ref = 1.2
z_ref = 2
class InletVelocity():
    def __call__(self, x):
        values = np.zeros((2, x.shape[1]), dtype=PETSc.ScalarType)
        values[0] = u_ref * pow(x[1] / z_ref, 0.16)
        return values

# Inlet
print(facet_tags)
u_inlet = Function(V)
inlet_velocity = InletVelocity()
u_inlet.interpolate(inlet_velocity)
bcu_inflow = dirichletbc(u_inlet, locate_dofs_topological(V, 1, facet_tags))

# Walls
u_nonslip = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType)

dokken · February 27, 2024, 3:50am

Without having the mesh or error message with stack trace, I see no way of giving further guidance.

LGerard · February 27, 2024, 12:59pm

Traceback (most recent call last):
File "/home/gerard/Dropbox/FEniCSx/turbulent_model.py", line 44, in <module>
bcu_inflow = dirichletbc(u_inlet, locate_dofs_topological(V, 1, facet_tags))
File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/bcs.py", line 81, in locate_dofs_topological
_entities = np.asarray(entities, dtype=np.int32)
TypeError: int() argument must be a string, a bytes-like object or a real number, not 'MeshTags'

Two question, 1. In the line, is the “cell_tags” and “facet_tags” similar to “bundaries” and “markers” in FEniCS?

mesh, cell_tags, facet_tags = gmshio.read_from_msh(
    "mesh_v1.msh", MPI.COMM_WORLD, 0, gdim=2)

DOLFINx is designed to run in parallel. My cuestion is, there are any difference when I provide the mesh using gmsh and jus read it, compared to creating the mesh with the API and meshio? Note that I have a complex mesh.

flo1 · February 27, 2024, 1:36pm

Thanks ! Yes, I’ve already used the gmsh.model.mesh.embed (and tag vertices) command for the non-transfinite mesh and it works fine. But with transfinite mesh I got the same error. Normally, a node exists as I set the good parity for the number of vertical cells, but the position is slightly different. But, it’s OK, I’ve created a json file to get the wanted data and I might divide my geometry in smaller bits, to have this node that equates the mesh vertices.

flo1 · February 27, 2024, 2:27pm

If the type of facet_tags is mesh.MeshTagsMetaClass. You need to add this line before bcu_inflow

features = facet_tags.find(int(your tag number for this DC condition))
bcu_inflow = dirichletbc(u_inlet, locate_dofs_topological(V, 1, features))

Have a look to the tutorial I’ve posted in the first comment.

flo1 · February 27, 2024, 2:29pm

I’ve an other question. Is the order in topology.create_connectivity is important ? It’s to create the Agency List that we can find in topology.connectivity ?

dokken · February 27, 2024, 6:36pm

LGerard:

Traceback (most recent call last):
File "/home/gerard/Dropbox/FEniCSx/turbulent_model.py", line 44, in <module>
bcu_inflow = dirichletbc(u_inlet, locate_dofs_topological(V, 1, facet_tags))
File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/bcs.py", line 81, in locate_dofs_topological
_entities = np.asarray(entities, dtype=np.int32)
TypeError: int() argument must be a string, a bytes-like object or a real number, not 'MeshTags'

You should send in facets tags directly as it is a meshtags object. You should send in the set of indices in the meshtags object (facet_tags.indices).

Again, you have not provided a reproducible example, so this requires extra effort from my end to decipher your code. I will not provide further assistant if you do not do your best to make the problem you are experiencing reproducible.

These are not classes in dolfin, so i can only assume from context that you mean MeshFunction objects for cells and facets. And yes, to some extent they are equivalent. See: Transition from MeshFunction in Fenicsx - #2 by dokken for more context.

If you write the mesh to xdmf with meshio or similar, the reading in of the mesh will be more efficient, as one can use the fact that xdmf/h5 is a binary format, and read in the data in a partitioned fashion. It would also avoid recreating the mesh every time.

dokken · February 27, 2024, 6:37pm

You need to call topology.create_connectivity before calling certain functions in DOLFINx. Most functions will throw errors if these aren’t created when they are needed.

LGerard · February 28, 2024, 2:44pm

My apologies, Dokken. I didn’t know how to share the mesh I generated with Gmsh. I created a GitHub repository to share it.

https://github.com/LGerard-gg/FEniCSx

I’m using the Gmsh software because the study region will be complex. However, at the moment, it’s a rectangle.

https://jsdokken.com/src/tutorial_gmsh.html

I relied on this tutorial to extract mesh data, although I’m not quite sure if cell_tags and facet_tags are what I expect. With the ‘gmsh.read_from_msh’ in the line

mesh, cell_tags, facet_tags = gmshio.read_from_msh(
    "mesh_v1.msh", MPI.COMM_WORLD, 0, gdim=2)

I marked the boundaries with Gmsh, so I wanted to see that labeling with an xdmf file, but I encountered an error related to the geometry. Now, however, it’s a different error.

LGerard · February 28, 2024, 2:44pm

This is my code, with the mesh shared in the previous message.

import gmsh
import numpy as np
import dolfinx
from dolfinx.cpp.mesh import Geometry_float32, Geometry_float64
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.fem import (Constant, Function, FunctionSpace,
                         assemble_scalar, dirichletbc, form, locate_dofs_topological, set_bc, locate_dofs_geometrical)
from dolfinx.io import (VTXWriter, distribute_entity_data, gmshio, XDMFFile)
from ufl import (FacetNormal, FiniteElement, Identity, Measure, TestFunction, TrialFunction, VectorElement,
                 as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym)

mesh, cell_tags, facet_tags = gmshio.read_from_msh(
        "mesh_v1.msh", MPI.COMM_WORLD, 0, gdim=2)

gmsh.initialize()
if MPI.COMM_WORLD.rank == 0:
    gmsh.model.add("Mesh from file")
    gmsh.merge("mesh_v1.msh")
output = gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=2)
gmsh.finalize()

v_cg2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
s_cg1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, v_cg2)
Q = FunctionSpace(mesh, s_cg1)

u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)

u_ref = 1.2
z_ref = 2
class InletVelocity():
    def __call__(self, x):
        values = np.zeros((2, x.shape[1]), dtype=PETSc.ScalarType)
        values[0] = u_ref * pow(x[1] / z_ref, 0.16)
        return values

# Inlet
print(facet_tags)
u_inlet = Function(V)
inlet_velocity = InletVelocity()
u_inlet.interpolate(inlet_velocity)
features = facet_tags.find(1)
bcu_inflow = dirichletbc(u_inlet, locate_dofs_topological(V, 1, features))
print(features)
# Walls
u_nonslip = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType)

with dolfinx.io.XDMFFile(mesh.comm, "ft.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(facet_tags, dolfinx.cpp.mesh.Geometry_float32, str, str = '/Xdmf/Domain')

dokken · March 20, 2024, 7:32pm

LGerard:

with dolfinx.io.XDMFFile(mesh.comm, "ft.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(facet_tags, dolfinx.cpp.mesh.Geometry_float32, str, str = '/Xdmf/Domain')

Simply change this to

with dolfinx.io.XDMFFile(mesh.comm, "ft.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(
        facet_tags, mesh.geometry)

Topic		Replies	Views
Issue reading the mesh tags imported from gmsh General mesh , boundary-conditions , dolfinx , mpi , gmsh	7	162	July 9, 2024
Cell Tags in the Wrong Order mesh	3	720	July 18, 2022
GMSH Transfinite Volume & Fenics I/0 I/O	6	2609	March 10, 2022
Seg fault when using mesh from gmshio.read_from_msh General	6	22	October 8, 2024
Read meshtag from xdmf created by <mesh generation with gmsh> I/O mesh , dolfinx , gmsh	2	62	July 20, 2024

Tranfinite gmsh mesh and read_meshtags lead to unexpected crash

Related topics