Mesh from a list of points and connectivity

mesh
1

Hi everyone,

I’m a beginner in FEniCSx, and I want to import a mesh that is stored in two lists, one with the points coordinates and one with the connectivity of the cells.

I managed to do that by first generating a .xmdf file using meshio and then importing the .xdmf file with dolfinx.io. Here is an example of what I do

import meshio
from dolfinx.io import XDMFFile
from mpi4py import MPI

points = [
    [-1.0, -1.0],
    [ 0.0, -1.0],
    [ 1.0, -1.0],
    [-1.0,  0.0],
    [ 0.0,  0.0],
    [ 1.0,  0.0],
    [-1.0,  1.0],
    [ 0.0,  1.0],
    [ 1.0,  1.0],
]

aux_cells = [
    [0,1,3,4],
    [1,2,4,5],
    [3,4,6,7],
    [4,5,7,8],
]

cells = [("quad", aux_cells)]

mesh = meshio.Mesh(points,cells)
meshio.write_points_cells("test.xdmf", points, cells)

with XDMFFile(MPI.COMM_WORLD, "test.xdmf", "r") as xdmf:
    msh = xdmf.read_mesh(name="Grid")

That works fine, but I wonder if there’s a more simple and direct approach.

2

If you want something that works in serial, see:

Consider the following if you want to work in parallel:

from mpi4py import MPI
import numpy as np
import dolfinx
import ufl

if MPI.COMM_WORLD.rank == 0:
    points = np.array([
        [-1.0, -1.0],
        [0.0, -1.0],
        [1.0, -1.0],
        [-1.0,  0.0],
        [0.0,  0.0],
        [1.0,  0.0],
        [-1.0,  1.0],
        [0.0,  1.0],
        [1.0,  1.0],
    ], dtype=np.float64)
    aux_cells = np.array([
        [0, 1, 3, 4],
        [1, 2, 4, 5],
        [3, 4, 6, 7],
        [4, 5, 7, 8],
    ], dtype=np.int64)

else:
    points = np.empty((0, 2), dtype=np.float64)
    aux_cells = np.empty((0, 4), dtype=np.int64)
domain = ufl.Mesh(ufl.VectorElement(
    "Lagrange", "quadrilateral", 1, dim=points.shape[1]))
mesh = dolfinx.mesh.create_mesh(MPI.COMM_WORLD, aux_cells, points, domain)
print(MPI.COMM_WORLD.rank, mesh.topology.index_map(mesh.topology.dim).size_local,
      mesh.topology.index_map(mesh.topology.dim).size_global)

returning:

root@dokken-XPS-15-9560:/fenics/shared# python3 mwe_hirsch.py 
0 4 4
root@dokken-XPS-15-9560:/fenics/shared# mpirun -n 2 python3 mwe_hirsch.py 
0 2 4
1 2 4
1 Like
3

Thanks Dokken. Just implemented the second way and it works great.

4

@dokken I think VectorElement is not defined anymore in the new UFL. What is the updated way of calling

domain = ufl.Mesh(ufl.VectorElement("Lagrange", cell_type.name, 2))

?

Thanks!

5 
import basix.ufl
domain = ufl.Mesh(basix.ufl.element("Lagrange", cell_type, degree, basix.LagrangeVariant.equispaced, shape=(gdim, ), gdim=gdim))

where cell_type is a string describing what cell you are using, degree the order of the geometry, and gdim the geometrical dimension of the mesh

1 Like
6

You are the best. Thanks for the quick answer! :slight_smile: