diff --git a/FIAT/discontinuous_pc.py b/FIAT/discontinuous_pc.py
index dbb56690..8973a58f 100644
--- a/FIAT/discontinuous_pc.py
+++ b/FIAT/discontinuous_pc.py
@@ -7,29 +7,17 @@
 # Modified by David A. Ham (david.ham@imperial.ac.uk), 2018
 
 from FIAT import finite_element, polynomial_set, dual_set, functional
-from FIAT.reference_element import (Point,
-                                    DefaultLine,
-                                    UFCInterval,
-                                    UFCQuadrilateral,
-                                    UFCHexahedron,
-                                    UFCTriangle,
-                                    UFCTetrahedron,
-                                    make_affine_mapping,
-                                    flatten_reference_cube)
+from FIAT.reference_element import (make_affine_mapping,
+                                    flatten_reference_cube,
+                                    cell_to_simplex)
 from FIAT.P0 import P0Dual
 import numpy as np
 
-hypercube_simplex_map = {Point(): Point(),
-                         DefaultLine(): DefaultLine(),
-                         UFCInterval(): UFCInterval(),
-                         UFCQuadrilateral(): UFCTriangle(),
-                         UFCHexahedron(): UFCTetrahedron()}
-
 
 class DPC0(finite_element.CiarletElement):
     def __init__(self, ref_el):
         flat_el = flatten_reference_cube(ref_el)
-        poly_set = polynomial_set.ONPolynomialSet(hypercube_simplex_map[flat_el], 0)
+        poly_set = polynomial_set.ONPolynomialSet(cell_to_simplex(flat_el), 0)
         dual = P0Dual(ref_el)
         # Implement entity_permutations when we handle that for HigherOrderDPC.
         # Currently, orientation_tuples in P0Dual(ref_el).entity_permutations
@@ -58,7 +46,7 @@ def __init__(self, ref_el, flat_el, degree):
 
         # Change coordinates here.
         # Vertices of the simplex corresponding to the reference element.
-        v_simplex = hypercube_simplex_map[flat_el].get_vertices()
+        v_simplex = cell_to_simplex(flat_el).get_vertices()
         # Vertices of the reference element.
         v_hypercube = flat_el.get_vertices()
         # For the mapping, first two vertices are unchanged in all dimensions.
@@ -74,12 +62,12 @@ def __init__(self, ref_el, flat_el, degree):
 
         # make nodes by getting points
         # need to do this dimension-by-dimension, facet-by-facet
-        top = hypercube_simplex_map[flat_el].get_topology()
+        top = cell_to_simplex(flat_el).get_topology()
 
         cur = 0
         for dim in sorted(top):
             for entity in sorted(top[dim]):
-                pts_cur = hypercube_simplex_map[flat_el].make_points(dim, entity, degree)
+                pts_cur = cell_to_simplex(flat_el).make_points(dim, entity, degree)
                 pts_cur = [tuple(np.matmul(A, np.array(x)) + b) for x in pts_cur]
                 nodes_cur = [functional.PointEvaluation(flat_el, x)
                              for x in pts_cur]
@@ -102,7 +90,7 @@ class HigherOrderDPC(finite_element.CiarletElement):
 
     def __init__(self, ref_el, degree):
         flat_el = flatten_reference_cube(ref_el)
