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.