diff --git a/FIAT/expansions.py b/FIAT/expansions.py
index 3049d44f..5beebc23 100644
--- a/FIAT/expansions.py
+++ b/FIAT/expansions.py
@@ -92,7 +92,7 @@ def dubiner_recurrence(dim, n, order, ref_pts, Jinv, scale, variant=None):
     dX = pad_jacobian(Jinv, pad_dim)
 
     phi0 = numpy.array([sum((ref_pts[i] - ref_pts[i] for i in range(dim)), 0.0)])
-    results = [numpy.zeros((num_members,) + (dim,)*k + phi0.shape[1:], dtype=phi0.dtype)
+    results = [numpy.zeros((num_members,) + (len(dX[0]),)*k + phi0.shape[1:], dtype=phi0.dtype)
                for k in range(order+1)]
 
     phi, dphi, ddphi = results + [None] * (2-order)
@@ -263,22 +263,23 @@ def __new__(cls, *args, **kwargs):
         reference element."""
         if cls is not ExpansionSet:
             return super().__new__(cls)
+        ref_el = args[0]
+        shape = ref_el.get_shape()
         try:
-            ref_el = args[0]
             expansion_set = {
                 reference_element.POINT: PointExpansionSet,
                 reference_element.LINE: LineExpansionSet,
                 reference_element.TRIANGLE: TriangleExpansionSet,
                 reference_element.TETRAHEDRON: TetrahedronExpansionSet,
-            }[ref_el.get_shape()]
-            return expansion_set(*args, **kwargs)
+            }[shape]
         except KeyError:
-            raise ValueError("Invalid reference element type.")
+            raise ValueError(f"Invalid reference element type {type(ref_el).__name__}.")
+        return expansion_set(*args, **kwargs)
 
     def __init__(self, ref_el, scale=None, variant=None):
         self.ref_el = ref_el
         self.variant = variant
-        sd = ref_el.get_spatial_dimension()
+        sd = ref_el.get_topological_dimension()
         top = ref_el.get_topology()
         base_ref_el = reference_element.default_simplex(sd)
         base_verts = base_ref_el.get_vertices()
@@ -303,7 +304,7 @@ def reconstruct(self, ref_el=None, scale=None, variant=None):
 
     def get_scale(self, n, cell=0):
         scale = self.scale
-        sd = self.ref_el.get_spatial_dimension()
+        sd = self.ref_el.get_topological_dimension()
         if isinstance(scale, str):
             vol = self.ref_el.volume_of_subcomplex(sd, cell)
             scale = scale.lower()
@@ -335,11 +336,12 @@ def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None):
         ref_pts = numpy.add(numpy.dot(pts, A.T), b).T
         Jinv = A if direction is None else numpy.dot(A, direction)[:, None]
         sd = self.ref_el.get_spatial_dimension()
+        tdim = self.ref_el.get_topological_dimension()
         scale = self.get_scale(n, cell=cell)
-        phi = dubiner_recurrence(sd, n, lorder, ref_pts, Jinv,
+        phi = dubiner_recurrence(tdim, n, lorder, ref_pts, Jinv,
                                  scale, variant=self.variant)
         if self.continuity == "C0":
-            phi = C0_basis(sd, n, phi)
+            phi = C0_basis(tdim, n, phi)
 
         # Pack linearly independent components into a dictionary
         result = {(0,) * sd: numpy.asarray(phi[0])}
@@ -418,7 +420,8 @@ def tabulate_normal_jumps(self, n, ref_pts, facet, order=0):
         :returns: a numpy array of tabulations of normal derivative jumps.
         """
         sd = self.ref_el.get_spatial_dimension()
-        transform = self.ref_el.get_entity_transform(sd-1, facet)
+        tdim = self.ref_el.get_topological_dimension()
+        transform = self.ref_el.get_entity_transform(tdim-1, facet)
         pts = transform(ref_pts)
         cell_point_map = compute_cell_point_map(self.ref_el, pts, unique=False)
         cell_node_map = self.get_cell_node_map(n)
@@ -507,7 +510,7 @@ def get_dmats(self, degree, cell=0):
         if degree == 0:
             return cache.setdefault(key, numpy.zeros((self.ref_el.get_spatial_dimension(), 1, 1), "d"))
 
-        D = self.ref_el.get_dimension()
+        D = self.ref_el.get_topological_dimension()
         top = self.ref_el.get_topology()
         verts = self.ref_el.get_vertices_of_subcomplex(top[D][cell])
         pts = reference_element.make_lattice(verts, degree, variant="gl")
@@ -556,7 +559,7 @@ def __eq__(self, other):
 class PointExpansionSet(ExpansionSet):
     """Evaluates the point basis on a point reference element."""
     def __init__(self, ref_el, **kwargs):
-        if ref_el.get_spatial_dimension() != 0:
+        if ref_el.get_topological_dimension() != 0:
             raise ValueError("Must have a point")
         super().__init__(ref_el, **kwargs)
 
@@ -570,8 +573,8 @@ def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None):
 class LineExpansionSet(ExpansionSet):
     """Evaluates the Legendre basis on a line reference element."""
     def __init__(self, ref_el, **kwargs):
-        if ref_el.get_spatial_dimension() != 1:
-            raise Exception("Must have a line")
+        if ref_el.get_topological_dimension() != 1:
+            raise ValueError("Must have a line")
         super().__init__(ref_el, **kwargs)
 
     def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None):
@@ -600,16 +603,16 @@ class TriangleExpansionSet(ExpansionSet):
     """Evaluates the orthonormal Dubiner basis on a triangular
     reference element."""
     def __init__(self, ref_el, **kwargs):
-        if ref_el.get_spatial_dimension() != 2:
-            raise Exception("Must have a triangle")
+        if ref_el.get_topological_dimension() != 2:
+            raise ValueError("Must have a triangle")
         super().__init__(ref_el, **kwargs)
 
 
 class TetrahedronExpansionSet(ExpansionSet):
     """Collapsed orthonormal polynomial expansion on a tetrahedron."""
     def __init__(self, ref_el, **kwargs):
-        if ref_el.get_spatial_dimension() != 3:
-            raise Exception("Must be a tetrahedron")
+        if ref_el.get_topological_dimension() != 3:
+            raise ValueError("Must be a tetrahedron")
         super().__init__(ref_el, **kwargs)
 
 
@@ -627,7 +630,7 @@ def polynomial_dimension(ref_el, n, continuity=None):
     elif continuity == "C0":
         space_dimension = sum(math.comb(n - 1, dim) * len(top[dim]) for dim in top)
     else:
-        dim = ref_el.get_spatial_dimension()
+        dim = ref_el.get_topological_dimension()
         space_dimension = math.comb(n + dim, dim) * len(top[dim])
     return space_dimension
 
@@ -641,7 +644,7 @@ def polynomial_entity_ids(ref_el, n, continuity=None):
     :returns: a dict of dicts mapping dimension and entity id to basis functions.
     """
     top = ref_el.get_topology()
-    sd = ref_el.get_spatial_dimension()
+    sd = ref_el.get_topological_dimension()
     entity_ids = {}
     cur = 0
     for dim in sorted(top):
@@ -668,7 +671,7 @@ def polynomial_cell_node_map(ref_el, n, continuity=None):
     :returns: a numpy array mapping cell id to basis functions supported on that cell.
     """
     top = ref_el.get_topology()
-    sd = ref_el.get_spatial_dimension()
+    sd = ref_el.get_topological_dimension()
 
     entity_ids = polynomial_entity_ids(ref_el, n, continuity)
     ref_entity_ids = polynomial_entity_ids(ref_el.construct_subelement(sd), n, continuity)
@@ -697,7 +700,7 @@ def compute_cell_point_map(ref_el, pts, unique=True, tol=1E-12):
     :returns: a dict mapping cell id to the point ids nearest to that cell.
     """
     top = ref_el.get_topology()
-    sd = ref_el.get_spatial_dimension()
+    sd = ref_el.get_topological_dimension()
     if len(top[sd]) == 1:
         return {0: Ellipsis}
 
diff --git a/FIAT/finite_element.py b/FIAT/finite_element.py
index f8d4447b..2656c8e1 100644
--- a/FIAT/finite_element.py
+++ b/FIAT/finite_element.py
@@ -190,7 +190,7 @@ def tabulate(self, order, points, entity=None):
                      default cell-wise tabulation is performed.
         """
         if entity is None:
-            entity = (self.ref_el.get_spatial_dimension(), 0)
+            entity = (self.ref_el.get_topological_dimension(), 0)
 
         entity_dim, entity_id = entity
         transform = self.ref_el.get_entity_transform(entity_dim, entity_id)
diff --git a/FIAT/functional.py b/FIAT/functional.py
index 9a6dbdf6..acfc9b5a 100644
--- a/FIAT/functional.py
+++ b/FIAT/functional.py
@@ -361,7 +361,7 @@ class IntegralMomentOfNormalDerivative(IntegralMomentOfDerivative):
     def __init__(self, ref_el, facet_no, Q_face, f_at_qpts):
         n = ref_el.compute_normal(facet_no)
         # map points onto facet
-        sd = ref_el.get_spatial_dimension()
+        sd = ref_el.get_topological_dimension()
         Q = quadrature.FacetQuadratureRule(ref_el, sd-1, facet_no, Q_face, avg=True)
         super().__init__(ref_el, Q, f_at_qpts, n, nm="IntegralMomentOfNormalDerivative")
 
diff --git a/FIAT/johnson_mercier.py b/FIAT/johnson_mercier.py
index af70ff0a..3bf36310 100644
--- a/FIAT/johnson_mercier.py
+++ b/FIAT/johnson_mercier.py
@@ -13,13 +13,13 @@ def __init__(self, ref_complex, degree, variant=None, quad_scheme=None):
             raise ValueError(f"Johnson-Mercier does not have the {variant} variant")
         ref_el = ref_complex.get_parent()
         top = ref_el.get_topology()
-        sd = ref_el.get_spatial_dimension()
         entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)}
         nodes = []
 
         # Face dofs: bidirectional (nn and nt) moments
-        dim = sd - 1
-        R = numpy.array([[0, 1], [-1, 0]])
+        tdim = ref_el.get_topological_dimension()
+        dim = tdim - 1
+        n = list(map(ref_el.compute_scaled_normal, sorted(top[dim])))
         ref_facet = ref_el.construct_subelement(dim)
         Qref = parse_quadrature_scheme(ref_facet, 2*degree, quad_scheme)
         P = polynomial_set.ONPolynomialSet(ref_facet, degree)
@@ -28,9 +28,9 @@ def __init__(self, ref_complex, degree, variant=None, quad_scheme=None):
             cur = len(nodes)
             Q = FacetQuadratureRule(ref_el, dim, f, Qref, avg=True)
             thats = ref_el.compute_tangents(dim, f)
-            if sd == 2:
+            if tdim == 2:
                 # Face moments of sigma.n against n P1 and t P1
-                nhat = numpy.dot(R, *thats)
+                nhat = n[f]
                 components = (nhat, *thats)
             else:
                 # Face moments of sigma.n against n P1 and (n x t_j) P1
@@ -44,14 +44,14 @@ def __init__(self, ref_complex, degree, variant=None, quad_scheme=None):
 
         cur = len(nodes)
         # Interior dofs: moments for each independent component
-        n = list(map(ref_el.compute_scaled_normal, sorted(top[sd-1])))
         Q = parse_quadrature_scheme(ref_complex, 2*degree-1, quad_scheme)
         P = polynomial_set.ONPolynomialSet(ref_el, degree-1, scale="L2 piola")
+        sd = ref_el.get_spatial_dimension()
         phis = P.tabulate(Q.get_points())[(0,) * sd]
         nodes.extend(TensorBidirectionalIntegralMoment(ref_el, n[i+1], n[j+1], Q, phi)
-                     for phi in phis for i in range(sd) for j in range(i, sd))
+                     for phi in phis for i in range(tdim) for j in range(i, tdim))
 
-        entity_ids[sd][0].extend(range(cur, len(nodes)))
+        entity_ids[tdim][0].extend(range(cur, len(nodes)))
 
         super().__init__(nodes, ref_el, entity_ids)
 
diff --git a/FIAT/macro.py b/FIAT/macro.py
index 2f2795b1..537a1aad 100644
--- a/FIAT/macro.py
+++ b/FIAT/macro.py
@@ -30,15 +30,18 @@ def xy_to_bary(verts, pts, result=None):
     verts = numpy.asarray(verts)
     pts = numpy.asarray(pts)
     npts = pts.shape[0]
-    sdim = verts.shape[1]
+    nverts = verts.shape[0]
 
-    mat = numpy.vstack((verts.T, numpy.ones((1, sdim+1))))
+    A = numpy.vstack((verts.T, numpy.ones((1, nverts))))
     b = numpy.vstack((pts.T, numpy.ones((1, npts))))
-    foo = numpy.linalg.solve(mat, b)
+    if A.shape[0] == A.shape[1]:
+        bary = numpy.linalg.solve(A, b)
+    else:
+        bary, *_ = numpy.linalg.lstsq(A, b)
     if result is None:
