How to import mixed triangle and quad meshes from gmsh to fenics

zhuwenguo · June 18, 2020, 10:44am

Hello,
I need to import meshes generated by gmsh to fenics, the mesh consists both quadrilateral and triangle elements. is it possible to use meshio to convert such a mesh into XDMF format?

thanks in advance!
wenguo.

dokken · June 18, 2020, 10:58am

There is no Mixed element support in fenics. The aim is to support this in the future in dolfinx.

Note that the support for quarilateral/hexahedral elements in dolfin is quite limited, and i strongly advice you to move to dolfinx if you require such elements.

Using the following arguments in gmsh will force all cells to be quadrilateral:

Recombine Surface {:};
Mesh.RecombinationAlgorithm = 2;

zhuwenguo · June 18, 2020, 11:25am

Hi, dokken
Thank you for your reply. The reason why I need to use the quad mesh is that I use the quad mesh within the region of layers where is the fabric composite materials. So I need the orientation angle of the quad element. After import the mesh into dolfin, I am planning to use the mesheditor to reconstructe the mesh which split the quad element into triangle elements. and since the geometry is not regular so it is not a good choice to use quad mesh only.
If meshio can conver mixed mesh into xdmf format, can solve my problem.
wenguo.

dokken · June 18, 2020, 11:41am

Meshio can extract mesh data, and you can modify it as you please, but the modifications has to be done before loading the mesh into dolfin via xdmf.

zhuwenguo · June 18, 2020, 11:56am

Hi dokken,
What does you mean “modifications”? you mean modify the meshio source code or made some changes on this? https://fenicsproject.discourse.group/t/transitioning-from-mesh-xml-to-mesh-xdmf-from-dolfin-convert-to-meshio/412

dokken · June 18, 2020, 12:06pm

As you would like to split all of your quadrilateral cells into triangles, you have to do this after reading the msh file into Python using meshio. Then you can get your quarilateral cells, and split them into triangles before writing them to xdmf.

zhuwenguo · June 18, 2020, 5:28pm

Hi,dokken
You’re right, I need to split the quad before convert to xdmf format. but it needs to modify the Mesh object. I think it is better to use meshio to read the mesh form gmsh .msh file and then with all the information stored in meshio::Mesh object, I can use dolfin.MeshEditor() to construct the whole mesh without converting into XDMF. and also split the quad at this stage.
And I have a problem about meshio’s data structure, what does the “field_data”, “cell_sets” mean?
I generated the following mixed mesh with python code, in the mesh there are two physical surface “quad” and “triangle” corresponding to the upper rectangle and lower rectangle respectively ,and two physical curve “top” and “bottom” correspoing to top line and bottom line respectively:

import gmsh
gmsh.initialize()
gmsh.option.setNumber(“General.Terminal”, 1)
gmsh.model.add(“demo”)

lc = 0.2
gmsh.model.geo.addPoint(-0.5, -1, 0, lc, 1)
gmsh.model.geo.addPoint(0.5, -1, 0, lc, 2)
gmsh.model.geo.addPoint(-0.5, 0, 0, lc, 3)
gmsh.model.geo.addPoint(0.5, 0, 0, lc, 4)
gmsh.model.geo.addPoint(-0.5, 1, 0, lc, 5)
gmsh.model.geo.addPoint(0.5, 1, 0, lc, 6)

gmsh.model.geo.addLine(1, 2, 1)
gmsh.model.geo.addLine(2, 4, 2)
gmsh.model.geo.addLine(4, 3, 3)
gmsh.model.geo.addLine(1, 3, 4)
gmsh.model.geo.addLine(3, 5, 5)
gmsh.model.geo.addLine(5, 6, 6)
gmsh.model.geo.addLine(4, 6, 7)

gmsh.model.geo.addCurveLoop([1, 2, 3, -4], 1)
gmsh.model.geo.addCurveLoop([-3, 7, -6, -5], 2)
gmsh.model.geo.addPlaneSurface([1], 1)
gmsh.model.geo.addPlaneSurface([2], 2)
gmsh.model.geo.mesh.setTransfiniteCurve(1, 3)
gmsh.model.geo.mesh.setTransfiniteCurve(2, 3)
gmsh.model.geo.mesh.setTransfiniteCurve(3, 3)
gmsh.model.geo.mesh.setTransfiniteCurve(4, 3)
gmsh.model.geo.mesh.setTransfiniteCurve(5, 3)
gmsh.model.geo.mesh.setTransfiniteCurve(6, 3)
gmsh.model.geo.mesh.setTransfiniteCurve(7, 3)
gmsh.model.geo.mesh.setTransfiniteSurface(1)
gmsh.model.geo.mesh.setRecombine(2, 1)

