Generation of second-order triangular meshes in legacy FEniCS

General
,
1

Hello everyone,

I am interested in using second-order triangular meshes in legacy FEniCS, but such an application seems to be only found in the UnitDiscMesh function. I have two questions regarding this:

  1. Is it possible to create a second-order triangular mesh using two numpy arrays of vertex and cell through MeshEditor()? I have found that this is implemented in UnitDiscMesh.cpp, but it appears that the Python version does not have a method like editor.add_entity_point() to add edge-based points. When I attempt to create a mesh using only editor.add_vertex() and editor.add_cell(), the resulting mesh seems invalid. The MWE is as follows:
import numpy as np
from fenics import *

tdim, gdim, degree = 2, 2, 2
mesh, editor = Mesh(), MeshEditor()
vertex = np.array([[0, 0], [1, 0], [0, 1],
                   [0.5, 0], [0, 0.5], [1, 1]])
cell = np.array([[0, 1, 2, 3, 4, 5]])

editor.open(mesh, "triangle", tdim, gdim, degree)
editor.init_vertices(len(vertex))
for i in range(len(vertex)):
    editor.add_vertex(i, vertex[i, :])
editor.init_cells(len(cell))
for i in range(len(cell)):
    editor.add_cell(i, cell[i, :])
editor.close()

with XDMFFile("text.xdmf") as outfile:
    outfile.write(mesh)
  1. Also, I can first generate a second-order triangular mesh using Gmsh with circle.geo and then convert circle.msh to circle.xdmf using meshio. However, the mesh data read through XDMFFile appears to be corrupted, and the difference between circle_out.xdmf and circle.xdmf can be seen in ParaView. The MWE is as follows:

This is my circle.geo file:

SetFactory("OpenCASCADE");

a = 0;
b = 0;
r = 1;
N = 6;

Circle(1) = {a, b, 0, r, 0, 2*Pi};
Curve Loop(1) = {1};
Plane Surface(1) = {1};
Physical Surface(1) = {1};

Transfinite Curve {1} = N+1 Using Progression 1;

This is my circle.msh file:

$MeshFormat
2.2 0 8
$EndMeshFormat
$Nodes
19
1 1 0 0
2 0.4999999999999988 0.8660254037844394 0
3 -0.5000000000000021 0.8660254037844374 0
4 -1 -1.20980699416795e-15 0
5 -0.4999999999999981 -0.8660254037844397 0
6 0.5000000000000009 -0.8660254037844382 0
7 0.8660254037844383 0.5000000000000007 0
8 -1.937171126634163e-15 1 0
9 -0.8660254037844396 0.4999999999999984 0
10 -0.8660254037844377 -0.5000000000000017 0
11 1.592665886326894e-15 -1 0
12 0.8660254037844388 -0.4999999999999997 0
13 -4.077392162158101e-16 -1.425974687551479e-16 0
14 0.4999999999999997 -7.129873437757396e-17 0
15 0.2500000000000003 -0.4330127018922191 0
16 0.2499999999999992 0.4330127018922196 0
17 -0.2499999999999993 -0.4330127018922199 0
18 -0.5000000000000002 -6.762022314615491e-16 0
19 -0.2500000000000012 0.4330127018922186 0
$EndNodes
$Elements
6
1 9 2 1 1 1 13 6 14 15 12
2 9 2 1 1 2 13 1 16 14 7
3 9 2 1 1 5 13 4 17 18 10
4 9 2 1 1 6 13 5 15 17 11
5 9 2 1 1 3 13 2 19 16 8
6 9 2 1 1 4 13 3 18 19 9
$EndElements

This is my Python code:

import meshio

def create_mesh(mesh, cell_type, data_name="name_to_read", 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={data_name:[cell_data]})
    return out_mesh

mesh_from_file = meshio.read("circle.msh")
mesh_triangle = create_mesh(mesh_from_file, "triangle6", prune_z=True)
meshio.write("circle.xdmf", mesh_triangle)

from fenics import *

mesh = Mesh()
with XDMFFile("circle.xdmf") as infile:
    infile.read(mesh)