-        return numpy.copy(foo.T)
+        result = bary.T.copy()
     else:
-        result[:, :] = foo[:, :].T
+        numpy.copyto(result, bary.T)
     return result
 
 
@@ -129,7 +132,7 @@ def __init__(self, parent, vertices, topology):
         self._child_to_parent = child_to_parent
         self._parent_to_children = parent_to_children
 
-        sd = parent.get_spatial_dimension()
+        sd = parent.get_topological_dimension()
         inv_top = invert_cell_topology(topology)
 
         # dict mapping cells to their boundary facets for each dimension,
@@ -186,7 +189,7 @@ def construct_subelement(self, dimension):
         return self.get_parent().construct_subelement(dimension)
 
     def get_facet_element(self):
-        dimension = self.get_spatial_dimension()
+        dimension = self.get_topological_dimension()
         return self.construct_subelement(dimension - 1)
 
     def is_macrocell(self):
@@ -256,7 +259,7 @@ class PowellSabinSplit(SplitSimplicialComplex):
     """
     def __init__(self, ref_el, dimension=1):
         self.split_dimension = dimension
-        sd = ref_el.get_spatial_dimension()
+        sd = ref_el.get_topological_dimension()
         top = ref_el.get_topology()
         connectivity = ref_el.get_connectivity()
         new_verts = list(ref_el.get_vertices())
@@ -293,7 +296,7 @@ def construct_subcomplex(self, dimension):
         """Constructs the reference subcomplex of the parent complex
         specified by subcomplex dimension.
         """
-        if dimension == self.get_dimension():
+        if dimension == self.get_topological_dimension():
             return self
         parent = self.get_parent_complex()
         subcomplex = parent.construct_subcomplex(dimension)
@@ -316,7 +319,7 @@ def __new__(cls, ref_el):
             return ref_el._split_cache.setdefault(cls, self)
 
     def __init__(self, ref_el):
-        sd = ref_el.get_spatial_dimension()
+        sd = ref_el.get_topological_dimension()
         super().__init__(ref_el, dimension=sd)
 
 
@@ -332,7 +335,7 @@ def __new__(cls, ref_el):
             return ref_el._split_cache.setdefault(cls, self)
 
     def __init__(self, ref_el):
-        sd = ref_el.get_spatial_dimension()
+        sd = ref_el.get_topological_dimension()
         super().__init__(ref_el, dimension=sd-1)
 
 
@@ -389,7 +392,7 @@ class MacroQuadratureRule(QuadratureRule):
     :returns: A quadrature rule on the sub entities of the simplicial complex.
     """
     def __init__(self, ref_el, Q_ref, parent_facets=None):
-        parent_dim = Q_ref.ref_el.get_spatial_dimension()
+        parent_dim = Q_ref.ref_el.get_topological_dimension()
         if parent_facets is not None:
             parent_to_children = ref_el.get_parent_to_children()
             facets = []
@@ -409,7 +412,7 @@ def __init__(self, ref_el, Q_ref, parent_facets=None):
 
         # Collapse repeated points if any of them lie on facets
         atol = 1E-10
-        sd = ref_el.get_spatial_dimension()
+        sd = ref_el.get_topological_dimension()
         top = ref_el.get_topology()
         for cell in top[sd]:
             bary = ref_el.compute_barycentric_coordinates(pts, entity=(sd, cell))
@@ -450,7 +453,7 @@ def __init__(self, ref_el, degree, order=1, vorder=None, shape=(), **kwargs):
         if not isinstance(order, (int, dict)):
             raise TypeError(f"'order' must be either an int or dict, not {type(order).__name__}")
 
-        sd = ref_el.get_spatial_dimension()
+        sd = ref_el.get_topological_dimension()
         if isinstance(order, int):
             order = {sd-1: dict.fromkeys(ref_el.get_interior_facets(sd-1), order)}
         if vorder is not None:
@@ -524,10 +527,17 @@ def hdiv_conforming_coefficients(U, order=0):
     ref_el = U.get_reference_element()
     coeffs = U.get_coeffs()
     shape = U.get_shape()
+
+    sd = ref_el.get_topological_dimension()
+    if sd != ref_el.get_spatial_dimension():
+        parent = ref_el.get_parent() or ref_el
+        thats = parent.compute_tangents(sd, 0)
+        mapping = "double contravariant piola" if len(shape) == 2 else "contravariant piola"
+        coeffs = pullback(coeffs, mapping, J=thats.T)
+        shape = coeffs.shape[1:-1]
+
     expansion_set = U.get_expansion_set()
     k = 1 if expansion_set.continuity == "C0" else 0
-
-    sd = ref_el.get_spatial_dimension()
     facet_el = ref_el.construct_subelement(sd-1)
 
     phi_deg = 0 if sd == 1 else degree - k
@@ -565,7 +575,7 @@ class HDivPolynomialSet(polynomial_set.PolynomialSet):
     :kwarg scale: The scale for the underlying ExpansionSet.
     """
     def __init__(self, ref_el, degree, order=0, **kwargs):
-        sd = ref_el.get_spatial_dimension()
+        sd = ref_el.get_topological_dimension()
         U = polynomial_set.ONPolynomialSet(ref_el, degree, shape=(sd,), **kwargs)
         coeffs = hdiv_conforming_coefficients(U, order=order)
         super().__init__(ref_el, degree, degree, U.expansion_set, coeffs)
@@ -613,9 +623,15 @@ def pullback(phi, mapping, J=None, Jinv=None, Jdet=None):
     if J is None:
         J = numpy.linalg.pinv(Jinv)
     if Jinv is None:
-        Jinv = numpy.linalg.pinv(J)
+        if J.shape[0] == J.shape[1]:
+            Jinv = numpy.linalg.inv(J)
+        else:
+            Jinv = numpy.linalg.pinv(J)
     if Jdet is None:
-        Jdet = numpy.linalg.det(J)
+        if J.shape[0] == J.shape[1]:
+            Jdet = numpy.linalg.det(J)
+        else:
+            Jdet = numpy.linalg.det(J.T @ J)**0.5
     F1 = Jinv.T
     F2 = J / Jdet
 
@@ -643,7 +659,7 @@ class MacroPolynomialSet(polynomial_set.PolynomialSet):
     """
     def __init__(self, ref_el, element):
         top = ref_el.get_topology()
-        sd = ref_el.get_spatial_dimension()
+        sd = ref_el.get_topological_dimension()
 
         mapping, = set(element.mapping())
         base_ref_el = element.get_reference_element()
diff --git a/FIAT/mardal_tai_winther.py b/FIAT/mardal_tai_winther.py
index 54e65cdb..be05a388 100644
--- a/FIAT/mardal_tai_winther.py
+++ b/FIAT/mardal_tai_winther.py
@@ -15,6 +15,7 @@
 from FIAT.nedelec import Nedelec
 from FIAT.quadrature import FacetQuadratureRule
 from FIAT.quadrature_schemes import create_quadrature
+from FIAT.reference_element import make_affine_mapping
 
 
 def curl(tabulation):
@@ -37,20 +38,25 @@ def curl(tabulation):
 
 def MardalTaiWintherSpace(ref_el, order=1):
     """Construct the MTW space BDM(order) + curl(B [P1]^d)."""
+    orig_ref_el = ref_el
     sd = ref_el.get_spatial_dimension()
-    k = sd + 1
+    tdim = ref_el.get_topological_dimension()
+    if tdim != sd:
+        ref_el = ref_el.construct_subelement(tdim)
+
+    k = tdim + 1
     assert order < k
-    # [Pk]^d = vector polynomials of degree k = sd+1
-    Pk = polynomial_set.ONPolynomialSet(ref_el, k, shape=(sd,), scale="orthonormal")
+    # [Pk]^d = vector polynomials of degree k = tdim+1
+    Pk = polynomial_set.ONPolynomialSet(ref_el, k, shape=(tdim,), scale="orthonormal")
 
     # Grab BDM(order) = [P_order]^d from [Pk]^d
     dimP1 = expansions.polynomial_dimension(ref_el, order)
     dimPk = expansions.polynomial_dimension(ref_el, k)
-    ids = [i+dimPk*j for i in range(dimP1) for j in range(sd)]
+    ids = [i+dimPk*j for i in range(dimP1) for j in range(tdim)]
     BDM = Pk.take(ids)
 
     # Project curl(B [P1]^d) into [Pk]^d
-    shape = () if sd == 2 else ((sd*(sd-1))//2,)
+    shape = () if tdim == 2 else ((tdim*(tdim-1))//2,)
     BP1 = polynomial_set.make_bubbles(ref_el, k+1, shape=shape)
 
     Q = create_quadrature(ref_el, 2*k)
@@ -60,21 +66,31 @@ def MardalTaiWintherSpace(ref_el, order=1):
     BP1_at_qpts = BP1.tabulate(qpts, 1)
 
     inner = lambda u, v, qwts: numpy.tensordot(u, numpy.multiply(v, qwts), axes=(range(1, u.ndim),)*2)
-    C = inner(curl(BP1_at_qpts), Pk_at_qpts[(0,)*sd], qwts)
+    C = inner(curl(BP1_at_qpts), Pk_at_qpts[(0,)*tdim], qwts)
     coeffs = numpy.tensordot(C, Pk.get_coeffs(), axes=(1, 0))
     curlBP1 = polynomial_set.PolynomialSet(ref_el, k, k, Pk.get_expansion_set(), coeffs)
 
-    return polynomial_set.polynomial_set_union_normalized(BDM, curlBP1)
+    MTW = polynomial_set.polynomial_set_union_normalized(BDM, curlBP1)
+
+    if tdim != sd:
+        # Manifold case
+        J, b = make_affine_mapping(ref_el.vertices, orig_ref_el.vertices)
+        coeffs = MTW.get_coeffs()
+        coeffs = numpy.tensordot(J, coeffs.transpose(1, 0, 2), (1, 0)).transpose(1, 0, 2)
+        expansion_set = Pk.get_expansion_set().reconstruct(ref_el=orig_ref_el)
+        MTW = polynomial_set.PolynomialSet(orig_ref_el, k, k, expansion_set, coeffs)
+    return MTW
 
 
 class MardalTaiWintherDual(dual_set.DualSet):
     """Degrees of freedom for Mardal-Tai-Winther elements."""
     def __init__(self, ref_el, order, quad_scheme):
+        tdim = ref_el.get_topological_dimension()
         sd = ref_el.get_spatial_dimension()
         top = ref_el.get_topology()
         entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top}
         nodes = []
-        degree = sd + 1
+        degree = tdim + 1
 
         # On each facet, let n be its normal.  We need to integrate
         # u.n against a Dubiner basis for P1
@@ -82,33 +98,35 @@ def __init__(self, ref_el, order, quad_scheme):
         ref_facet = ref_el.get_facet_element()
         Q = parse_quadrature_scheme(ref_facet, degree+order, quad_scheme)
 
+        fdim = ref_facet.get_spatial_dimension()
         P1 = polynomial_set.ONPolynomialSet(ref_facet, order)
-        P1_at_qpts = P1.tabulate(Q.get_points())[(0,)*(sd - 1)]
-        if sd == 2:
+        P1_at_qpts = P1.tabulate(Q.get_points())[(0,)*fdim]
+        if tdim == 2:
             # For 2D just take the constant
             RT_at_qpts = P1_at_qpts[:1, None, :]
         else:
             # Basis for lowest-order RT [(1, 0), (0, 1), (x, y)]
-            RT_at_qpts = numpy.zeros((3, sd-1, P1_at_qpts.shape[-1]))
+            RT_at_qpts = numpy.zeros((3, tdim-1, P1_at_qpts.shape[-1]))
             RT_at_qpts[0, 0, :] = P1_at_qpts[0, None, :]
             RT_at_qpts[1, 1, :] = P1_at_qpts[0, None, :]
             RT_at_qpts[2, 0, :] = P1_at_qpts[1, None, :]
             RT_at_qpts[2, 1, :] = P1_at_qpts[2, None, :]
 
-        for f in sorted(top[sd-1]):
+        for f in sorted(top[tdim-1]):
             cur = len(nodes)
             n = ref_el.compute_scaled_normal(f)
-            Qf = FacetQuadratureRule(ref_el, sd-1, f, Q, avg=True)
+            Qf = FacetQuadratureRule(ref_el, tdim-1, f, Q, avg=True)
             # Normal moments against P_{order}
             nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, numpy.outer(n, phi)) for phi in P1_at_qpts)
             # Map the RT basis into the facet
+
             Jf = Qf.jacobian()
             phis = numpy.tensordot(Jf, RT_at_qpts.transpose(1, 0, 2), (1, 0)).transpose(1, 0, 2)
-            if sd == 3:
+            if tdim == 3:
                 # Moments against cross(n, RT)
                 phis = numpy.cross(n[None, :, None], phis, axis=1)
             nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, phi) for phi in phis)
-            entity_ids[sd-1][f].extend(range(cur, len(nodes)))
+            entity_ids[tdim-1][f].extend(range(cur, len(nodes)))
 
         # Interior nodes: moments against Nedelec(order-1)
         if order > 1:
@@ -117,7 +135,7 @@ def __init__(self, ref_el, order, quad_scheme):
             phis = Ned.tabulate(0, Q.get_points())[(0,) * sd]
             cur = len(nodes)
             nodes.extend(FrobeniusIntegralMoment(ref_el, Q, phi) for phi in phis)
-            entity_ids[sd][0] = list(range(cur, len(nodes)))
+            entity_ids[tdim][0] = list(range(cur, len(nodes)))
 
         super().__init__(nodes, ref_el, entity_ids)
 
@@ -126,7 +144,7 @@ class MardalTaiWinther(finite_element.CiarletElement):
     """The definition of the Mardal-Tai-Winther element.
     """
     def __init__(self, ref_el, order=1, quad_scheme=None):
-        sd = ref_el.get_spatial_dimension()
+        sd = ref_el.get_topological_dimension()
         if sd not in (2, 3):
             raise ValueError(f"{type(self).__name__} only defined in dimension 2 and 3.")
         if not ref_el.is_simplex():
diff --git a/FIAT/morley.py b/FIAT/morley.py
index af09c52c..d8a38356 100644
--- a/FIAT/morley.py
+++ b/FIAT/morley.py
@@ -19,7 +19,7 @@ class MorleyDualSet(dual_set.DualSet):
 
     def __init__(self, ref_el, degree):
         top = ref_el.get_topology()
-        sd = ref_el.get_spatial_dimension()
+        sd = ref_el.get_topological_dimension()
         entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top}
         nodes = []
 
diff --git a/FIAT/polynomial_set.py b/FIAT/polynomial_set.py
index b3aee67f..bbc45864 100644
--- a/FIAT/polynomial_set.py
+++ b/FIAT/polynomial_set.py
@@ -224,7 +224,7 @@ class ONSymTensorPolynomialSet(PolynomialSet):
     def __init__(self, ref_el, degree, size=None, **kwargs):
         expansion_set = expansions.ExpansionSet(ref_el, **kwargs)
 
-        sd = ref_el.get_spatial_dimension()
+        sd = ref_el.get_topological_dimension()
         if size is None:
             size = sd
 
@@ -256,7 +256,7 @@ class TracelessTensorPolynomialSet(PolynomialSet):
     def __init__(self, ref_el, degree, size=None, **kwargs):
         expansion_set = expansions.ExpansionSet(ref_el, **kwargs)
 
-        sd = ref_el.get_spatial_dimension()
+        sd = ref_el.get_topological_dimension()
         if size is None:
             size = sd
 
@@ -286,11 +286,11 @@ def make_bubbles(ref_el, degree, codim=0, shape=(), scale="L2 piola"):
     """Construct a polynomial set with codim bubbles up to the given degree.
     """
     poly_set = ONPolynomialSet(ref_el, degree, shape=shape, scale=scale, variant="bubble")
-    if ref_el.get_spatial_dimension() == 0:
+    if ref_el.get_topological_dimension() == 0:
         return poly_set
 
     entity_ids = expansions.polynomial_entity_ids(ref_el, degree, continuity="C0")
-    sd = ref_el.get_spatial_dimension()
+    sd = ref_el.get_topological_dimension()
     dim = sd - codim
     indices = list(chain(*entity_ids[dim].values()))
     if shape != ():
diff --git a/FIAT/quadrature.py b/FIAT/quadrature.py
index 7cf3159c..4a8c8c9d 100644
--- a/FIAT/quadrature.py
+++ b/FIAT/quadrature.py
@@ -174,7 +174,7 @@ class CollapsedQuadratureSimplexRule(QuadratureRule):
     from the hypercube to the simplex."""
 
     def __init__(self, ref_el, m):
-        dim = ref_el.get_spatial_dimension()
+        dim = ref_el.get_topological_dimension()
         Ref1 = reference_element.default_simplex(dim)
         pts_ref, wts_ref = simplexgausslegendre(dim, m)
         pts, wts = map_quadrature(pts_ref, wts_ref, Ref1, ref_el)
diff --git a/FIAT/quadrature_schemes.py b/FIAT/quadrature_schemes.py
index ca0c478a..20ac7682 100644
--- a/FIAT/quadrature_schemes.py
+++ b/FIAT/quadrature_schemes.py
@@ -69,7 +69,7 @@ def create_quadrature(ref_el, degree, scheme="default", entity=None):
         return FacetQuadratureRule(ref_el, dimension, entity_id, Q_ref)
 
     if ref_el.is_macrocell():
-        dimension = ref_el.get_dimension()
+        dimension = ref_el.get_topological_dimension()
         sub_el = ref_el.construct_subelement(dimension)
         Q_ref = create_quadrature(sub_el, degree, scheme=scheme)
         return MacroQuadratureRule(ref_el, Q_ref)
@@ -119,7 +119,7 @@ def _fiat_scheme(ref_el, degree):
 def _kmv_lump_scheme(ref_el, degree):
     """Specialized quadrature schemes for P < 6 for KMV simplical elements."""
 
-    sd = ref_el.get_spatial_dimension()
+    sd = ref_el.get_topological_dimension()
     if sd == 1:
         num_points = degree + 1
         return GaussLobattoLegendreQuadratureLineRule(ref_el, num_points)
@@ -334,7 +334,7 @@ def xg_scheme(ref_el, degree):
         Mathematics with Applications, vol. 59, no. 2, pp. 663-676, 2010.
         http://dx.doi.org/10.1016/j.camwa.2009.10.027
     """
-    dim = ref_el.get_spatial_dimension()
+    dim = ref_el.get_topological_dimension()
     if dim == 2:
         from FIAT.xg_quad_data import triangle_table as table
     elif dim == 3:
diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py
index 87794b50..b896059e 100644
--- a/FIAT/reference_element.py
+++ b/FIAT/reference_element.py
@@ -218,6 +218,10 @@ def get_topology(self):
         and each value is a dictionary mapping."""
         return self.topology
 
+    def get_topological_dimension(self):
+        """Returns the topological dimension of the element."""
+        return max(self.topology)
+
     def get_connectivity(self):
         """Returns a dictionary encoding the connectivity of the element.
 
@@ -399,8 +403,8 @@ def __init__(self, shape, vertices, topology):
     def compute_normal(self, facet_i, cell=None):
         """Returns the unit normal vector to facet i of codimension 1."""
 
-        t = self.get_topology()
-        sd = self.get_spatial_dimension()
+        top = self.get_topology()
+        sd = self.get_topological_dimension()
 
         # To handle simplicial complex case:
         # Find a subcell of which facet_i is on the boundary