gmsh.model.addPhysicalGroup(2, [1], 1)
gmsh.model.setPhysicalName(2, 1, “quad”)

gmsh.model.addPhysicalGroup(2, [2], 2)
gmsh.model.setPhysicalName(2, 2, “triangle”)

gmsh.model.addPhysicalGroup(1, [6], 3)
gmsh.model.setPhysicalName(1, 3, “top”)

gmsh.model.addPhysicalGroup(1, [1], 4)
gmsh.model.setPhysicalName(1, 4, “bottom”)

gmsh.model.geo.synchronize()
gmsh.model.mesh.generate(2)
gmsh.write(“demo.msh”)
gmsh.fltk.run()
gmsh.finalize()

and then use meshio to get the mesh data, with dict

import meshio
from dolfin import *
msh = meshio.read("demo.msh")
print(type(msh))
for cell in msh.cells:
    if cell.type == "triangle":
        triangle_cells = cell.data
    elif cell.type == "quad":
        quad_cells = cell.data
    elif  cell.type == "line":
        line_cells = cell.data
    elif  cell.type == "tetra":
        tetra_cells = cell.data

print(msh.__dict__)

I get (I have ommited the output for “points”, “cells”, “point_Data”, “cell_data”)

‘field_data’: {‘top’: array([3, 1]), ‘bottom’: array([4, 1]), ‘quad’: array([1, 2]), ‘triangle’: array([2, 2])}, ‘point_sets’: {}, ‘cell_sets’: {‘top’: [array(, dtype=uint64), array([0, 1], dtype=uint64), array(, dtype=uint64), array(, dtype=uint64)], ‘bottom’: [array([0, 1], dtype=uint64), array(, dtype=uint64), array(, dtype=uint64), array(, dtype=uint64)], ‘quad’: [array(, dtype=uint64), array(, dtype=uint64), array([0, 1, 2, 3], dtype=uint64), array(, dtype=uint64)], ‘triangle’: [array(, dtype=uint64), array(, dtype=uint64), array(, dtype=uint64), array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15],
dtype=uint64)]}, ‘gmsh_periodic’: None, ‘info’: None}

demo1168×915 70.8 KB

dokken · June 18, 2020, 6:38pm

Your guess is as good as mine. The only things required to create a mesh in dolfin is:

The points of the mesh (msh.points)
The connectivity if each cell (msh.cells)
The physical data tags (msh.cell_data)

See for instance: Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio

zhuwenguo · June 19, 2020, 11:07am

Hi dokken,
After I transformed the mesh into dolfin, I need to apply a DirichletBC to a set of nodes which have been grouped into a physical group in Gmsh, a subdomain is needed for the DirichletBC. so my problem is how to create a subdomain based on the mesh information I have got in meshio without converting .msh to xdmf?

DirichletBC(U, Constant(0.0), subdomain, method='pointwise')

Thank you.
wenguo

dokken · June 19, 2020, 11:15am

You need to create a MeshValueCollection for the entities of interest.

mesh = UnitSquareMesh(10,10)
facet_collection = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
facet_collection.set_value(a, b) # Can use set_values as well
facet_function = cpp.mesh.MeshFunctionSizet(mesh, facet_collection)

zhuwenguo · June 19, 2020, 6:47pm

Hi, dokken
I am sorry, I still confused about the reason to create a MeshValueCollection(“size_t”, mesh, mesh.topology().dim()-1), For the demo above, I assign the top and bottom edges of the geometry as physical curve. and I want to apply pointwise Dirichlet boundary condition on the nodes on these two edges. according to the construction of DirichletBC, it seem that a subdomain is needed, right?

And, The physical curve information I got from meshio “cells” is as following:

