MeshEditor equivalent

dolfinx
1

Hi all

I would like to know if there’s a way to create a 1D mesh directly from a list of vertices in dolfinx?

In dolfin I used to do:

        import fenics as f
        vertices = [1, 2, 3, 4]
        nb_points = len(vertices)
        nb_cells = nb_points - 1
        editor = f.MeshEditor()
        mesh = f.Mesh()
        editor.open(mesh, "interval", 1, 1)  # top. and geom. dimension are both 1
        editor.init_vertices(nb_points)  # number of vertices
        editor.init_cells(nb_cells)     # number of cells
        for i in range(0, nb_points):
            editor.add_vertex(i, np.array([vertices[i]]))
        for j in range(0, nb_cells):
            editor.add_cell(j, np.array([j, j+1]))
        editor.close()
2

In dolfinx, you can simply send in the mesh geometry and topology, as for instance shown in: dolfinx/test_cell.py at 8c61d6f632eb20470631522381c070e812221662 · FEniCS/dolfinx · GitHub

3 Likes
3

Thanks @dokken

Unfortunately, the following doesn’t work

from dolfinx.fem import FunctionSpace
from dolfinx.mesh import create_mesh
from mpi4py import MPI
import ufl


gdim, shape, degree = 2, "triangle", 1
cell = ufl.Cell(shape, geometric_dimension=gdim)
domain = ufl.Mesh(ufl.VectorElement("Lagrange", cell, degree))
x = [[0., 0., 0.], [0., 1., 0.], [1., 1., 0.]]
cells = [[0, 1, 2]]
mesh = create_mesh(MPI.COMM_WORLD, cells, x, domain)

V = FunctionSpace(mesh, ("P", 1))

Produces:

Traceback (most recent call last):
  File "/home/test_generate_mesh.py", line 14, in <module>
    V = FunctionSpace(mesh, ("P", 1))
  File "/usr/local/dolfinx-real/lib/python3.8/dist-packages/dolfinx/fem/function.py", line 440, in __init__
    super().__init__(mesh.ufl_domain(), ufl_element)
  File "/usr/local/lib/python3.9/dist-packages/ufl/functionspace.py", line 45, in __init__
    error("Non-matching cell of finite element and domain.")
  File "/usr/local/lib/python3.9/dist-packages/ufl/log.py", line 158, in error
    raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Non-matching cell of finite element and domain.
4

This is because you send in three-dimensional coordinates, while the VectorElement defaults to a 2D mesh (not a manifold in 3D).
With some slight alterations to your code:

from dolfinx.fem import FunctionSpace
from dolfinx.mesh import create_mesh
from mpi4py import MPI
import ufl
import numpy as np

gdim, shape, degree = 2, "triangle", 1
cell = ufl.Cell(shape, geometric_dimension=gdim)
domain = ufl.Mesh(ufl.VectorElement("Lagrange", cell, degree))
x = np.array([[0., 0., 0.], [0., 1., 0.], [1., 1., 0.]])
cells = np.array([[0, 1, 2]], dtype=np.int64)
mesh = create_mesh(MPI.COMM_WORLD, cells, x[:, :gdim], domain)

V = FunctionSpace(mesh, ("P", 1))
print(V.tabulate_dof_coordinates(), mesh.geometry.dim)

the code runs without error and prints:

[[0. 0. 0.]
 [0. 1. 0.]
 [1. 1. 0.]] 2
5

Ok. I’m mainly interested in creating 1D meshes from a list of vertices. Is there a simple way to do that from this snipset?

6

Simply change the cell type and geometrical dimension, as follows:

from dolfinx.fem import FunctionSpace
from dolfinx.mesh import create_mesh
from mpi4py import MPI
import ufl
import numpy as np

gdim, shape, degree = 1, "interval", 1
cell = ufl.Cell(shape, geometric_dimension=gdim)
domain = ufl.Mesh(ufl.VectorElement("Lagrange", cell, degree))
x = np.array([[0], [0.2], [0.6], [1], [1.3]])
cells = np.array([[0, 1], [1, 2], [2, 3], [3, 4]], dtype=np.int64)
mesh = create_mesh(MPI.COMM_WORLD, cells, x, domain)

V = FunctionSpace(mesh, ("P", 1))
print(V.tabulate_dof_coordinates(), mesh.geometry.dim)
print(mesh.topology.index_map(mesh.topology.dim).size_local)
7

Hi @dokken,

Using domain = ufl.Mesh(ufl.VectorElement("Lagrange", cell, degree)) produces a lot of UFL warnings along the way (when restricting the dx measurement eg. dx(1)):

DeprecationWarning: Converting elements created in UFL to Basix elements is deprecated. You should create the elements directly using basix.ufl.element instead

Do you know what would be the equivalent for this in basix?

Thanks,
Rem

EDIT: nvm managed to find the basix equivalent, it was more straightforward than I thought!

degree = 1
domain = ufl.Mesh(
            basix.ufl.element(
                basix.ElementFamily.P,
                "interval",
                degree,
                shape=(1,)
            )
        )
8

You can even do

degree = 1
domain = ufl.Mesh(
            basix.ufl.element(
                "Lagrange",
                "interval",
                degree,
                shape=(1,)
            )
        )
1 Like