Msh to xml for quad elements

bagla0 · March 16, 2023, 9:10pm

Hi FEniCS family,

I am trying to convert .msh file to xml with physical region and facets region using command
dolfin-convert test.msh mesh.xml

Converting from Gmsh format (.msh, .gmsh) to DOLFIN XML format
*** Unable to find cells of supported type.

The mesh contains quad elements with order 2. How can I get the command work to find quad cells ? I tried using meshio, but, that not able to find right recourse of code.

Any help is greatly appreciated.

Kindly find the attached mesh.

dokken · March 16, 2023, 9:58pm

There is no proper support for quad elements in legacy Dolfin. See dolfinx

For such support.

This is for instance used in
https://jsdokken.com/dolfinx-tutorial/chapter2/ns_code2.html#loading-mesh-and-boundary-markers

dokken · March 16, 2023, 10:14pm

As you have only supplied a picture of the mesh, not the actual geo file generating the msh-file, I cannot help you further.

In general, dolfin-convert has been deprecated for years, and I’ve written many posts on this forum on how to use meshio to convert meshes for legacy dolfin.
You could also find this at: Mesh generation and conversion with GMSH and PYGMSH | Jørgen S. Dokken

bagla0 · March 16, 2023, 10:33pm

Please find the geo file used to form .msh file. I am getting nan values for all analysis. The major problem, I am facing is due to physical facets. The present case shows physical facets at the square boundary edges only.

What line region should be included in Physical Curve ?

Kindly suggest me ways to resolve the physical facets. i tried a lot but, not getting the issue.
r = 0.282095;