@@ -409,28 +413,22 @@ def compute_normal(self, facet_i, cell=None):
         if cell is None:
             cell = next(k for k, facets in enumerate(self.connectivity[(sd, sd-1)])
                         if facet_i in facets)
-        verts = numpy.asarray(self.get_vertices_of_subcomplex(t[sd][cell]))
+
         # Interval case
-        if self.get_shape() == LINE:
-            v_i = t[1][cell].index(t[0][facet_i][0])
+        if sd == 1:
+            verts = numpy.array(self.get_vertices_of_subcomplex(top[sd][cell]))
+            v_i = top[sd][cell].index(top[sd-1][facet_i][0])
             n = verts[v_i] - verts[[1, 0][v_i]]
             return n / numpy.linalg.norm(n)
 
-        # vectors from vertex 0 to each other vertex.
-        vert_vecs_from_v0 = verts[1:, :] - verts[:1, :]
-
-        (u, s, _) = numpy.linalg.svd(vert_vecs_from_v0)
+        # Orthogonalize the set of vectors that span the cell
+        cell_span = self.compute_tangents(sd, cell)
+        (_, s, vt) = numpy.linalg.svd(cell_span)
         rank = len([si for si in s if si > 1.e-10])
-
-        # this is the set of vectors that span the simplex
-        spanu = u[:, :rank]
-
-        vert_coords_of_facet = \
-            self.get_vertices_of_subcomplex(t[sd-1][facet_i])
+        cell_space = numpy.transpose(vt[:rank, :])
 
         # now I find everything normal to the facet.
-        vcf = numpy.asarray(vert_coords_of_facet)
-        facet_span = vcf[1:, :] - vcf[:1, :]
+        facet_span = self.compute_tangents(sd-1, facet_i)
         (_, sf, vft) = numpy.linalg.svd(facet_span)
 
         # now get the null space from vft
@@ -438,24 +436,22 @@ def compute_normal(self, facet_i, cell=None):
         facet_normal_space = numpy.transpose(vft[rankfacet:, :])
 
         # now, I have to compute the intersection of
-        # facet_span with facet_normal_space
-        foo = linalg_subspace_intersection(facet_normal_space, spanu)
-
-        num_cols = foo.shape[1]
-
-        if num_cols != 1:
-            raise Exception("barf in normal computation")
+        # cell_space with facet_normal_space
+        normal = linalg_subspace_intersection(facet_normal_space, cell_space)
+        if normal.shape[1] != 1:
+            raise RuntimeError(f"Found a normal space of dimension {normal.shape[1]}")
 
         # now need to get the correct sign
         # get a vector in the direction
-        nfoo = foo[:, 0]
+        normal = normal.flatten()
 
         # what is the vertex not in the facet?
-        verts_set = set(t[sd][cell])
-        verts_facet = set(t[sd - 1][facet_i])
-        verts_diff = verts_set.difference(verts_facet)
+        verts_set = set(top[sd][cell])
+        verts_facet = set(top[sd - 1][facet_i])
+        verts_diff = verts_set - verts_facet
         if len(verts_diff) != 1:
-            raise Exception("barf in normal computation: getting sign")
+            raise ValueError(f"Facet {facet_i} is not contained by the cell {cell}. "
+                             "Unable to compute sign of outward-facing normal.")
         vert_off = verts_diff.pop()
         vert_on = verts_facet.pop()
 
@@ -463,10 +459,9 @@ def compute_normal(self, facet_i, cell=None):
         v_to_facet = numpy.array(self.vertices[vert_on]) \
             - numpy.array(self.vertices[vert_off])
 
-        if numpy.dot(v_to_facet, nfoo) > 0.0:
-            return nfoo
-        else:
-            return -nfoo
+        if numpy.dot(v_to_facet, normal) < 0.0:
+            normal *= -1
+        return normal
 
     def compute_tangents(self, dim, i):
         """Computes tangents in any dimension based on differences
@@ -540,7 +535,7 @@ def make_points(self, dim, entity_id, order, variant=None, interior=1):
     def volume(self):
         """Computes the volume of the simplicial complex in the appropriate
         dimensional measure."""
-        sd = self.get_spatial_dimension()
+        sd = self.get_topological_dimension()
         return sum(self.volume_of_subcomplex(sd, k)
                    for k in self.topology[sd])
 
@@ -551,8 +546,12 @@ def volume_of_subcomplex(self, dim, facet_no):
     def compute_scaled_normal(self, facet_i):
         """Returns the unit normal to facet_i of scaled by the
         volume of that facet."""
-        dim = self.get_spatial_dimension()
-        if dim == 2:
+        dim = self.get_topological_dimension()
+        if self.get_spatial_dimension() == 3 and dim == 2:
+            R = numpy.cross(*self.compute_tangents(dim, 0))
+            R *= -1/numpy.linalg.norm(R)
+            return numpy.cross(R, *self.compute_tangents(dim-1, facet_i))
+        elif dim == 2:
             n, = self.compute_tangents(dim-1, facet_i)
             n[0], n[1] = n[1], -n[0]
             return n
@@ -575,7 +574,7 @@ def get_entity_transform(self, dim, entity):
         :arg entity: entity number (integer)
         """
         topology = self.get_topology()
-        celldim = self.get_spatial_dimension()
+        celldim = self.get_topological_dimension()
         codim = celldim - dim
         if dim == 0:
             # Special case vertices.
@@ -587,7 +586,7 @@ def get_entity_transform(self, dim, entity):
             return lambda x: x
         else:
             subcell = self.construct_subelement(dim)
-            subdim = subcell.get_spatial_dimension()
+            subdim = subcell.get_topological_dimension()
             assert subdim == celldim - codim
 
             # Entity vertices in entity space.
@@ -618,10 +617,10 @@ def compute_barycentric_coordinates(self, points, entity=None, rescale=False):
         if len(points) == 0:
             return points
         if entity is None:
-            entity = (self.get_spatial_dimension(), 0)
+            entity = (self.get_topological_dimension(), 0)
         entity_dim, entity_id = entity
         top = self.get_topology()
-        sd = self.get_spatial_dimension()
+        sd = self.get_topological_dimension()
 
         # get a subcell containing the entity and the restriction indices of the entity
         indices = slice(None)
@@ -922,7 +921,7 @@ def cell_orientation_reflection_map(self):
         return make_cell_orientation_reflection_map_simplex(self.get_dimension())
 
     def get_facet_element(self):
-        dimension = self.get_spatial_dimension()
+        dimension = self.get_topological_dimension()
         return self.construct_subelement(dimension - 1)
 
 