with XDMFFile("circle_out.xdmf") as outfile:
    outfile.write(mesh)

Alternatively, are there other ways to use second-order triangular meshes in legacy FEniCS besides UnitDiscMesh? I am aware that FEniCSx has advantages in supporting higher-order meshes, but due to certain historical reasons, I am currently still using legacy FEniCS. I would greatly appreciate any assistance.

2

This question has been asked several times, see for instance Is it possible to program (from scratch) a 2nd-order tetrahedron in FeniCs Legacy?

IMHO, the short answer is: while it is certainly possible to support to some extent higher order meshes in legacy dolfin, the most likely scenario is that you will have to implement it on your own.

3

Thank you very much for your response.

I have reviewed this answer as well as similar discussions previously, but I still find myself unclear about certain implementation details in Is it possible to program (from scratch) a 2nd-order tetrahedron in FeniCs Legacy?. Specifically, I am unsure how to implement the I/O interface following the lines of dolfinx.io.gmshio and dolfinx.plot.vtk_mesh to read higher-order meshes in legacy FEniCS. Additionally, I am having some difficulty understanding the statement “ordered in the UFC fashion, not VTK or Gmsh node ordering”. Could you please provide further clarification on this matter?

I would also like to inquire if it is possible to accomplish this with concise Python code instead of C++ code, and if there are more applications of higher-order meshes in legacy FEniCS that might be relevant to my inquiry, such as Support for curved edges.

I would greatly appreciate any further guidance or insights you might be able to provide. :slightly_smiling_face:

4

For a higher order grid you have to have some specific order you traverse node. For instance, for a higher order (2nd order triangle) degrees of freedom appears on the faces of each element. You need to loop through these in a consistent order.
VTK, GMSH and DOLFINx operate with different ways of doing this.
See for instance

image
image180×180 4.45 KB

From
DefElement Defilement that shows the dolfinx ordering, while VTK uses
image

from https://www.kitware.com/modeling-arbitrary-order-lagrange-finite-elements-in-the-visualization-toolkit/
perm_vtk in dolfinx maps the transformation from one to the other.

This has to go through mesheditor, as it is the only way of constructing meshes in legacy dolfin.

I would strongly suggest to move to dolfinx, as you would have to add these feature to legacy dolfin yourself (as we as volunteers don’t have time to maintain and update a deprecated software).

There has in addition been made several other improvements to dolfinx, supporting arbitrary order meshes (legacy dolfin is restricted to second order).

Also the code compilation of such forms are way better in dolfinx.

5

Thank you so much for your sincere response! I plan to complete my legacy FEniCS code as much as possible before fully transitioning to FEniCSx. I don’t have any further questions on this matter for now. Wish you a wonderful day!

6

I believe I have finally achieved the customized generation of second-order triangular meshes with the limited support of legacy FEniCS, thanks to the hidden built-in function p_refine as used in ringleb_example.py and ringleb_mesh.py. This function appears to be specifically useful for converting straight-edged triangular meshes into second-order ones, which is exactly what I needed.

To assist those who may require this, I would like to present a brief example below.

from fenics import *

n, degree, gdim = 1, 1, 2

mesh0 = UnitDiscMesh.create(MPI.comm_world, n, degree, gdim)
mesh = p_refine(mesh0)
mesh.init(mesh.topology().dim()-1, mesh.topology().dim())
for i in range(mesh.num_facets()):
    if Facet(mesh, i).exterior():
        midpoint = Facet(mesh, i).midpoint().array()[:2]
        mesh.coordinates()[mesh.num_vertices() + i] = midpoint/sqrt(sum(midpoint**2))

# which is equivalent to the following code
mesh = UnitDiscMesh.create(MPI.comm_world, n, degree+1, gdim)

I hope this information is helpful to some people. Once again, my heartfelt thanks to everyone!

7

But this only makes a second order mesh if you have a circular mesh (as it assumes that the new coordinate should be at radius=1).

It would not work for a generalized mesh, as you don’t know how to curve the midpoint of the facet.

Anyhow, if this is sufficient for your purpose, good for you :slight_smile: