Issue with quadrilateral mesh orientation when importing SALOME mesh into FEniCS 2019.1

mesh
1

Hi all,

I created my quadrilateral mesh using SALOME and then imported it into my code. After generating the .xdmf and .h5 files, I encountered an issue when trying to use them: there seems to be a mismatch between the orientation in the .xdmf file and what FEniCS 2019 expects.

Here is the error message I get:

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 order quadrilateral cell.
 Reason:  Cell is not orderable.
 Where:   This error was encountered inside QuadrilateralCell.cpp.
 Process: 0

DOLFIN version: 2019.1.0
Git changeset:  ec7503cef5f328d4927ba1050abe55c3519ac4e9

Could you please help me solve this problem?
I’m also attaching my conversion script used to transform the .med file into .xdmf/.h5 format:

 msh = meshio.read(".../Mesh_Tun_EURO_FLiBe_NEW_quadr.med")
 print("CELLS:", msh.cells)
 print("CELL DATA:", msh.cell_data_dict.keys())
 print("CELL SETS:", msh.cell_sets)
 print("FIELD DATA:", msh.field_data)
 
 def check_unorderable_quads(points, quads):
     bad_cells = [  ]
     for i, q in enumerate(quads):
         pts = points[q, :2]
         area = 0.5 * (
             (pts[0,0]*pts[1,1] + pts[1,0]*pts[2,1] +
              pts[2,0]*pts[3,1] + pts[3,0]*pts[0,1]) -
             (pts[1,0]*pts[0,1] + pts[2,0]*pts[1,1] +
              pts[3,0]*pts[2,1] + pts[0,0]*pts[3,1])
         )
         if abs(area) < 1e-12 or np.isnan(area):
             bad_cells.append(i)
     return bad_cells
 
 quads = msh.cells_dict["quad"]
 bad = check_unorderable_quads(msh.points, quads)
 
 def fix_quad_order(points, quads):
     
     fixed_quads = []
     inverted = 0
     for q in quads:
         pts = points[q, :2]
         # Calcola il baricentro
         cx, cy = pts.mean(axis=0)
         angles = np.arctan2(pts[:, 1] - cy, pts[:, 0] - cx)
         order = np.argsort(angles) 
         q_new = q[order]
         area = (
             (pts[1,0] - pts[0,0]) * (pts[1,1] + pts[0,1]) +
             (pts[2,0] - pts[1,0]) * (pts[2,1] + pts[1,1]) +
             (pts[3,0] - pts[2,0]) * (pts[3,1] + pts[2,1]) +
             (pts[0,0] - pts[3,0]) * (pts[0,1] + pts[3,1])
         )
         if area > 0:  
             q_new = q_new[[0, 3, 2, 1]]
             inverted += 1
         q_new = q_new[[0, 1, 2, 3]]
         fixed_quads.append(q_new)
     return np.array(fixed_quads, dtype=int)
 
 
 def convert_med_to_xdmf(
     med_file,
     cell_file="mesh_domains.xdmf",
     facet_file="mesh_boundaries.xdmf",
     cell_type="quad",
     facet_type="line",
 ):
     """
     Converte un file MED in XDMF correggendo l'ordine dei vertici
     per compatibilità con FEniCS 2019.1 (quadrilateri antiorari)."""
     
     msh = meshio.read(med_file)
 
     correspondance_dict = msh.cell_tags
     cell_data_types = msh.cell_data_dict["cell_tags"].keys()
 
     for mesh_block in msh.cells:
         if mesh_block.type == cell_type:
             quads = mesh_block.data
             quads_fixed = fix_quad_order(msh.points, quads)
 
             meshio.write_points_cells(
                 cell_file,
                 msh.points,
                 [meshio.CellBlock(cell_type, quads_fixed)],
                 cell_data={"f": [-1 * msh.cell_data_dict["cell_tags"][cell_type]]},
             )
 
         elif mesh_block.type == facet_type:
             meshio.write_points_cells(
                 facet_file,
                 msh.points,
                 [mesh_block],
                 cell_data={"f": [-1 * msh.cell_data_dict["cell_tags"][facet_type]]},  )
     return correspondance_dict, cell_data_types
 correspondance_dict, cell_data_types = convert_med_to_xdmf(
     "...Mesh_Tun_EURO_FLiBe_NEW_quadr.med",
     cell_file=".../mesh_domains.xdmf",
     facet_file=".../mesh_boundaries.xdmf"
 )
 mesh = meshio.read(".../mesh_domains.xdmf")
 print(mesh)
 
mesh_bnd = meshio.read(".../mesh_boundaries.xdmf")
 print(mesh_bnd)

Thank you very much for your help!

Best regards,
Andrea.

2

The support for quadrilateral meshes are very limited in FEniCS 2019.1.0.

I would strongly advice you to move to DOLFINx, as the support for arbitrary order triangle, tetrahedral, quadrilateral and hexahedral elements have been greatly enhanced.

1 Like