@@ -1033,7 +1032,11 @@ def __init__(self):
     def compute_normal(self, i):
         "UFC consistent normal"
         t = self.compute_tangents(1, i)[0]
-        n = numpy.array((t[1], -t[0]))
+        if self.get_spatial_dimension() == 2:
+            n = numpy.array((t[1], -t[0]))
+        else:
+            R = -numpy.cross(*self.compute_tangents(2, 0))
+            n = numpy.cross(R, t)
         return n / numpy.linalg.norm(n)
 
 
@@ -1609,38 +1612,20 @@ def __init__(self):
 
 def make_affine_mapping(xs, ys):
     """Constructs (A,b) such that x --> A * x + b is the affine
-    mapping from the simplex defined by xs to the simplex defined by ys."""
-
-    dim_x = len(xs[0])
-    dim_y = len(ys[0])
-
-    if len(xs) != len(ys):
-        raise Exception("")
-
-    # find A in R^{dim_y,dim_x}, b in R^{dim_y} such that
-    # A xs[i] + b = ys[i] for all i
-
-    mat = numpy.zeros((dim_x * dim_y + dim_y, dim_x * dim_y + dim_y), "d")
-    rhs = numpy.zeros((dim_x * dim_y + dim_y,), "d")
-
-    # loop over points
-    for i in range(len(xs)):
-        # loop over components of each A * point + b
-        for j in range(dim_y):
-            row_cur = i * dim_y + j
-            col_start = dim_x * j
-            col_finish = col_start + dim_x
-            mat[row_cur, col_start:col_finish] = numpy.array(xs[i])
-            rhs[row_cur] = ys[i][j]
-            # need to get terms related to b
-            mat[row_cur, dim_y * dim_x + j] = 1.0
-
-    sol = numpy.linalg.solve(mat, rhs)
-
-    A = numpy.reshape(sol[:dim_x * dim_y], (dim_y, dim_x))
-    b = sol[dim_x * dim_y:]
-
-    return A, b
+    mapping from the simplex defined by xs to the simplex defined by ys.
+    Uses least-squares when the simplex dimension does not match the spatial
+    dimension.
+    """
+    X = numpy.asarray(xs)
+    Y = numpy.asarray(ys)
+    DX = X[1:] - X[:1]
+    DY = Y[1:] - Y[:1]
+    if DX.shape[0] == DX.shape[1]:
+        AT = numpy.linalg.solve(DX, DY)
+    else:
+        AT, *_ = numpy.linalg.lstsq(DX, DY)
+    b = Y[0] - numpy.dot(X[0], AT)
+    return AT.T, b
 
 
 def ufc_hypercube(spatial_dim):
diff --git a/finat/morley.py b/finat/morley.py
index ca3e89e8..aafefa8c 100644
--- a/finat/morley.py
+++ b/finat/morley.py
@@ -9,8 +9,8 @@
 
 
 def morley_transform(cell, J, detJ, face):
-    adjugate = lambda A: ListTensor([[A[1, 1], -1*A[1, 0]], [-1*A[0, 1], A[0, 0]]])
-    sd = cell.get_spatial_dimension()
+    adjugate = lambda A: ListTensor([[A[1, 1], -A[1, 0]], [-A[0, 1], A[0, 0]]])
+    sd = cell.get_topological_dimension()
     thats = cell.compute_tangents(sd-1, face)
     nhat = numpy.cross(*thats)
     ahat = numpy.linalg.norm(nhat)
@@ -38,7 +38,7 @@ def __init__(self, cell, degree=2):
         super().__init__(FIAT.Morley(cell, degree=degree))
 
     def basis_transformation(self, coordinate_mapping):
-        sd = self.cell.get_spatial_dimension()
+        sd = self.cell.get_topological_dimension()
         top = self.cell.get_topology()
         # Jacobians at barycenter
         bary, = self.cell.make_points(sd, 0, sd+1)
diff --git a/finat/mtw.py b/finat/mtw.py
index 3e79c5d0..3750d221 100644
--- a/finat/mtw.py
+++ b/finat/mtw.py
@@ -17,7 +17,7 @@ def __init__(self, cell, order=1):
         super().__init__(FIAT.MardalTaiWinther(cell, order=order))
 
     def basis_transformation(self, coordinate_mapping):
-        sd = self.cell.get_spatial_dimension()
+        sd = self.cell.get_topological_dimension()
         bary, = self.cell.make_points(sd, 0, sd+1)
         J = coordinate_mapping.jacobian_at(bary)
         detJ = coordinate_mapping.detJ_at(bary)
diff --git a/finat/piola_mapped.py b/finat/piola_mapped.py
index 447237e2..19444ad3 100644
--- a/finat/piola_mapped.py
+++ b/finat/piola_mapped.py
@@ -10,8 +10,8 @@
 def piola_inverse(fiat_cell, J, detJ):
     """Return the basis transformation of evaluation at a point.
     This simply inverts the Piola transform inv(J / detJ) = adj(J)."""
-    sd = fiat_cell.get_spatial_dimension()
-    Jnp = numpy.array([[J[i, j] for j in range(sd)] for i in range(sd)])
+    gdim, tdim = J.shape
+    Jnp = numpy.array([[J[i, j] for j in range(tdim)] for i in range(gdim)])
     return adjugate(Jnp)
 
 
@@ -27,7 +27,7 @@ def normal_tangential_edge_transform(fiat_cell, J, detJ, f):
     alpha = Jn @ Jt
     beta = Jt @ Jt
     # Compute the last row of inv([[1, 0], [alpha/detJ, beta/detJ]])
-    return (-1 * alpha / beta, detJ / beta)
+    return (-alpha / beta, detJ / beta)
 
 
 def normal_tangential_face_transform(fiat_cell, J, detJ, f):
@@ -65,7 +65,7 @@ def normal_tangential_transform(fiat_cell, J, detJ, f):
         Bnt is the numpy.ndarray of normal-tangential coefficients, and
         Btt is the tangential-tangential coefficient.
     """
-    if fiat_cell.get_spatial_dimension() == 2:
+    if fiat_cell.get_topological_dimension() == 2:
         return normal_tangential_edge_transform(fiat_cell, J, detJ, f)
     else:
         return normal_tangential_face_transform(fiat_cell, J, detJ, f)
@@ -81,7 +81,7 @@ def __init__(self, fiat_element):
 
         # On each facet we expect the normal dof followed by the tangential ones
         # The tangential dofs should be numbered last, and are constrained to be zero