-        poly_set = polynomial_set.ONPolynomialSet(hypercube_simplex_map[flat_el], degree)
+        poly_set = polynomial_set.ONPolynomialSet(cell_to_simplex(flat_el), degree)
         dual = DPCDualSet(ref_el, flat_el, degree)
         formdegree = flat_el.get_spatial_dimension()  # n-form
         super().__init__(poly_set=poly_set,
diff --git a/FIAT/orientation_utils.py b/FIAT/orientation_utils.py
index 2852f33f..92cc68db 100644
--- a/FIAT/orientation_utils.py
+++ b/FIAT/orientation_utils.py
@@ -97,7 +97,7 @@ def _make_axis_perms_tensorproduct(cells, dim):
     This is the single sources of extrinsic orientations and corresponding axis permutations.
 
     """
-    from FIAT.reference_element import UFCInterval
+    from FIAT.reference_element import LINE
 
     # Handle extrinsic orientations.
     # This is complex and we need to think to make this function more general.
@@ -116,8 +116,8 @@ def _make_axis_perms_tensorproduct(cells, dim):
         #          dim == (2, 1) ->
         #          triangle x interval (1 possible extrinsic orientation).
         axis_perms = (tuple(range(nprod)), )  # Identity: no permutations
-    elif len(set(cells)) == 1 and isinstance(cells[0], UFCInterval):
-        # Tensor product of intervals.
+    elif len(set(cells)) == 1 and cells[0].get_shape() == LINE:
+        # Tensor product of intervals (any 1D reference cell implementation)
         # Example: interval x interval x interval x interval
         #          dim == (0, 1, 1, 1) ->
         #          point x interval x interval x interval  (1! * 3! possible extrinsic orientations).
diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py
index 3427e93c..aae357bd 100644
--- a/FIAT/reference_element.py
+++ b/FIAT/reference_element.py
@@ -133,7 +133,7 @@ class Cell:
     """Abstract class for a reference cell.  Provides accessors for
     geometry (vertex coordinates) as well as topology (orderings of
     vertices that make up edges, faces, etc."""
-    def __init__(self, shape, vertices, topology):
+    def __init__(self, shape, vertices, topology, sub_entities=None):
         """The constructor takes a shape code, the physical vertices expressed
         as a list of tuples of numbers, and the topology of a cell.
 
@@ -145,24 +145,25 @@ def __init__(self, shape, vertices, topology):
         self.vertices = vertices
         self.topology = topology
 
-        # Given the topology, work out for each entity in the cell,
-        # which other entities it contains.
-        self.sub_entities = {}
-        for dim, entities in topology.items():
-            self.sub_entities[dim] = {}
-
-            for e, v in entities.items():
-                vertices = frozenset(v)
-                sub_entities = []
-
-                for dim_, entities_ in topology.items():
-                    for e_, vertices_ in entities_.items():
-                        if vertices.issuperset(vertices_):
-                            sub_entities.append((dim_, e_))
-
-                # Sort for the sake of determinism and by UFC conventions
-                self.sub_entities[dim][e] = sorted(sub_entities)
-
+        if sub_entities is not None:
+            self.sub_entities = sub_entities
+        else:
+            # If sub entity list not provided
+            # Given the topology, work out for each entity in the cell,
+            # which other entities it contains.
+            self.sub_entities = {}
+            for dim, entities in topology.items():
+                self.sub_entities[dim] = {}
+
+                for e, v in entities.items():
+                    vertices = frozenset(v)
+                    sub_entities = []
+                    for dim_, entities_ in topology.items():
+                        for e_, vertices_ in entities_.items():
+                            if vertices.issuperset(vertices_):
+                                sub_entities.append((dim_, e_))
+
+                    self.sub_entities[dim][e] = sorted(list(sub_entities))
         # Build super-entity dictionary by inverting the sub-entity dictionary
         self.super_entities = {dim: {entity: [] for entity in topology[dim]} for dim in topology}
         for dim0 in topology:
@@ -183,7 +184,6 @@ def __init__(self, shape, vertices, topology):
                     neighbors = children if dim1 < dim0 else parents
                     d01_entities = tuple(e for d, e in neighbors if d == dim1)
                     self.connectivity[(dim0, dim1)].append(d01_entities)
-
         # Dictionary with derived cells
         self._split_cache = {}
 
@@ -387,14 +387,14 @@ class SimplicialComplex(Cell):
 
     This consists of list of vertex locations and a topology map defining facets.
     """
-    def __init__(self, shape, vertices, topology):
+    def __init__(self, shape, vertices, topology, sub_entities=None):
         # Make sure that every facet has the right number of vertices to be
         # a simplex.
         for dim in topology:
             for entity in topology[dim]:
                 assert len(topology[dim][entity]) == dim + 1
 
-        super().__init__(shape, vertices, topology)
+        super().__init__(shape, vertices, topology, sub_entities)
 
     def compute_normal(self, facet_i, cell=None):
         """Returns the unit normal vector to facet i of codimension 1."""
@@ -528,7 +528,7 @@ def make_points(self, dim, entity_id, order, variant=None, interior=1):
         facet of dimension dim.  Order indicates how many points to
         include in each direction."""
         if dim == 0:
-            return (self.get_vertices()[entity_id], )
+            return (self.get_vertices()[self.get_topology()[dim][entity_id][0]],)
         elif 0 < dim <= self.get_spatial_dimension():
             entity_verts = \
                 self.get_vertices_of_subcomplex(
@@ -1155,7 +1155,7 @@ def compute_normal(self, i):
 class TensorProductCell(Cell):
     """A cell that is the product of FIAT cells."""
 
-    def __init__(self, *cells):
+    def __init__(self, *cells, sub_entities=None):
         # Vertices
         vertices = tuple(tuple(chain(*coords))
                          for coords in product(*[cell.get_vertices()
@@ -1178,7 +1178,7 @@ def __init__(self, *cells):
             topology[dim] = dict(enumerate(topology[dim][key]
                                            for key in sorted(topology[dim])))
 
-        super().__init__(TENSORPRODUCT, vertices, topology)
+        super().__init__(TENSORPRODUCT, vertices, topology, sub_entities)
         self.cells = tuple(cells)
 
     def __repr__(self):
@@ -1421,7 +1421,7 @@ def is_macrocell(self):
 class Hypercube(Cell):
     """Abstract class for a reference hypercube"""
 
-    def __init__(self, dimension, product):
+    def __init__(self, dimension, product, sub_entities=None):
         self.dimension = dimension
         self.shape = hypercube_shapes[dimension]
 
@@ -1429,7 +1429,7 @@ def __init__(self, dimension, product):
         verts = product.get_vertices()
         topology = flatten_entities(pt)
 
-        super().__init__(self.shape, verts, topology)
+        super().__init__(self.shape, verts, topology, sub_entities)
 
         self.product = product
         self.unflattening_map = compute_unflattening_map(pt)
@@ -1872,3 +1872,10 @@ def max_complex(complexes):
         return max_cell
     else:
         raise ValueError("Cannot find the maximal complex")
+
+
+def cell_to_simplex(cell):
+    if cell.is_simplex():
+        return cell
+    else:
+        return ufc_simplex(cell.get_dimension())
diff --git a/finat/__init__.py b/finat/__init__.py
index bde435a1..34ea55b8 100644
--- a/finat/__init__.py
+++ b/finat/__init__.py
@@ -8,7 +8,8 @@
                             Lagrange, Real, Serendipity,                                     # noqa: F401
                             TrimmedSerendipityCurl, TrimmedSerendipityDiv,                   # noqa: F401
                             TrimmedSerendipityEdge, TrimmedSerendipityFace,                  # noqa: F401
-                            Nedelec, NedelecSecondKind, RaviartThomas, Regge)                # noqa: F401
+                            Nedelec, NedelecSecondKind, RaviartThomas, Regge,                # noqa: F401
+                            FuseElement)                                                     # noqa: F401
 from .spectral import (GaussLobattoLegendre, GaussLegendre, KongMulderVeldhuizen,            # noqa: F401
                        Legendre, IntegratedLegendre,                                         # noqa: F401
                        FDMLagrange, FDMQuadrature, FDMDiscontinuousLagrange,                 # noqa: F401
diff --git a/finat/element_factory.py b/finat/element_factory.py
index 8563aa60..b5d5715f 100644
--- a/finat/element_factory.py
+++ b/finat/element_factory.py
@@ -109,15 +109,24 @@
 element is supported, but must be handled specially because it doesn't
 have a direct FInAT equivalent."""
 
+hexahedron_tpc = ufl.TensorProductCell(ufl.interval, ufl.interval, ufl.interval)
+quadrilateral_tpc = ufl.TensorProductCell(ufl.interval, ufl.interval)
+
 
 @cache
 def as_fiat_cell(cell):
     """Convert a ufl cell to a FIAT cell.
 
     :arg cell: the :class:`ufl.Cell` to convert."""
+    if isinstance(cell, str):
+        cell = finat.ufl.as_cell(cell)
     if not isinstance(cell, ufl.AbstractCell):
         raise ValueError("Expecting a UFL Cell")
-    return ufc_cell(cell)
+
+    if hasattr(cell, "to_fiat"):
+        return cell.to_fiat()
+    else:
+        return ufc_cell(cell)
 
 
 @singledispatch
@@ -325,23 +334,28 @@ def convert_restrictedelement(element, **kwargs):
     return finat.RestrictedElement(finat_elem, element.restriction_domain()), deps
 
 
-hexahedron_tpc = ufl.TensorProductCell(ufl.interval, ufl.interval, ufl.interval)
-quadrilateral_tpc = ufl.TensorProductCell(ufl.interval, ufl.interval)
+@convert.register(finat.ufl.FuseElement)
+def convert_fuse_element(element, **kwargs):
+    return finat.fiat_elements.FuseElement(element.triple), set()
+
+
 _cache = weakref.WeakKeyDictionary()
 
 
-def create_element(ufl_element, shape_innermost=True, shift_axes=0, restriction=None):
+def create_element(ufl_element, shape_innermost=True, shift_axes=0, restriction=None, cell_backend=finat.ufl.CellBackend.FIAT):
     """Create a FInAT element (suitable for tabulating with) given a UFL element.
 
     :arg ufl_element: The UFL element to create a FInAT element from.
     :arg shape_innermost: Vector/tensor indices come after basis function indices
     :arg restriction: cell restriction in interior facet integrals
                       (only for runtime tabulated elements)
+    :arg cell_backend: Enum determining the default cell type to use.
     """
     finat_element, deps = _create_element(ufl_element,
                                           shape_innermost=shape_innermost,
                                           shift_axes=shift_axes,
-                                          restriction=restriction)
+                                          restriction=restriction,
+                                          cell_backend=cell_backend)
     return finat_element
 
 
diff --git a/finat/fiat_elements.py b/finat/fiat_elements.py
index 22408cf0..f52c55ca 100644
--- a/finat/fiat_elements.py
+++ b/finat/fiat_elements.py
@@ -439,3 +439,8 @@ def __init__(self, cell, degree, **kwargs):
 class NedelecSecondKind(VectorFiatElement):
     def __init__(self, cell, degree, **kwargs):
         super().__init__(FIAT.NedelecSecondKind(cell, degree, **kwargs))
+
+
+class FuseElement(FiatElement):
+    def __init__(self, triple):
+        super(FuseElement, self).__init__(triple.to_fiat())
diff --git a/finat/ufl/__init__.py b/finat/ufl/__init__.py
index 21a7d13d..9da7f2db 100644
--- a/finat/ufl/__init__.py
+++ b/finat/ufl/__init__.py
@@ -14,8 +14,9 @@
 from finat.ufl.brokenelement import BrokenElement  # noqa: F401
 from finat.ufl.enrichedelement import EnrichedElement, NodalEnrichedElement  # noqa: F401
 from finat.ufl.finiteelement import FiniteElement  # noqa: F401
-from finat.ufl.finiteelementbase import FiniteElementBase  # noqa: F401
+from finat.ufl.finiteelementbase import FiniteElementBase, as_cell, CellBackend  # noqa: F401
 from finat.ufl.hdivcurl import HCurlElement, HDivElement, WithMapping, HDiv, HCurl  # noqa: F401
 from finat.ufl.mixedelement import MixedElement, TensorElement, VectorElement  # noqa: F401
 from finat.ufl.restrictedelement import RestrictedElement  # noqa: F401
 from finat.ufl.tensorproductelement import TensorProductElement  # noqa: F401
+from finat.ufl.fuseelement import FuseElement  # noqa: F401
diff --git a/finat/ufl/elementlist.py b/finat/ufl/elementlist.py
index 7bd949a3..865ad00f 100644
--- a/finat/ufl/elementlist.py
+++ b/finat/ufl/elementlist.py
@@ -18,6 +18,7 @@
 # Modified by Pablo Brubeck, 2024
 
 import warnings
+import inspect
 
 from numpy import asarray
 
@@ -416,6 +417,10 @@ def canonical_element_description(family, cell, order, form_degree):
             raise ValueError("Need dimension to handle element aliases.")
         (family, order) = aliases[family](family, tdim, order, form_degree)
 
+    # Check that we don't have a raw FUSE object
+    if not isinstance(family, str) and 'fuse' in inspect.getmodule(family).__name__:
+        raise ValueError("Received unconverted FUSE triple - did you forget to call to_ufl()?")
+
     # Check that the element family exists
     if family not in ufl_elements:
         raise ValueError(f"Unknown finite element '{family}'.")
diff --git a/finat/ufl/finiteelement.py b/finat/ufl/finiteelement.py
index dade48b0..bc2a0bdb 100644
--- a/finat/ufl/finiteelement.py
+++ b/finat/ufl/finiteelement.py
@@ -11,9 +11,9 @@
 # Modified by Massimiliano Leoni, 2016
 # Modified by Matthew Scroggs, 2023
 
-from ufl.cell import TensorProductCell, as_cell
+from ufl.cell import TensorProductCell
 from finat.ufl.elementlist import canonical_element_description, simplices
-from finat.ufl.finiteelementbase import FiniteElementBase
+from finat.ufl.finiteelementbase import FiniteElementBase, as_cell
 from ufl.utils.formatting import istr
 
 
diff --git a/finat/ufl/finiteelementbase.py b/finat/ufl/finiteelementbase.py
index 48d0c24a..1dcac489 100644
--- a/finat/ufl/finiteelementbase.py
+++ b/finat/ufl/finiteelementbase.py
@@ -14,9 +14,10 @@
 from abc import abstractmethod, abstractproperty
 from hashlib import md5
 from typing import Sequence, Union
+from enum import Enum
 
 from ufl import pullback
-from ufl.cell import AbstractCell, as_cell
+from ufl.cell import AbstractCell, as_cell as as_cell_ufl
 from ufl.finiteelement import AbstractFiniteElement
 from ufl.utils.sequences import product
 
@@ -291,3 +292,21 @@ def pullback(self):
             return supported_pullbacks[self.mapping()]
         except KeyError:
             raise ValueError(f"Unsupported mapping: {self.mapping()}")
+
+
+class CellBackend(Enum):
+    FIAT = 1
+    FUSE = 2
+
+
+def as_cell(cell: AbstractCell | str | tuple[AbstractCell, ...], cell_backend: CellBackend = CellBackend.FIAT) -> AbstractCell:
+    if isinstance(cell, str) and cell_backend == CellBackend.FUSE:
+        try:
+            import fuse
+        except ModuleNotFoundError as exc:
+            raise ModuleNotFoundError(
+                "FUSE cell construction requires the optional 'fuse' dependency. "
+            ) from exc
+        return fuse.constructCellComplex(cell)
+    else:
+        return as_cell_ufl(cell)
diff --git a/finat/ufl/fuseelement.py b/finat/ufl/fuseelement.py
new file mode 100644
index 00000000..d46837cf
--- /dev/null
+++ b/finat/ufl/fuseelement.py
@@ -0,0 +1,49 @@
+"""Element."""
+# -*- coding: utf-8 -*-
+# Copyright (C) 2025 India Marsden
+#
+# SPDX-License-Identifier:    LGPL-3.0-or-later
+
+from finat.ufl.finiteelementbase import FiniteElementBase
+
+
+class FuseElement(FiniteElementBase):
+    """
+    A finite element defined using FUSE.
+
+    :arg triple: An ElementTriple object defined with FUSE
+    :arg cell: Optional (defaults to triple.cell) The cell the element is defined on
+
+    """
+
+    def __init__(self, triple, cell=None):
+        try:
+            import fuse
+        except ModuleNotFoundError as exc:
+            raise ModuleNotFoundError(
+                "FUSE element creation requires the optional 'fuse' dependency. "
+            ) from exc
+        assert isinstance(triple, fuse.ElementTriple)
+        self.triple = triple
+        if not cell:
+            cell = self.triple.cell.to_ufl()
+
+        degree = self.triple.degree
+        super().__init__("FUSE", cell, degree, None, triple.get_value_shape())
+
+    def __repr__(self):
+        return repr(self.triple)
+
+    def __str__(self):
+        return f"<FuseElem on {self.triple.cell}>"
+
+    @property
+    def sobolev_space(self):
+        return self.triple.spaces[2].ufl_sobolev_space(self.triple.form_degree,
+                                                       self.triple.cell.dim())
+
+    def mapping(self):
+        return self.triple.spaces[2].mapping()
+
+    def reconstruct(self, family=None, cell=None, degree=None, quad_scheme=None, variant=None):
+        return FuseElement(self.triple, cell=cell)
diff --git a/finat/ufl/mixedelement.py b/finat/ufl/mixedelement.py
index 3e07c4fa..1c984b74 100644
--- a/finat/ufl/mixedelement.py
+++ b/finat/ufl/mixedelement.py
@@ -13,10 +13,10 @@
 
 import numpy as np
 
-from ufl.cell import CellSequence, as_cell
+from ufl.cell import CellSequence
 from ufl.domain import MeshSequence
 from finat.ufl.finiteelement import FiniteElement
-from finat.ufl.finiteelementbase import FiniteElementBase, shifted_sub_degrees
+from finat.ufl.finiteelementbase import FiniteElementBase, shifted_sub_degrees, as_cell
 from ufl.permutation import compute_indices
 from ufl.pullback import MixedPullback, SymmetricPullback
 from ufl.utils.indexflattening import flatten_multiindex, shape_to_strides, unflatten_index
diff --git a/finat/ufl/tensorproductelement.py b/finat/ufl/tensorproductelement.py
index 9a9a9828..29896f20 100644
--- a/finat/ufl/tensorproductelement.py
+++ b/finat/ufl/tensorproductelement.py
@@ -13,8 +13,9 @@
 
 from itertools import chain
 
-from ufl.cell import TensorProductCell, as_cell
-from finat.ufl.finiteelementbase import FiniteElementBase, shifted_sub_degrees
+from ufl.cell import TensorProductCell
+from finat.ufl.finiteelementbase import FiniteElementBase, shifted_sub_degrees, as_cell
+
 from ufl.sobolevspace import DirectionalSobolevSpace
 
 
diff --git a/pyproject.toml b/pyproject.toml
index e2668d68..96f46ad7 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -34,6 +34,9 @@ doc = [
   "sphinx",
 ]
 test = ["pytest"]
+fuse = [
+  "fuse-element @ git+https://github.com/firedrakeproject/fuse.git",
+]
 
 [tool.setuptools]
 packages = ["FIAT", "finat", "finat.ufl", "gem"]
diff --git a/test/FIAT/unit/test_reference_element.py b/test/FIAT/unit/test_reference_element.py
index d166ba6e..0ce63002 100644
--- a/test/FIAT/unit/test_reference_element.py
+++ b/test/FIAT/unit/test_reference_element.py
@@ -22,6 +22,7 @@
 from FIAT.reference_element import UFCInterval, UFCTriangle, UFCTetrahedron
 from FIAT.reference_element import Point, TensorProductCell, UFCQuadrilateral, UFCHexahedron
 from FIAT.reference_element import is_ufc, is_hypercube, default_simplex, flatten_reference_cube, Hypercube
+from FIAT.reference_element import Cell
 
 point = Point()
 interval = UFCInterval()
@@ -86,6 +87,51 @@ def test_ufc_connectivity_Dx(cell):
         assert connectivity[0] == tuple(range(len(connectivity[0])))
 
 
+def test_explicit_sub_entities_reproduces_default():
+    """An explicit sub_entities mapping equal to the one Cell would
+    compute automatically must produce identical super_entities and
+    connectivity."""
+    topology = {0: {0: (0,), 1: (1,)}, 1: {0: (0, 1)}}
+    vertices = ((0.0,), (1.0,))
+
+    reference = Cell("interval", vertices, topology)
+    explicit = Cell("interval", vertices, topology,
+                    sub_entities=reference.sub_entities)
+
+    assert explicit.sub_entities == reference.sub_entities
+    assert explicit.super_entities == reference.super_entities
+    assert explicit.connectivity == reference.connectivity
+
+
+def test_explicit_sub_entities_custom_mapping():
+    """A deliberately different (but self-consistent) explicit
+    sub_entities mapping should be used as-is, and super_entities /
+    connectivity should be derived from it rather than recomputed
+    from the topology."""
+    topology = {0: {0: (0,), 1: (1,)}, 1: {0: (0, 1)}}
+    vertices = ((0.0,), (1.0,))
+
+    # Deliberately omit the "each entity is a sub-entity of itself"
+    # entries that the automatic computation would include.
+    custom_sub_entities = {
+        0: {0: [], 1: []},
+        1: {0: [(0, 0), (0, 1)]},
+    }
+
+    cell = Cell("interval", vertices, topology,
+                sub_entities=custom_sub_entities)
+
+    assert cell.sub_entities == custom_sub_entities
+    assert cell.super_entities == {
+        0: {0: [(1, 0)], 1: [(1, 0)]},
+        1: {0: []},
+    }
+    assert cell.connectivity[(1, 0)] == [(0, 1)]
+    assert cell.connectivity[(0, 1)] == [(0,), (0,)]
+    assert cell.connectivity[(1, 1)] == [()]
+    assert cell.connectivity[(0, 0)] == [(), ()]
+
+
 @pytest.mark.parametrize(('cell', 'volume'),
                          [pytest.param(point, 1, marks=pytest.mark.xfail(conditional=sys.version_info < (3, 6))),
                           (interval, 1),