'cells': [CellBlock(type='line', data=array([[0, 6],[6, 1]])), CellBlock(type='line', data=array([[ 4, 11],[11, 5]]))

It gives the nodes index on the boundary, if I am not wrong MeshValueCollection(“size_t”, mesh, mesh.topology().dim()-1) for 2d problem it is for 1-d entities, how to assign nodes values on 1-d mvc?

dokken · June 19, 2020, 10:03pm

You can use mesh functions, instead of subdomains in a Dirichlet condition to apply boundary conditions. I do not see the Need for a pointwise condition if you know the facets you have dofs on where you would like to assign the boundary condition. See for instance:
https://bitbucket.org/fenics-project/dolfin/raw/f0a90b8fffc958ed612484cd3465e59994b695f9/python/demo/documented/stokes-taylor-hood/demo_stokes-taylor-hood.py
Where subdomains is a mesh function that is employed in the Dirichlet condition

zhuwenguo · June 20, 2020, 7:48am

Hi dokken,

Since in my problem, I need to solve a linear system like this:

 bc.apply(A, rhs)
 solve(A, U.vector(), rhs)

So I need to set the DirichletBC method as “pointwise”.

dokken · June 20, 2020, 7:55am

You can still use the constructer i was pointing to, changing the method to pointwise

zhuwenguo · June 20, 2020, 8:14am

Hi, dokken
In my 2d problem, I used a mixedfunction, I need two kind of pointwise DirichletBC:

for some reason, I need to set the subfunction of some region to be zeros, so I need to construct a 2d meshfunction.
to assign known values to the part of the boundary.

UF3_ELEMENT = VectorElement("CG", mesh.ufl_cell(), degree, 3)
UF3 = FunctionSpace(mesh, UF3_ELEMENT)

VF1_ELEMENT = FiniteElement("CG", mesh.ufl_cell(), degree=degree)
VF1 = FunctionSpace(mesh, VF1_ELEMENT)

UF3VF1_ELEMENT = MixedElement(UF3_ELEMENT, VF1_ELEMENT)
UF3VF1 = FunctionSpace(mesh, UF3VF1_ELEMENT)

volumes = MeshFunction("size_t", mesh, mesh.topology().dim())
bc = DirichletBC(UF3VF1.sub(1), Constant(0.0), volumes, 1, method='pointwise'))
bc.apply(A, rhs)
solve(A, U.vector(), rhs)

But I get the following error:

*** Error:   Unable to create Dirichlet boundary condition.
*** Reason:  User MeshFunction is not a facet MeshFunction (dimension is wrong).
*** Where:   This error was encountered inside DirichletBC.cpp.

dokken · June 20, 2020, 11:10am

Since you are trying to deactivate all the dofs in a given cell, you need to map these regions to a mesh function of facets.
However, for me to be of proper help it would be beneficial if you could create a minimal example such a that I can get a proper idea of what you want to achieve.
as far as I can understand, you would like to do something like this:

from dolfin import *

mesh = UnitSquareMesh(10,10)
V = FunctionSpace(mesh, "CG", 2)
v = Function(V)
cf = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
cf.array()[:5] = 1

ff = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
for cell in cells(mesh):
    if cf.array()[cell.index()] == 1:
        facets = cell.entities(mesh.topology().dim()-1)
        for facet in facets:
            ff.array()[facet] = 1

File("cf.pvd").write(cf)
bc = DirichletBC(V, Constant(2), ff, 1)
bc.apply(v.vector())
XDMFFile("v.xdmf").write_checkpoint(v, "v", 0)

If you have a look at cf.xdmf and v.xdmf, you will observe that the dofs corresponding to each cell has been marked by the boundary condition

zhuwenguo · June 20, 2020, 12:49pm

Hi dokken,

I am sorry. Here I create a minimal demo. For a geometry combined by two squares vertically. I want to assign two Dirichlet boundary conditions.

The subfunction UV.sub(1) of the mixed function UV on the upper square to be zero;
At the same time assign 1.0 and -1.0 to the top and bottom edges of the lower square’s UV.sub(1).

demo02652×842 80 KB

The mesh file “demo.msh” is generated by:

lc = 0.5;
Point(1) = {-0.5, -1, 0, lc};
Point(2) = {0.5, -1, 0, lc};
Point(3) = {-0.5, 0, 0, lc};
Point(4) = {0.5, 0, 0, lc};
Point(5) = {-0.5, 1, 0, lc};
Point(6) = {0.5, 1, 0, lc};
Line(1) = {1, 2};
Line(2) = {2, 4};
Line(3) = {4, 3};
Line(4) = {1, 3};
Line(5) = {3, 5};
Line(6) = {5, 6};
Line(7) = {4, 6};
Curve Loop(1) = {4, -3, -2, -1};
Plane Surface(1) = {1};

Curve Loop(2) = {5, 6, -7, 3};
Plane Surface(2) = {2};

Physical Surface("passive") = {2};
Physical Surface("active") = {1};

Physical Curve("active_upper") = {3};
Physical Curve("active_bottom") = {1};

Then I imported the mesh using meshio, and reconstructed the mesh using mesheditor of dolfin. In order to apply the boundary condition I create 2 meshfunctions volumes and top_bottom_edge.

import numpy as np
import meshio
from dolfin import *

msh = meshio.read("demo.msh")

# cells from more than one cellBlocks
triangle_cells = np.array([None])
line_cells     = np.array([None])
triangle_data = np.array([None])
line_data     = np.array([None])

for cell in msh.cells:
    if cell.type == "triangle":
        if triangle_cells.all() == None:
            triangle_cells = cell.data
        else:
            triangle_cells = np.concatenate((triangle_cells,cell.data))

    elif  cell.type == "line":
        if line_cells.all() == None:
            line_cells = cell.data
        else:
            line_cells = np.concatenate((line_cells,cell.data))
    else:
        break

