Support vertex-only 'meshes'

[jorgensd]

Describe new/missing feature

For mixed dimensional problems, it would be sensible to have 0D (vertex meshes). However, in the current implementation, there are many assumptions (such as the build_local_dual_graph, not clear what we should do in this case) that does not work with vertex meshes.

Changes that is required (and clear): dolfinx/cpp/dolfinx/mesh/graphbuild.cpp

constexpr mesh::CellType get_cell_type(int num_vertices, int tdim)
{
  switch (tdim)
  {
  case 0:
    return mesh::CellType::point;

dolfinx/python/dolfinx/mesh.py

_uflcell_to_dolfinxcell = {
    "vertex": CellType.point,

Suggestion user interface

from mpi4py import MPI
import ufl
import dolfinx
import numpy as np
points = np.array([[0, 0], [0, 1], [1, 0], [1, 1]], dtype=np.float64)
vertex = np.array([[0], [1], [2], [3]], dtype=np.int64)
ufl_tri = ufl.Mesh(ufl.VectorElement("DG", ufl.vertex, 0, dim=2))

def serial_partitioner(comm, n, m, topo):
    dest = np.zeros(topo.num_nodes, dtype=np.int32)

    return dolfinx.cpp.graph.AdjacencyList_int32(dest)

dolfinx.log.set_log_level(dolfinx.log.LogLevel.INFO)
tri_mesh = dolfinx.mesh.create_mesh(
    MPI.COMM_WORLD, vertex, points, ufl_tri, partitioner=serial_partitioner)

[garth-wells]

I think simpler would be to allow partitioner=None.

Another approach would be to create Topology and Geometry, and call the Mesh constructor directly. Some interface simplifications might be possible. Do we support a 'geometry' element for points?

[jorgensd]

I've got an implementation of this only requires a single line change in basix:

diff --git a/python/basix/ufl.py b/python/basix/ufl.py
index 18c98c7..0d7c421 100644
--- a/python/basix/ufl.py
+++ b/python/basix/ufl.py
@@ -103,6 +103,8 @@ class _ElementBase(_AbstractFiniteElement):
     ):
         """Initialise the element."""
         self._repr = repr
+        if cellname == "point":
+            cellname = "vertex"
         self._cellname = cellname
         self._reference_value_shape = reference_value_shape
         self._degree = degree

as the ufl CellName for a point is vertex.

MWE:

# Create a mesh consisting of points only
# Author: Jørgen S. Dokken
# SPDX-License-Identifier: MIT

from mpi4py import MPI
import basix.ufl

import dolfinx
import numpy as np
import numpy.typing as npt
import ufl

def create_point_mesh(comm: MPI.Intracomm, points:npt.NDArray[np.float32]|npt.NDArray[np.float64])->dolfinx.mesh.Mesh:
    """
    Create a mesh consisting of points only.

    Note:
        No nodes are shared between processes.

    Args:
        comm: MPI communicator to create the mesh on.
        points: Points local to the process in the mesh.
    """
    # Create mesh topology
    cells = np.arange(points.shape[0], dtype=np.int32).reshape(-1, 1)
    topology = dolfinx.cpp.mesh.Topology(MPI.COMM_WORLD, dolfinx.mesh.CellType.point)
    num_nodes_local = cells.shape[0]
    imap = dolfinx.common.IndexMap(MPI.COMM_WORLD, num_nodes_local)
    local_range = imap.local_range[0]
    igi = np.arange(num_nodes_local, dtype=np.int64)+local_range
    topology.set_index_map(0, imap)
    topology.set_connectivity(dolfinx.graph.adjacencylist(cells), 0, 0)

    # Create mesh geometry
    e = basix.ufl.element("Lagrange", "point", 0, shape=(points.shape[1],))
    c_el = dolfinx.fem.coordinate_element(e.basix_element)
    geometry = dolfinx.mesh.create_geometry(imap, cells, c_el._cpp_object, nodes,  igi)

    # Create DOLFINx mesh
    if points.dtype == np.float64:
        cpp_mesh = dolfinx.cpp.mesh.Mesh_float64(comm, topology, geometry._cpp_object)
    elif points.dtype == np.float32:
        cpp_mesh = dolfinx.cpp.mesh.Mesh_float32(comm, topology, geometry._cpp_object)
    else:
        raise RuntimeError(f"Unsupported dtype for mesh {points.dtype}")
    # Wrap as Python object
    return dolfinx.mesh.Mesh(cpp_mesh, domain = ufl.Mesh(e))

if __name__ == "__main__":
    np.random.seed(MPI.COMM_WORLD.rank)
    nodes = np.random.rand(100, 3)

    mesh = create_point_mesh(MPI.COMM_WORLD, nodes)
    print(mesh.topology.index_map(0).size_local, mesh.topology.index_map(0).size_global)
    print(mesh.geometry.x)
    V = dolfinx.fem.functionspace(mesh, ("DG", 0, (mesh.geometry.dim,)))
    u = dolfinx.fem.Function(V)
    u.interpolate(lambda x: (x[0],x[1],x[2]))

    with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "point_mesh.xdmf", "w") as xdmf:
        xdmf.write_mesh(mesh)
        xdmf.write_function(u)

Making it work for sub-meshes would require some extra work though (as we would need to make this construction in the create_submesh part of the code, and modify the sub-index map to have the appropriate ghosts (shouldn't be too hard though, as the sub-index map is all we need).

[jorgensd]

We are fairly close to this now, as the only thing needed is:

# Create a mesh consisting of points only
# Author: Jørgen S. Dokken
# SPDX-License-Identifier: MIT

from mpi4py import MPI
import basix.ufl

import dolfinx
import numpy as np
import numpy.typing as npt
import ufl

def create_point_mesh(comm: MPI.Intracomm, points:npt.NDArray[np.float32]|npt.NDArray[np.float64])->dolfinx.mesh.Mesh:
    """
    Create a mesh consisting of points only.

    Note:
        No nodes are shared between processes.

    Args:
        comm: MPI communicator to create the mesh on.
        points: Points local to the process in the mesh.
    """
    # Create mesh topology
    cells = np.arange(points.shape[0], dtype=np.int32).reshape(-1, 1)
    topology = dolfinx.cpp.mesh.Topology(MPI.COMM_WORLD, dolfinx.mesh.CellType.point)
    num_nodes_local = cells.shape[0]
    imap = dolfinx.common.IndexMap(MPI.COMM_WORLD, num_nodes_local)
    local_range = imap.local_range[0]
    igi = np.arange(num_nodes_local, dtype=np.int64)+local_range
    topology.set_index_map(0, imap)
    topology.set_connectivity(dolfinx.graph.adjacencylist(cells), 0, 0)

    # Create mesh geometry
    e = basix.ufl.element("Lagrange", "point", 0, shape=(points.shape[1],))
    c_el = dolfinx.fem.coordinate_element(e.basix_element)
    geometry = dolfinx.mesh.create_geometry(imap, cells, c_el._cpp_object, points,  igi)

    # Create DOLFINx mesh
    if points.dtype == np.float64:
        cpp_mesh = dolfinx.cpp.mesh.Mesh_float64(comm, topology, geometry._cpp_object)
    elif points.dtype == np.float32:
        cpp_mesh = dolfinx.cpp.mesh.Mesh_float32(comm, topology, geometry._cpp_object)
    else:
        raise RuntimeError(f"Unsupported dtype for mesh {points.dtype}")
    # Wrap as Python object
    return dolfinx.mesh.Mesh(cpp_mesh, domain = ufl.Mesh(e))

Could we add this constructor to the mesh generation routines?