Point(1) = {-0.5, -0.5, 0, 0.02};
Point(2) = {0.5, -0.5, 0, 0.02};
Point(3) = {0.5, 0.5, 0, 0.02};
Point(4) = {-0.5, 0.5, 0, 0.02};
Point(5) = {-0.5, 0, 0, 0.02};
Point(6) = {0, -0.5, 0, 0.02};
Point(7) = {0.5, 0, 0, 0.02};
Point(8) = {0, 0.5, 0, 0.02};
Point(9) = {0, 0, 0, 0.02};
Point(10) = {0, 0.282095, 0, 0.02};
Point(11) = {0, -0.282095, 0, 0.02};
Point(12) = {0.282095, 0, 0, 0.02};
Point(13) = {-0.282095, 0, 0, 0.02};
Point(14) = {-r/2Sqrt(2),r/2Sqrt(2), 0, 0.02};
Point(15) = {-r/2Sqrt(2), -r/2Sqrt(2), 0, 0.02};
Point(16) = {r/2Sqrt(2), -r/2Sqrt(2), 0, 0.02};
Point(17) = {r/2Sqrt(2), r/2Sqrt(2), 0, 0.02};
Point(18) = {-r/2Sqrt(2)/2, r/2Sqrt(2)/2, 0, 0.02};
Point(19) = {-r/2Sqrt(2)/2, -r/2Sqrt(2)/2, 0, 0.02};
Point(20) = {r/2Sqrt(2)/2, -r/2Sqrt(2)/2, 0, 0.02};
Point(21) = {r/2Sqrt(2)/2, r/2Sqrt(2)/2, 0, 0.02};
Point(22) = {-r/2Sqrt(2)/2, 0, 0, 0.02};
Point(23) = {0, -r/2Sqrt(2)/2, 0, 0.02};
Point(24) = {r/2Sqrt(2)/2, 0, 0, 0.02};
Point(25) = {0, r/2Sqrt(2)/2, 0, 0.02};
Line(3) = {1, 6};
Line(4) = {6, 2};
Line(5) = {2, 7};
Line(6) = {7, 3};
Line(7) = {3, 8};
Line(8) = {8, 4};
Line(9) = {4, 5};
Line(10) = {5, 1};
Line(11) = {18, 22};
Line(12) = {22, 19};
Line(13) = {19, 23};
Line(14) = {23, 20};
Line(15) = {20, 24};
Line(16) = {24, 21};
Line(17) = {21, 25};
Line(18) = {25, 18};
Line(19) = {18, 14};
Line(20) = {14, 4};
Line(21) = {22, 13};
Line(22) = {13, 5};
Line(23) = {19, 15};
Line(24) = {15, 1};
Line(25) = {23, 11};
Line(26) = {11, 6};
Line(27) = {20, 16};
Line(28) = {16, 2};
Line(29) = {24, 12};
Line(30) = {12, 7};
Line(31) = {21, 17};
Line(32) = {17, 3};
Line(33) = {25, 10};
Line(34) = {10, 8};
Circle(35) = {10, 9, 14};
Circle(36) = {14, 9, 13};
Circle(37) = {13, 9, 15};
Circle(38) = {15, 9, 11};
Circle(39) = {11, 9, 16};
Circle(40) = {16, 9, 12};
Circle(41) = {12, 9, 17};
Circle(42) = {17, 9, 10};
Line Loop(43) = {8, -20, -35, 34};
Plane Surface(44) = {43};
Line Loop(45) = {9, -22, -36, 20};
Plane Surface(46) = {45};
Line Loop(47) = {10, -24, -37, 22};
Plane Surface(48) = {47};
Line Loop(49) = {3, -26, -38, 24};
Plane Surface(50) = {49};
Line Loop(51) = {4, -28, -39, 26};
Plane Surface(52) = {51};
Line Loop(53) = {5, -30, -40, 28};
Plane Surface(54) = {53};
Line Loop(55) = {6, -32, -41, 30};
Plane Surface(56) = {55};
Line Loop(57) = {7, -34, -42, 32};
Plane Surface(58) = {57};
Line Loop(59) = {35, -19, -18, 33};
Plane Surface(60) = {59};
Line Loop(61) = {36, -21, -11, 19};
Plane Surface(62) = {61};
Line Loop(63) = {37, -23, -12, 21};
Plane Surface(64) = {63};
Line Loop(65) = {38, -25, -13, 23};
Plane Surface(66) = {65};
Line Loop(67) = {39, -27, -14, 25};
Plane Surface(68) = {67};
Line Loop(69) = {40, -29, -15, 27};
Plane Surface(70) = {69};
Line Loop(71) = {41, -31, -16, 29};
Plane Surface(72) = {71};
Line Loop(73) = {42, -33, -17, 31};
Plane Surface(74) = {73};
Line(77) = {9, 22};
Line(78) = {9, 25};
Line(79) = {9, 23};
Line(80) = {9, 24};
Line Loop(81) = {18, 11, -77, 78};
Plane Surface(82) = {81};
Line Loop(83) = {17, -78, 80, 16};
Plane Surface(84) = {83};
Line Loop(85) = {12, 13, -79, 77};
Plane Surface(86) = {85};
Line Loop(87) = {14, 15, -80, 79};
Plane Surface(88) = {87};
Physical Surface(2) = {82, 86, 88, 84, 60, 62, 64, 66, 68, 70, 72, 74};
Physical Surface(1) = {44, 46, 48, 50, 52, 54, 56, 58};
Transfinite Line { 9, 10, 3, 4, 5, 6, 7, 8, 35, 36, 37, 38, 39, 40, 41, 42, 18, 11, 12, 13, 14, 15, 16, 17, 78, 77, 79, 80} = 3 Using Progression 1;
Transfinite Line {33, 19, 21, 23, 25, 27, 29, 31} = 4 Using Progression 1;
Transfinite Line {20, 22, 24, 26, 28, 30, 32, 34} = 3 Using Progression 1;
Mesh.ElementOrder = 2;

