Unable to read .msh mesh file which contains only rectangular elements

Hi all,
I have created a 2d mesh in GMSH with only rectangular elements. I want to read it in FEniCS.
I am following the code from https://jsdokken.com/converted_files/tutorial_pygmsh.html
The mesh has only line elements, no triangular and no tetrahedral elements.
But when I plot the mesh in FEniCS, only boundary cells are appearing, interior cells have vanished.
I am not understanding why is this happening.
I am attaching the code snippet and snapshot of the mesh.

from fenics import *
import meshio
mesh_from_file = meshio.read("rectangle.msh")

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)
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    if prune_z:
        out_mesh.prune_z_0()
    return out_mesh

line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)
meshio.write("facet_mesh.xdmf", line_mesh)

#triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True)
#meshio.write("mesh.xdmf", triangle_mesh)

mesh = Mesh()
with XDMFFile("facet_mesh.xdmf") as infile:
    infile.read(mesh)
    mvc = MeshValueCollection("size_t", mesh, 2)

plot(mesh)

You are stating the your mesh has rectangular elements, therefore, you should save a mesh with “quad” elements, not line elements. The line elements are only the facet markers that you had created in your mesh. Please note that support for quadrilateral elements are very limited in dolfin, and it is recommended to use dolfinx if you want to use such elements.
I have created a tutorial for dolfinx at: The FEniCSx tutorial — FEniCSx tutorial

1 Like

I have tried with ‘quad’ and ‘quadrilateral’ element, but getting this error.

line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)
meshio.write("facet_mesh.xdmf", line_mesh)

#triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True)
#meshio.write("mesh.xdmf", triangle_mesh)

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

ValueError                                Traceback (most recent call last)
<ipython-input-16-b3b44afb1659> in <module>
      5 #meshio.write("mesh.xdmf", triangle_mesh)
      6 
----> 7 quad_mesh = create_mesh(mesh_from_file, "quadrilateral", prune_z=True)
      8 meshio.write("mesh.xdmf", quad_mesh)

<ipython-input-15-90f22cb4147e> in create_mesh(mesh, cell_type, prune_z)
      2 def create_mesh(mesh, cell_type, prune_z=False):
      3     cells = mesh.get_cells_type(cell_type)
----> 4     cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
      5     out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
      6     if prune_z:

~/.local/lib/python3.6/site-packages/meshio/_mesh.py in get_cell_data(self, name, cell_type)
    226     def get_cell_data(self, name: str, cell_type: str):
    227         return np.concatenate(
--> 228             [d for c, d in zip(self.cells, self.cell_data[name]) if c.type == cell_type]
    229         )
    230 

ValueError: need at least one array to concatenate

When I type ‘‘mesh_from_file.cells’’, I am getting following o/p.

mesh_from_file.cells

[<meshio CellBlock, type: line, num cells: 19>,
 <meshio CellBlock, type: line, num cells: 39>,
 <meshio CellBlock, type: line, num cells: 19>,
 <meshio CellBlock, type: line, num cells: 39>]

Actually, there are 39*19= 741 cells in my mesh. Why here it is showing less, also it is showing only line elements, no quadrilateral elements.

This is probably because you have not specified a physical surface in you Gmsh script. Without the geo file i cant do much more.

I think, I had specified a physical surface in GMSH script.
It is not giving me permission to upload .geo file and .msh file.

