Msh to xdmf (quadrilateral_9 not recognized

General
1

Hello,

I would like to convert a msh file generated with gmsh python API into two xdmf files (one for the full mesh and one for the facets). I use FEniCS (not FEniCSx). My script is organized as follow:

  1. I generate a msh file using gmsh python API. The mesh generation follows test problem 1. At the end I have a msh file called mesh.msh.

  2. Then I try to create a mesh & facet objects using tutorial 2. Inside the Create_mesh function, the mesh object I would like to return is called mesh and is a Mesh instance, the facet object is called facet and is a MeshFunctionSizet instance.

The whole script is written just below:

"""
This script takes only care about the gmsh mesh creation
"""

############################################
#------------ Mesh Creation ---------------#
############################################
import gmsh
import time
import sys
import numpy as np
import math
from dolfin import *
import meshio

# cylinder radius
r = 0.05

# length of the fluid domain
L = 2.2

# hight
h = 0.41

# coordinates of the center of the cylinder
c_x, c_y = [0.2], [0.2]

# dimension of the mesh
gdim = 2

# Calculating the time that mesh generation will take
tic = time.perf_counter()

gmsh.initialize()

gmsh.model.add("Schaffer_cylinder")

# target mesh size for any point
lc = 0.2
# lc = 0.75

# extra fine 
# factors = {"cylinders": 0.025, "cylinders_far": 0.6, "wake_start": 0.05, "wake_end": 0.5}

# very_fine (with l=1 --> 8791)
# factors = {"cylinders": 0.045, "cylinders_far": 1, "wake_start": 0.05, "wake_end": 0.25}

# fine (with lc = 0.35 = 52246)
# factors = {"cylinders": 0.045, "cylinders_far": 0.75, "wake_start": 0.055, "wake_end": 0.15}

# # mirror comparison
# factors = {"cylinders": 0.045, "cylinders_far": 1, "wake_start": 0.055, "wake_end": 0.35}

# test settings
factors = {"cylinders": 0.02, "cylinders_far": 0.2, "wake_start": 0.02, "wake_end": 0.2}

# standard
# factors = {"cylinders": 0.07, "cylinders_far": 1, "wake_start": 0.08, "wake_end": 0.4}

# coarse
# factors = {"cylinders": 0.1, "cylinders_far": 1, "wake_start": 0.12, "wake_end": 0.4}

# very coarse
# factors = {"cylinders": 0.14, "cylinders_far": 1, "wake_start": 0.19, "wake_end": 0.4}

############################################
#---------- Geometry creation -------------#
############################################

# Rectangle
rectangle = gmsh.model.occ.addRectangle(0, 0, 0, L, h, tag=1)

# Cylinder
obstacle = gmsh.model.occ.addDisk(c_x[0], c_y[0], 0, r, r)

# Fluid domain
fluid = gmsh.model.occ.cut([(gdim, rectangle)], [(gdim, obstacle)]) 
gmsh.model.occ.synchronize()




############################################
#--------- Tags & Physical Groups ---------#
############################################

fluid_marker = 1
volumes = gmsh.model.getEntities(dim=gdim)
assert (len(volumes) == 1)
gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], tag=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 = [], [], [], []
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, tag=wall_marker)
gmsh.model.setPhysicalName(1,wall_marker, "walls")
gmsh.model.addPhysicalGroup(1, inflow, tag=inlet_marker)
gmsh.model.setPhysicalName(1, inlet_marker, "inlet")
gmsh.model.addPhysicalGroup(1, outflow, tag=outlet_marker)
gmsh.model.setPhysicalName(1, outlet_marker, "outlet")
gmsh.model.addPhysicalGroup(1, obstacle, tag=obstacle_marker)
gmsh.model.setPhysicalName(1, obstacle_marker, "cylinder")



############################################
#---------------- Fields ------------------#
############################################

res_min = r / 5
distance_field = gmsh.model.mesh.field.add("Distance")
gmsh.model.mesh.field.setNumbers(distance_field, "EdgesList", obstacle)
threshold_field = gmsh.model.mesh.field.add("Threshold")
gmsh.model.mesh.field.setNumber(threshold_field, "IField", distance_field)
gmsh.model.mesh.field.setNumber(threshold_field, "LcMin", res_min)
gmsh.model.mesh.field.setNumber(threshold_field, "LcMax", 0.15 * h)
gmsh.model.mesh.field.setNumber(threshold_field, "DistMin", r)
gmsh.model.mesh.field.setNumber(threshold_field, "DistMax", 2 * h)
min_field = gmsh.model.mesh.field.add("Min")
gmsh.model.mesh.field.setNumbers(min_field, "FieldsList", [threshold_field])
gmsh.model.mesh.field.setAsBackgroundMesh(min_field)


############################################
#------------ Mesh Generation -------------#
############################################

gmsh.option.setNumber("Mesh.Algorithm", 8)
gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2)
gmsh.option.setNumber("Mesh.RecombineAll", 1)
gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.setOrder(2)
gmsh.model.mesh.optimize("Netgen")