for key in msh.cell_data_dict["gmsh:physical"].keys():
    if key == "triangle":
        if triangle_data.all() == None:
            triangle_data = msh.cell_data_dict["gmsh:physical"][key]
        else:
            triangle_data = np.concatenate((triangle_data,msh.cell_data_dict["gmsh:physical"][key]))

    elif  key == "line":
        if line_data.all() == None:
            line_data = msh.cell_data_dict["gmsh:physical"][key]
        else:
            line_data = np.concatenate((line_data,msh.cell_data_dict["gmsh:physical"][key]))
    else:
        break

# reconstruct mesh using MeshEditor.
points_num   = msh.points.shape[0]
triangle_num = triangle_data.shape[0]

dolfin_mesh = Mesh()
me = MeshEditor()
me.open(dolfin_mesh, "triangle", 2, 2)
me.init_vertices(points_num)
for n in range(points_num):
    me.add_vertex(n, [msh.points[n,0], msh.points[n,1]])

me.init_cells(triangle_num)
for c in range(triangle_num):
    me.add_cell(c, [triangle_cells[c,0], triangle_cells[c,1], triangle_cells[c,2]])

me.close()

# Define function

U_ELEMENT = FiniteElement("CG", dolfin_mesh.ufl_cell(), 1)
U = FunctionSpace(dolfin_mesh, U_ELEMENT)
V_ELEMENT = FiniteElement("CG", dolfin_mesh.ufl_cell(), 1)
V = FunctionSpace(dolfin_mesh, V_ELEMENT)
UV_ELEMENT = MixedElement(U_ELEMENT, V_ELEMENT)
UV = FunctionSpace(dolfin_mesh, UV_ELEMENT)
uv = Function(UV)

# 1. Define meshfunction to assign the subfunction v of the upper squre(whose volumes[upper_squre] = 1) to be zero.

volumes = MeshFunction("size_t", dolfin_mesh, dolfin_mesh.topology().dim())
volumes.set_all(0)
for cell in cells(dolfin_mesh):
    i = cell.index()
    volumes[i] = triangle_data[i]

bc = DirichletBC(UV.sub(1), Constant(0.0), volumes, 1, method='pointwise')
bc.apply(uv.vector())

v = Function(V, name = "v")
assign(v, uv.sub(1))
volumesFile = XDMFFile('volumes.xdmf')
volumesFile.write(v)

# 2. Define meshfunction to assign the subfunction v of the top edge of lower squre to be 1.0, and the bottome to be -1.

top_bottom_edge = MeshFunction("size_t", dolfin_mesh, dolfin_mesh.topology().dim()-2)
top_bottom_edge.set_all(0)
for i in range(line_cells.shape[0]):
    p1 = line_cells[i,0]
    p2 = line_cells[i,1]

    top_bottom_edge.set_value(p1, line_data[i])
    top_bottom_edge.set_value(p2, line_data[i])
bcs = []
bcs.append(DirichletBC(UV.sub(1), Constant(1.0), top_bottom_edge, 1, method='pointwise'))
bcs.append(DirichletBC(UV.sub(1), Constant(-1.0), top_bottom_edge, 2, method='pointwise'))
[bc.apply(uv.vector()) for bc in bcs]

v = Function(V, name = "v")
assign(v, uv.sub(1))
top_bottom_edgeFile = XDMFFile('top_bottom_edge.xdmf')
top_bottom_edgeFile.write(v)

dokken · June 20, 2020, 12:51pm

@zhuwenguo can you confirm if what i suggested in the previous post could work for your problem ?

zhuwenguo · June 20, 2020, 12:58pm

Hi dokken,
Thank you for your help, still have trouble on apply boundary conditions to specific edges of cells. in your code

it seems that all the edges of these cells are assigned bc, but for my problem, i need to apply bc to those edges that are specified in the gmsh through physical curves.

Thank you for you kindness.

dokken · June 20, 2020, 1:33pm

Then you need to further filter out the facets in some way, for instance by checking if it is an exterior or interior facet, or checking geometrical properties of the facet. If the facets are assign through physical Curves in Gmsh, i do not see the reason for wanting to use a cell function in your Dirichlet condition. You should then create a meshfunction for facets, and add values to those corresponding to the info from Gmsh.

Topic		Replies	Views
Hybrid mesh for PML	14	1786	July 31, 2020
Importing XMDF in fenics for mesh I/O	5	577	September 29, 2022
Tetrahedral and Triangular Meshes Imported from Gmsh in FEniCSx dolfinx mesh , dolfinx , submesh	1	280	February 6, 2024
.xdmf mixed type - problem with generating mesh mesh	3	397	August 10, 2023
Gmsh, Meshio, Dolfin mesh	5	1416	June 3, 2023

How to import mixed triangle and quad meshes from gmsh to fenics

Related topics