You can add it by encapsulating it with 3x `
Such as shown in Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio - #62 by dokken

Okay. I am trying here. Please tell me if you are getting from this.

SetFactory("OpenCASCADE");
//+
Point(1) = {0, -0.1, 0, 1.0};
//+
Point(2) = {0, 1, 0, 1.0};
//+
Point(3) = {2, 1, 0, 1.0};
//+
Point(4) = {2, 0, 0, 1.0};
//+
Recursive Delete {
  Point{1}; 
}
//+
Point(5) = {0, 0, 0, 1.0};
//+
Line(1) = {2, 5};
//+
Line(2) = {5, 4};
//+
Line(3) = {4, 3};
//+
Line(4) = {3, 2};
//+
Curve Loop(1) = {4, 1, 2, 3};
//+
Plane Surface(1) = {1};
//+
Physical Curve("inlet", 5) = {1};
//+
Physical Curve("outlet", 6) = {3};
//+
Physical Curve("top", 7) = {4};
//+
Physical Curve("bottom", 8) = {2};
//+
Transfinite Surface {1} = {2, 3, 4, 5};
//+
Transfinite Curve {1, 3} = 20 Using Progression 1;
//+
Transfinite Curve {4, 2} = 40 Using Progression 1;
//+
Recombine Surface {1};

By running gmsh -2 on this file and then opening the resulting msh file in gmsh, you will observe that there are no volume elements in the mesh.
With the following modifications:

SetFactory("OpenCASCADE");
//+
Point(1) = {0, -0.1, 0, 1.0};
//+
Point(2) = {0, 1, 0, 1.0};
//+
Point(3) = {2, 1, 0, 1.0};
//+
Point(4) = {2, 0, 0, 1.0};
//+
Recursive Delete {
  Point{1}; 
}
//+
Point(5) = {0, 0, 0, 1.0};
//+
Line(1) = {2, 5};
//+
Line(2) = {5, 4};
//+
Line(3) = {4, 3};
//+
Line(4) = {3, 2};
//+
Curve Loop(1) = {4, 1, 2, 3};
//+
Plane Surface(1) = {1};
//+
Transfinite Surface {1} = {2, 3, 4, 5};
//+
Transfinite Curve {1, 3} = 20 Using Progression 1;
//+
Transfinite Curve {4, 2} = 40 Using Progression 1;
//+
Recombine Surface {1};
//+
Physical Curve("inlet", 5) = {1};
//+
Physical Curve("outlet", 6) = {3};
//+
Physical Curve("top", 7) = {4};
//+
Physical Curve("bottom", 8) = {2};
//+
Physical Surface(9) = {1};

you can convert your mesh with the following commands

import numpy
import meshio
mesh_from_file = meshio.read("transf_gmsh.msh")

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)
    out_mesh = meshio.Mesh(points=mesh.points, cells={
                           cell_type: cells}, cell_data={"name_to_read": [cell_data]})
    if prune_z:
        out_mesh.prune_z_0()
    return out_mesh


line_mesh = create_mesh(mesh_from_file, "quad", prune_z=True)
meshio.write("mesh_quad.xdmf", line_mesh)
line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)
meshio.write("facet_mesh.xdmf", line_mesh)

but as I said, you should use FEniCSx if you want to work with quad meshes

2 Likes

Yes, sir, I will use FEniCSx when solving the actual problem.
Now I am only practicing.
I have used your file and the given code.
Now, I am getting this error.

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

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)
    out_mesh = meshio.Mesh(points=mesh.points, cells={
                           cell_type: cells}, cell_data={"name_to_read": [cell_data]})
    if prune_z:
        out_mesh.prune_z_0()
    return out_mesh

line_mesh = create_mesh(mesh_from_file, "quad", prune_z=True)
meshio.write("mesh_quad.xdmf", line_mesh)
line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)
meshio.write("facet_mesh.xdmf", line_mesh)

It is writing the file “mesh_quad.xdmf”
But when I am reading that mesh, I am getting following error.

mesh = Mesh()
with XDMFFile("mesh_quad.xdmf") as infile:
    infile.read(mesh)
    mvc = MeshValueCollection("size_t", mesh, 2)

RuntimeError                              Traceback (most recent call last)
<ipython-input-22-b5af187c9096> in <module>
      1 mesh = Mesh()
      2 with XDMFFile("mesh_quad.xdmf") as infile:
----> 3     infile.read(mesh)
      4     mvc = MeshValueCollection("size_t", mesh, 2)

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:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------

As I said in the previous posts, there is very limited support for quadrilaterals in dolfin, and this is why you cannot read in the xdmf file.

Okay, got it. The above issue is because of quadrilateral cells.
Now, I have created a mesh in a rectangular domain with purely triangular elements and I think I have successfully imported it into FEniCS too. After the plot(mesh) command it is showing the correct mesh.
I need to solve the problem for which I need to create subdomains for all four boundaries. I did the same thing with built-in UnitSquareMesh in FEniCS and got the correct answer. But when I am using the same code, it gives me an error while creating boundary parts for this mesh. I am attaching the code and mesh figure here.

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

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)
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    if prune_z:
        out_mesh.prune_z_0()
    return out_mesh

line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)
meshio.write("facet_mesh.xdmf", line_mesh)

triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True)
meshio.write("mesh.xdmf", triangle_mesh)

mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
    mvc = MeshValueCollection("size_t", mesh, 2)

# Define boundary segments for Neumann, Robin and Dirichlet conditions

# Create mesh function over cell facets
boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim()-1)

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-13-2980ce74ea92> in <module>
      2 
      3 # Create mesh function over cell facets
----> 4 boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim()-1)

AttributeError: module 'dolfin.mesh' has no attribute 'topology'

If there is any tutorial in which the mesh is imported and solved the problem by creating subdomains, please suggest me.

I have one more doubt, why do we create different meshes as facet_mesh and inner element mesh. What is the need for that?

You should really create your subdomains in the meshing tool, as you would otherwise get jagged interfaces, as the subdomains does not necessarily align with the elements if defined geometrically in dolfin.
For further help, please add the corresponding geo file (as I cannot reproduce any errors without the mesh.

The facet mesh is created such that we can read it into a mesh function (through a mesh value collection), so that we can define boundary conditions on these facets (without having to backward engineer these boundaries in dolfin).

1 Like

Ok. Thanks, sir for the help.
I am sending you the geometry file, mesh file, and .ipynb file on slack as it will be very much longer text here. Please tell me will it be ok for you otherwise, I will send it here.
The .ipynb file contains the code written for in-built UnitSquareMesh.
I want to write the same code and get the same result for the mesh created in GMSH. This is my requirement.

Please add the geo file here (encapsulated as previously). Similarly with the python code in the notebook. The reason for me requiring this is that this might help others with similar issues.

Here is the .geo file of my mesh.

// Gmsh project created on Wed Jun 16 11:22:49 2021
SetFactory("OpenCASCADE");
//+
Point(1) = {0, 0, 0, 3.0};
//+
Point(2) = {2, 0, 0, 3.0};
//+
Point(3) = {2, 1, 0, 3.0};
//+
Point(4) = {0, 1, 0, 3.0};
//+
Line(1) = {4, 1};
//+
Line(2) = {1, 2};
//+
Line(3) = {2, 3};
//+
Line(4) = {3, 4};
//+
Curve Loop(1) = {1, 2, 3, 4};
//+
Plane Surface(1) = {1};
//+
Physical Curve("fixed", 5) = {1};
//+
Physical Curve("free", 6) = {3};
//+
Physical Curve("top", 7) = {4};
//+
Physical Curve("bottom", 8) = {2};
//+
Physical Surface("fluid", 9) = {1};

This is the code which I have written for built-in Unitsquare mesh.

from fenics import *
import matplotlib.pyplot as plt

# Problem parameters
L = 1
F = 10
EA = 2
p = F/EA

# Defining mesh
mesh = UnitSquareMesh(2, 1)
V = FunctionSpace(mesh, 'P', 1)
plot(mesh)

# Define boundary segments for Neumann, Robin and Dirichlet conditions

# Create mesh function over cell facets
boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim()-1)

# Mark lower boundary facets as subdomain 0
class LowerNeumanBoundary(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14   # tolerance for coordinate comparisons
        return on_boundary and abs(x[1]) < tol
    
Gamma_R = LowerNeumanBoundary()
Gamma_R.mark(boundary_parts, 0)

# Mark upper boundary facets as subdomain 1
class UpperNeumannBoundary(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14   # tolerance for coordinate comparisons
        return on_boundary and abs(x[1] - 1) < tol

Gamma_N = UpperNeumannBoundary()
Gamma_N.mark(boundary_parts, 1)

# Mark left boundary as subdomain 2
class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14   # tolerance for coordinate comparisons
        return on_boundary and abs(x[0]) < tol

Gamma_0 = LeftBoundary()
Gamma_0.mark(boundary_parts, 2)

# Mark right boundary as subdomain 3
class RightBoundary(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14   # tolerance for coordinate comparisons
        return on_boundary and abs(x[0] - 1) < tol

Gamma_1 = RightBoundary()
Gamma_1.mark(boundary_parts, 3)

bc = DirichletBC(V, Constant(0), boundary_parts, 2)
ds = Measure("ds", subdomain_data=boundary_parts)

# Defining the problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0)
a = -inner(grad(u), grad(v))*dx
g = Expression('p', p=p, degree=1)
g1 = Constant(0)
L = f*v*dx - g*v*ds(3) + g1*v*ds(0) + g1*v*ds(1)

#problem = VariationalProblem(a, L, bc)
#u = problem.solve()

u = Function(V)
solve(a == L, u, bc)

plot(u)
u.vector().get_local()
array([ 0. ,  0. ,  2.5,  2.5,  5. ,  5. ])

This is the code that I want to write for the mesh created in GMSH, and want to get the similar result of what got in the above code.

from fenics import *
import matplotlib.pyplot as plt
pip install --user meshio

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

pip install --user h5py
import h5py
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)
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    if prune_z:
        out_mesh.prune_z_0()
    return out_mesh

line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)
meshio.write("facet_mesh.xdmf", line_mesh)

triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True)
meshio.write("mesh.xdmf", triangle_mesh)

mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
    mvc = MeshValueCollection("size_t", mesh, 2)

plot(mesh)
# Problem parameters
L = 2
F = 10
EA = 2
p = F/EA

# Define boundary segments for Neumann, Robin and Dirichlet conditions

# Create mesh function over cell facets
boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
AttributeError                            Traceback (most recent call last)
<ipython-input-13-2980ce74ea92> in <module>
      2 
      3 # Create mesh function over cell facets
----> 4 boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim()-1)

AttributeError: module 'dolfin.mesh' has no attribute 'topology'

I got the above error.

Running your code in a single cell works for me:

from fenics import *
import matplotlib.pyplot as plt

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

import h5py
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)
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    if prune_z:
        out_mesh.prune_z_0()
    return out_mesh

line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)
meshio.write("facet_mesh.xdmf", line_mesh)

triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True)
meshio.write("mesh.xdmf", triangle_mesh)

mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
    mvc = MeshValueCollection("size_t", mesh, 2)

plot(mesh)
# Problem parameters
L = 2
F = 10
EA = 2
p = F/EA

# Define boundary segments for Neumann, Robin and Dirichlet conditions
# Create mesh function over cell facets
boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim()-1)

Maybe it’s worth installing all dependencies and only then running the code in sequence, to make sure that all variables are in the correct state.

1 Like

Did you use the same mesh which I have given?
Did you not get the error for the last line

boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim()-1)```
If you have used different mesh, then I would request you please go through my .geo file once, and tell me if there is any mistake.
Thank you