-        sd = self.cell.get_spatial_dimension()
+        sd = self.cell.get_topological_dimension()
         reduced_dofs = deepcopy(self._element.entity_dofs())
         reduced_dim = 0
         cur = reduced_dofs[sd-1][0][0]
@@ -99,7 +99,7 @@ def space_dimension(self):
         return self._space_dimension
 
     def basis_transformation(self, coordinate_mapping):
-        sd = self.cell.get_spatial_dimension()
+        sd = self.cell.get_topological_dimension()
         bary, = self.cell.make_points(sd, 0, sd+1)
         J = coordinate_mapping.jacobian_at(bary)
         detJ = coordinate_mapping.detJ_at(bary)
diff --git a/test/finat/conftest.py b/test/finat/conftest.py
index c31ef735..1fab8dd3 100644
--- a/test/finat/conftest.py
+++ b/test/finat/conftest.py
@@ -20,7 +20,12 @@ def cell_size(self):
         return np.ones((len(self.ref_cell.vertices),))
 
     def detJ_at(self, point):
-        return gem.Literal(np.linalg.det(self.A))
+        shape = self.A.shape
+        if shape[0] == shape[1]:
+            detA = np.linalg.det(self.A)
+        else:
+            detA = np.linalg.det(self.A.T @ self.A)**0.5
+        return gem.Literal(detA)
 
     def jacobian_at(self, point):
         return gem.Literal(self.A)
@@ -32,14 +37,14 @@ def normalized_reference_edge_tangents(self):
                         for i in sorted(top[1])]))
 
     def reference_normals(self):
-        sd = self.ref_cell.get_spatial_dimension()
+        sd = self.phys_cell.get_topological_dimension()
         top = self.ref_cell.get_topology()
         return gem.Literal(
             np.asarray([self.ref_cell.compute_normal(i)
                         for i in sorted(top[sd-1])]))
 
     def physical_normals(self):
-        sd = self.phys_cell.get_spatial_dimension()
+        sd = self.phys_cell.get_topological_dimension()
         top = self.phys_cell.get_topology()
         return gem.Literal(
             np.asarray([self.phys_cell.compute_normal(i)
@@ -108,18 +113,19 @@ def ref_el():
 
 @pytest.fixture
 def phys_el():
-    K = {dim: FIAT.ufc_simplex(dim) for dim in (2, 3)}
+    K = {dim: FIAT.ufc_simplex(int(dim)) for dim in (2, 2.5, 3)}
     K[2].vertices = ((0.0, 0.1), (1.17, -0.09), (0.15, 1.84))
     K[3].vertices = ((0, 0, 0),
                      (1., 0.1, -0.37),
                      (0.01, 0.987, -.23),
                      (-0.1, -0.2, 1.38))
+    K[2.5].vertices = K[3].vertices[:-1]
     return K
 
 
 @pytest.fixture
 def ref_to_phys(ref_el, phys_el):
-    return {dim: MyMapping(ref_el[dim], phys_el[dim]) for dim in ref_el}
+    return {dim: MyMapping(ref_el[int(dim)], phys_el[dim]) for dim in phys_el}
 
 
 @pytest.fixture
diff --git a/test/finat/test_zany_mapping.py b/test/finat/test_zany_mapping.py
index 9220dca9..fe74d2ad 100644
--- a/test/finat/test_zany_mapping.py
+++ b/test/finat/test_zany_mapping.py
@@ -14,7 +14,7 @@ def make_unisolvent_points(element, interior=False):
     top = ref_complex.get_topology()
     pts = []
     if interior:
-        dim = ref_complex.get_spatial_dimension()
+        dim = ref_complex.get_topological_dimension()
         for entity in top[dim]:
             pts.extend(ref_complex.make_points(dim, entity, degree+dim+1, variant="gll"))
     else:
@@ -33,12 +33,13 @@ def check_zany_mapping(element, ref_to_phys, *args, **kwargs):
     ref_element = finat_element._element
     ref_cell = ref_element.get_reference_element()
     phys_cell = phys_element.get_reference_element()
-    sd = ref_cell.get_spatial_dimension()
 
     shape = ref_element.value_shape()
+    sd = ref_cell.get_spatial_dimension()
     ref_pts = make_unisolvent_points(ref_element, interior=True)
     ref_vals = ref_element.tabulate(0, ref_pts)[(0,)*sd]
 
+    sd = phys_cell.get_spatial_dimension()
     phys_pts = make_unisolvent_points(phys_element, interior=True)
     phys_vals = phys_element.tabulate(0, phys_pts)[(0,)*sd]
 
@@ -51,16 +52,31 @@ def check_zany_mapping(element, ref_to_phys, *args, **kwargs):
                                                           phys_cell.vertices)
         K = []
         if "covariant" in mapping:
-            K.append(np.linalg.inv(J).T)
+            if J.shape[0] == J.shape[1]:
+                JinvT = np.linalg.inv(J).T
+            else:
+                JinvT = J @ np.linalg.inv(J.T @ J)
+            K.append(JinvT)
         if "contravariant" in mapping:
-            K.append(J / np.linalg.det(J))
+            if J.shape[0] == J.shape[1]:
+                detJ = np.linalg.det(J)
+            else:
+                detJ = np.linalg.det(J.T @ J)**0.5
+            K.append(J / detJ)
 
         if len(shape) == 2:
             piola_map = lambda x: K[0] @ x @ K[-1].T
         else:
             piola_map = lambda x: K[0] @ x
 
-        ref_vals_piola = np.zeros(ref_vals.shape)
+        shp = list(ref_vals.shape)
+        if J.shape[0] != J.shape[1]:
+            if len(shape) == 2:
+                shp[1:-1] = [K[i].shape[0] for i in (0, -1)]
+            elif len(shape) == 1:
+                shp[1:-1] = [K[i].shape[0] for i in (0,)]
+
+        ref_vals_piola = np.zeros(shp)
         for i in range(ref_vals.shape[0]):
             for k in range(ref_vals.shape[-1]):
                 ref_vals_piola[i, ..., k] = piola_map(ref_vals[i, ..., k])
@@ -178,6 +194,16 @@ def test_piola(ref_to_phys, element, dimension):
     check_zany_mapping(element, ref_to_phys[dimension])
 
 
+@pytest.mark.parametrize("element", [
+    finat.MardalTaiWinther,
+    finat.JohnsonMercier,
+    finat.Morley,
+])
+def test_piola_manifold(ref_to_phys, element):
+    dimension = 2.5
+    check_zany_mapping(element, ref_to_phys[dimension])
+
+
 @pytest.mark.parametrize("dimension, element, degree", [
     (3, finat.MardalTaiWinther, 2),
     (3, finat.GuzmanNeilanFirstKindH1, 2),