Import second order mesh to FEnics using mesh.io

mesh
1

Hi all,

I want to import a second order mesh generated by GMSH into FEnics. I am using meshio to import the mesh. The gmsh .msh file* is in format of MeshFormat 2.2 0 8 . Since I want to use tetrahedron elements with 10 Nodes and triangle with 6 nodes, I try to write the code below.

I would be very thankful for any support :slight_smile:

Cheers
Ovel

*PS msh file below

msh = meshio.read(mesh_file) #select the mesh file

cell_element_order = "tetra10"
facet_element_order = "triangle6"
degree_FE = 1 # Finite Element degree

meshio.write("./Temp/mesh.xdmf", 
             meshio.Mesh(points=msh.points, 
                        cells={cell_element_order: msh.cells[cell_element_order]}, 
                        cell_data={cell_element_order: 
                                   {"name_to_read": msh.cell_data[cell_element_order]["gmsh:physical"]}}))

meshio.write("./Temp/mesh_boundary.xdmf", 
             meshio.Mesh(points=msh.points, 
                        cells={facet_element_order: msh.cells[facet_element_order]},
                        cell_data={facet_element_order: 
                                   {"name_to_read": msh.cell_data[facet_element_order]["gmsh:physical"]}}))
mesh = Mesh()

with XDMFFile("./Temp/mesh.xdmf") as infile:
    infile.read(mesh)
    mvc = MeshValueCollection("size_t", mesh, DIM) 
    
with XDMFFile("./Temp/mesh_boundary.xdmf") as infile:
    infile.read(mvc, "name_to_read")   
    mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

Unfortunately I got an error. Does anyone have any experience importing second order meshes?

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-14-5828977a4a77> in <module>
     19 
     20 with XDMFFile("./Temp/mesh.xdmf") as infile:
---> 21     infile.read(mesh)
     22     mvc = MeshValueCollection("size_t", mesh, DIM)
     23 

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 recognise cell type.
*** Reason:  Unknown value "tetrahedron_10".
*** Where:   This error was encountered inside XDMFFile.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  76bf12b24408e3ce2364eb22e62c802a92c90ff2
*** -------------------------------------------------------------------------

*Mesh file

$MeshFormat
2.2 0 8
$EndMeshFormat
$PhysicalNames
5
2 2 "ROOF_BACK_2D"
2 3 "FRONT_2D"
2 4 "LOUDSPEAKER_2D"
2 5 "VOLUMENK_RPER_1_1_2D"
3 1 "car_body_ansys_simplified_5cm_lo"
$EndPhysicalNames
$Nodes
90580
1 4.262164115905762 1.389446139335632 2.520000050645899e-012
2 4.227527141571045 1.391816854476929 0.1622721999883652
3 4.234670639038086 1.391465783119202 0.1213999465107918
4 4.243685722351074 1.39084804058075 0.08089913427829742
5 4.253204822540283 1.390129327774048 0.0405135490000248
6 4.194525241851807 1.372372269630432 0.5614960789680481
7 4.209649562835693 1.383093476295471 0.520927369594574
8 4.214374542236328 1.388261437416077 0.4765497148036957
9 4.215836048126221 1.389703273773193 0.431630551815033
10 4.215663433074951 1.390226721763611 0.3866600096225739
11 4.216094017028809 1.390782833099365 0.3416904211044312
12 4.21720552444458 1.391260027885437 0.29673171043396
....... * I cut some lines, because the file is quite long *
90578 2.592917442321777 0.6017194390296936 0.6876446008682251
90579 2.586379766464233 0.5731936693191528 0.6859179139137268
90580 2.834211587905884 0.6411230564117432 0.6843199133872986
$EndNodes
$Elements
71703
1 9 2 2 2 61 62 95 324 330 325
2 9 2 2 2 96 92 102 421 423 442
3 9 2 2 2 94 118 122 433 494 434
4 9 2 2 2 103 91 123 417 420 472
5 9 2 2 2 105 92 109 424 425 476
...
18062 9 2 5 5 2690 6278 4127 12824 19258 12823
18063 9 2 5 5 4127 6278 6803 19258 28293 19260
18064 9 2 5 5 4127 6803 4161 19260 19391 19257
18065 11 2 1 1 3324 10643 5843 40556 15558 26728 15555 43022 51032 90435
18066 11 2 1 1 3324 10122 10643 40556 43021 64400 15558 43022 90435 64398
18067 11 2 1 1 7895 10643 10641 10122 31028 36130 31027 56442 64399 64400
18068 11 2 1 1 4765 6939 37209 6948 21735 53935 48000 21736 53967 29485
18069 11 2 1 1 4187 37094 40718 8325 45676 73606 45680 19474 57834 57830
18070 11 2 1 1 4187 4188 40718 40665 19472 45683 45680 45679 90481 45682
...
71700 11 2 1 1 4772 6956 40895 40697 21770 54001 48029 48028 90489 54000
71701 11 2 1 1 4667 6898 38639 38592 21333 53843 47735 47734 83018 53841
71702 11 2 1 1 3844 36155 6660 6657 44644 53194 18033 18031 29030 53188
71703 11 2 1 1 39533 39581 39864 39943 88171 88404 88174 88175 89605 88405
$EndElements
2