gmsh.write("mesh.msh")


# gmsh.model.mesh.generate(gdim)
if '-nopopup' not in sys.argv:
    gmsh.fltk.run()

############################################
#---------- Creating xdmf Files -----------#
# ------- and mesh & facet objects --------#
############################################

def create_Mesh():
	"""
	Method that creates the Mesh object from a msh file. It also creates facet_tags.

	Returns
	-------
	the mesh object
	"""
	mesh = Mesh()

	# The mesh generated with gmsh's python API
	mesh_file = "mesh.msh"

	# xdmf file that will contain the mesh
	# xdmf_file = os.path.splitext(mesh_file)[0] + ".xdmf"
	xdmf_file = "full_mesh.xdmf"

	# xdmf file that will contain the 1D elements mesh (facets)
	# facet_file = os.path.splitext(mesh_file)[0] + "_facet.xdmf"
	facet_file = "facet_mesh.xdmf"

	# read mesh from msh file
	msh = meshio.read(mesh_file)

	# create the 1D element mesh
	# line_mesh = create_mesh(msh, "line", prune_z=True)
	line_mesh = create_mesh(msh, "line3", prune_z=True)
	meshio.write(facet_file, line_mesh)


	# TODO set 
	# create the xdmf mesh
	# triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
	# meshio.write(xdmf_file, triangle_mesh)
	quad_mesh = create_mesh(msh, "quad9", prune_z=True)
	meshio.write(xdmf_file, quad_mesh)

	# writing xdmf file into Mesh object
	with XDMFFile(xdmf_file) as infile:
		infile.read(mesh)

	# retrieving the facets
	mvc1 = MeshValueCollection("size_t", mesh, 1)
	with XDMFFile(facet_file) as infile:
		infile.read(mvc1, "name_to_read")
	facet = cpp.mesh.MeshFunctionSizet(mesh, mvc1)

	# reading mesh from file
	# self.parameters["mesh"] = mesh
	# self.parameters["facet_tags"] = facet

	return mesh, facet



def create_mesh(mesh, cell_type, data_name="name_to_read", prune_z=False):
	"""
	Function that create a meshio.Mesh object from msh mesh using meshio.

	Parameters
	----------
	mesh meshio.Mesh
		the mesh object on which the mesh will be written
	cell_type : str
		see meshio doc
	data_name : str
		see meshio doc
	prune_z: bool
		see meshio doc

	Returns
	-------
	meshio.Mesh
		the mesh created from the msh mesh file.
	"""
	# retrieving the cell type (line or triangles or quads)
	cells = mesh.get_cells_type(cell_type)

	# getting the cell data knowing physical groups of the msh file and the cell type
	cell_data = mesh.get_cell_data("gmsh:physical", cell_type)

	# retrieving the points corresponding to the cells (the z coordinate is pruned in the 2D case!)
	points = mesh.points[:, :2] if prune_z else mesh.points

	# creating the meshio.Mesh object from the previously defined parameters
	out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={
						   data_name: [cell_data]})
	return out_mesh


full_mesh, facet_mesh = create_Mesh()

The mesh (following test problem 1) is made of second order quadrilateral elements. In meshio, I checked, and the cells type are "line3" or "quad9". So I modified the tutorial 2, and changed "line" to "line3" and "quad" to "quad9".

The problem occurs when the class XDMFFile read the xdmf file created with meshio. The error specifies:

*** Error:   Unable to recognise cell type.
*** Reason:  Unknown value "quadrilateral_9".

Also note that if I comment the line gmsh.model.mesh.setOrder(2) and change back "line3" to "line" and "quad9" to "quad", then I get this error:

*** Error:   Unable to order quadrilateral cell.
*** Reason:  Cell is not orderable.
*** Where:   This error was encountered inside QuadrilateralCell.cpp.
*** Process: 0

The tutorial 2 is useful to me, but I must admit that I do not really know how it works. As a consequence, my modification do not work. Could someone help me adapting the tutorial 2 so that I get the mesh and facets (with line3 & quad9)? Because with a standard 2D triangle mesh, everything works fine.
Thanks and regards.

2

Legacy FEniCS does not support higher order quadrilateral cells from file.

3

I see. But now, commenting gmsh.model.mesh.setOrder(2), I get a mesh with elements of order 1 right ? Then how to avoid the Unable to order quadrilateral cell error ? Is tutorial 2 only fit for triangle elements ?

4

I think you will have a general issue reading in unordered quadrilateral meshes in legacy Dolfin, due to the issues described in 1.1 of: https://arxiv.org/pdf/2102.11901.pdf
as reading from XDMF doesn’t expose the turn of reordering parameter: Bitbucket
as shown in:
Cell node labeling - #2 by dokken

5

Ok so, as I understand, I should abandon the quadrilateral type mesh, and opt for a regular triangle mesh ?

6

If you want to use legacy fenics, yes.