Yes, I used your .geo file, and it works.

I have run my code again now and it worked fine. Maybe I had not installed all dependencies. Thanks for the help.

Hello,

I’m new to the forum and don’t know if this is the right place to ask my question, but I have a very similar problem regarding the reading of the .msh file. Following the tuto Mesh generation and conversion with GMSH and PYGMSH | Jørgen S. Dokken, I also have the error “ValueError: need at least one array to concatenate” (see picture below).

However, I did create a triangular mesh on GMSH.

Do you have any idea of the problem? Thanks a lot!

.geo file:

l_out = 1.0;
l_in = 0.1;
// Point of the domain
Point(1) = {0, 0, 0, l_out};
Point(2) = {10, 0, 0, l_out};
Point(3) = {10, 5-0.1, 0, l_in};
Point(4) = {10, 5+0.1, 0, l_in};
Point(5) = {10, 10, 0, l_out};
Point(6) = {0, 10, 0, l_out};
Point(7) = {0, 5+0.1, 0, l_in};
Point(8) = {0, 5-0.1, 0, l_in};
// Lines between points
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 5};
Line(5) = {5, 6};
Line(6) = {6, 7};
Line(7) = {7, 8};
Line(8) = {8, 1};
Line(9) = {8, 3};
Line(10) = {7, 4};
// Curve Loops
Curve Loop(1) = {1, 2, -9, 8}; // Line Loop
Curve Loop(2) = {4, 5, 6, 10};
Curve Loop(3) = {9, 3, -10, 7};
// Plane Surface
Plane Surface(1) = {1};
Plane Surface(2) = {2};
Plane Surface(3) = {3};
// Physical groups
// Lines
Physical Curve(“L_bottom”, 1) = {8, 1, 2};
Physical Curve(“L_top”, 2) = {4, 5, 6};
Physical Curve(“L_left”, 3) = {7};
Physical Curve(“L_right”, 4) = {3};
// Surfaces
Physical Surface(“Bottom”) = {1};
Physical Surface(“Top”) = {2};
Physical Surface(“Middle”) = {3};