Error obtaining dofs between vertex and cell in 2D mesh

dolfinx
1

Dear FENICS community. I am trying to import meshes from gmsh using meshio to then study a problem with discontinuous elements in Dolfinx.

Provided that I have a 2D mesh file called mesh.xdmf, I execute

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh(name="Grid")
    ct = xdmf.read_meshtags(mesh, name="Grid")
mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "mt.xdmf", "r") as xdmf:
    ft = xdmf.read_meshtags(mesh, name="Grid")

Following @dokken example in Define a boundary from DG elements, I ran the lines

def boundary(x):
    return np.full(x.shape[1], True)

vertex_to_cell = mesh.topology.connectivity(0, mesh.topology.dim)
boundary_vertices = dolfinx.mesh.locate_entities(mesh, 0, boundary)
boundary_cells = []

for vertex in boundary_vertices:
    boundary_cells.append(vertex_to_cell.links(vertex))
boundary_cells = np.hstack(boundary_cells)

But this error comes out

Traceback (most recent call last):....
    boundary_cells.append(vertex_to_cell.links(vertex))
AttributeError: 'NoneType' object has no attribute 'links'

Since mesh.topology.connectivity(0, mesh.topology.dim) returns None. However, if I let the Python console open after the error and rerun line, sometimes it returns the vertex_to_cell dofs. I don’t now how this is happening.

On the other hand, the line

facet_to_cell = mesh.topology.connectivity(1, mesh.topology.dim)

returns the facet-to-cell dofs properly. Moreover, mesh.topology.connectivity(0, 1_or_2_here) returns None. The only one that seems to work is mesh.topology.connectivity(0, 0), i.e., vertex-to-vertex dofs

2

For the sake of reproducibility, I’ve made a similar code using a built-in mesh that one can write/read in the single code.

Please note that whenever you ask for a connectivity(a, b), you should make sure it as been created. This can be done by adding

mesh.topology.create_connectivity(a, b)

before requesting it as

a_to_b = mesh.topology.connectivity(a, b)

The reason for such an approach is that for many codes, you do not need these connectivities, which are time and memory consuming.

import dolfinx
import dolfinx.io
from mpi4py import MPI
import numpy as np


def create_and_write_mesh(filename: str):
    mesh = dolfinx.UnitSquareMesh(MPI.COMM_WORLD, 10, 10)
    mesh.topology.create_entities(mesh.topology.dim - 1)
    # Create a facet tag that has a unique marker per facet
    num_local_facets = mesh.topology.index_map(mesh.topology.dim-1).size_local
    facet_indices = np.arange(num_local_facets, dtype=np.int32)
    ft = dolfinx.MeshTags(mesh, mesh.topology.dim-1,
                          facet_indices, facet_indices)
    ft.name = "Facetmarkers"
    # Create a cell tag that has a unique marker per cell
    num_local_cells = mesh.topology.index_map(mesh.topology.dim).size_local
    cell_indices = np.arange(num_local_cells, dtype=np.int32)
    ct = dolfinx.MeshTags(mesh, mesh.topology.dim,
                          cell_indices, cell_indices)
    ct.name = "Cellmarkers"

    with dolfinx.io.XDMFFile(mesh.mpi_comm(), filename, "w") as xdmf:
        xdmf.write_mesh(mesh)
        xdmf.write_meshtags(ct)
        mesh.topology.create_connectivity(
            mesh.topology.dim - 1, mesh.topology.dim)
        xdmf.write_meshtags(ft)


filename = "cube.xdmf"
create_and_write_mesh("cube.xdmf")

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, filename, "r") as xdmf:
    mesh = xdmf.read_mesh()
    ct = xdmf.read_meshtags(mesh, name="Cellmarkers")
    mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
    ft = xdmf.read_meshtags(mesh, name="Facetmarkers")


def boundary(x):
    return np.full(x.shape[1], True)


mesh.topology.create_connectivity(0, mesh.topology.dim)
vertex_to_cell = mesh.topology.connectivity(0, mesh.topology.dim)
boundary_vertices = dolfinx.mesh.locate_entities_boundary(mesh, 0, boundary)
boundary_cells = []

for vertex in boundary_vertices:
    boundary_cells.append(vertex_to_cell.links(vertex))
boundary_cells = np.hstack(boundary_cells)
print(boundary_cells)
3 Likes