Mesh.Points = 1;
Mesh.SurfaceFaces=1;
Mesh.SurfaceEdges=1;
Mesh.VolumeEdges=1;
Mesh.ColorCarousel=2;
Mesh.Points = 1;
Mesh.SurfaceFaces=1;
Mesh.SurfaceEdges=1;
Mesh.VolumeEdges=1;
Mesh.ColorCarousel=2;
//+
Physical Curve(1) = {9, 10};
//+
Physical Curve(2) = {8, 7};
//+
Physical Curve(3) = {6, 5};
//+
Physical Curve(4) = {3, 4};

bagla0 · March 18, 2023, 6:56pm

Thanks @dokken for your response.

I used the meshio as mentioned by you for my quad cells mesh. I am able to generate the xdmf files using line_mesh and quad_mesh commands.

import meshio
mesh_from_file = meshio.read("vf_60_quad.msh")
mesh_from_file

<meshio mesh object>
  Number of points: 209
  Number of cells:
    line: 32
    quad: 192
  Cell data: gmsh:physical, gmsh:geometrical

import numpy
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("facet_mesh.xdmf", line_mesh)

quad_mesh = create_mesh(mesh_from_file, "quad", prune_z=True)
meshio.write("mesh.xdmf", quad_mesh)

I am finding difficulty to use these data for the FEniCS input as it was done in xml file using method Mesh and MeshFunction.

I am again wanting to assure that, using meshio, we can get workable quad mesh for FEniCS/Dolfin.

Any help is greatly appreciated.

dokken · March 18, 2023, 7:07pm

As I said in the previous post

In DOLFNx, you can load these meshes with no problem.
I believe that you cannot use them with legacy DOLFIN (as I again said, the support isn’t there).

There are many post on how to read XDMFFiles into DOLFIN:
For instance:

Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio

from dolfin import * 
mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2) 
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
File("dolfinmesh.pvd").write(mesh)
File("dolfincellfunc.pvd").write(mf)

bagla0 · March 18, 2023, 7:36pm

Thanks @dokken for your response.

What I know about dolfinx is, it is part of FEniCSx with other updated libraries like FFC-x and I need to install FEniCSx seperately.
Moreover, FEniCSx has different syntax for similar required work than FEniCS.
I am sorry, but, it feels little confusing to me that,
Dolfin legacy is related to FEniCS where quad elements are not supported but,
FEniCSx can support quad elements using Dolfinx legacy (Is it what you want to say about Dolfinx ?)
(and if I want to use quad element .msh file (simply like triangle cells), I need to use FEniCSx.)
I also heard about firedrake which is related to FEniCS, Can you also give your comments on this?

dokken · March 18, 2023, 7:56pm

DOLFINx has a different API, which is more verbose, but more efficient than Legacy DOLFIN I most cases.
DOLFIN is not being actively maintained and updated with new features.

DOLFIN (legacy): Some initial support for quads and hexes were added to legacy dolfin, for instance with built in meshes. External meshes (coming form msh, xdmf and other formats were never fully supported). There is an open issue at: fenics-project / DOLFIN / issues / #997 - Quad/hex meshes need ordering check — Bitbucket which probably never will be resolved.
DOLFINx: Support for quadrilaterals throughout the code (even higher order quadrilateral cells). Thus to use an external mesh with any FEniCS code, you would need to use DOLFINx.

Firedrake spun out of FEniCS a long while back (and has a quite similar interface to FEniCS). It is fully written in Python. I am not going to add any comments regarding their functionality, as that is best asked to the firedrake developers (and I am biased as DOLFIN/DOLFINx developer), see: https://www.firedrakeproject.org/ for information about that project.

Topic		Replies	Views
Physical region and facet not read by meshconvert mesh	12	799	August 4, 2021
Convert quad element 2D at .msh in .xml mesh	1	425	August 3, 2021
Dolfin-convert for convert msh to xml file mesh	1	1140	May 3, 2020
Conversion from .msh to .xml	6	5435	May 13, 2019
Problem converting file.msh (from gmsh) to file.xml I/O	12	4382	February 5, 2019

Msh to xml for quad elements

Related topics