There are no support for higher order tetrahedrons, quadrilaterals or hexahedrons in dolfin.
However, in dolfinx, there is support up to 9th order triangles and quadrilaterals, 3rd order tetrahedrons and 2nd order hexahedrons.
See for instance Higher order mesh tests.
However, right now, we do not support higher order mesh functions, but its on the todo list.

3

Hi Dokken, thank you for the information. It is very helpful.

Best regards,
Ovel

4

Hi all,
I have a Gmsh .msh file with a first order mesh. Then, using dolfin-convert, .xml files for mesh, facets and physical regions are created. I can import these files in FEniCS and my problem is solved correctly.

My aim is to solve the problem using higher order elements. As a first approach, I changed the finite element order from 1 to 2 in my FEniCs code, mantaining the same three .xml files for mesh, facets and subdomains. Any error is printed and it seems to work properly but when I visualize my results some problems appear in boundaries, seeming that second order nodes on boundaries are not being well identified in my code.

I would like to know if it is possible to follow this procedure to solve a problem with second order elements, without the need of explicitly read a second order mesh in FEniCS and only introducing this behaviour through the finite element order.

Thanks in advance.

5

Second order finite element functions are only properly visualized in Paraview if you use
XDMFFile.write_checkpoint. If you do not use this function, the solution will be interpolated to a first order space for visualization.

1 Like
6

Thank you very much. I was not visualizing my solution properly, now second order elements are observed. But my problem in that boundary still remains. There a Neumann condition is applied, through a Lagrange multiplier in my variational formulation, and in Paraview it seems that this boundary condition is the only one that is not being properly imposed.

Could it be a problem of mesh connectivity? As FEniCS is receiving a first order mesh and then the problem is solved with a second order element, perhaps the new boundary nodes are not automatically identified.

7

FEniCS can use any order for the finite element used to solve the problem, on a linear mesh, as the mesh element and the finite element is decoupled (as opposed to some traditional FEM software).

You can verify this for a manufactured solution, such as shown here:
https://jorgensd.github.io/dolfinx-tutorial/chapter4/convergence.html
(Using the development version of FEniCS),
and here:
https://fenicsproject.org/pub/tutorial/html/._ftut1020.html
(subsection computing convergence rates, using FEniCS v 2016.2.0)

If your boundary condition is not properly imposed, you need to create a minimal working code example that reproduces the problem (and that anyone can run on their local system and get the same result).

8

@ims Hi Ims,

Could you please let me know how did you convert the finite element order from 1 to 2. I am solving a solid mechanics problem and I need to use a 2nd order element. I have generated my mesh using Gmsh which is 1st order. After loading in fenics, I want to convert that mesh into 2nd order. If I create a 2nd order element mesh in Gmsh then I can not import that into Fenics.

9

Hi,

I you create a 1st order mesh in Gmsh (.msh format) then, by using dolfin-convert, you can generate the proper .xml files to read in FEniCS. If you have solved your problem with 1st order elements this proceeding will be familiar. So, in order to read mesh and subdomains in FEniCS:

mesh = Mesh(‘mesh.xml’)
subdomains = MeshFunction(‘size_t’, mesh, ‘mesh_physical_region.xml’)
boundaries = MeshFunction(‘size_t’, mesh, ‘mesh_facet_region.xml’)

Then, if you want to solve a problem with 1st order Lagrange elements, a function space is defined as follows:
V = FunctionSpace(mesh, “CG”, 1)

But if you want 2nd order Lagrange elements, simply increase the order of the function space as
V = FunctionSpace(mesh, “CG”, 2)
and the rest of the code will remain unchanged.

As far as I know, FEniCS automatically does the conversion of your input 1st order Gmsh mesh to solve the problem with 2nd order elements. I am not an expert, so I can be wrong, but I think this has worked for me.

1 Like
10

In FEniCS the mesh and function space are decoupled, as @ims suggests. This means that you can have an arbitrary order space on a first order mesh.

There is also some limited support for higher order meshes (curved facets), which has been extended in dolfinx.

1 Like