diff --git a/README.md b/README.md
index 582398abd..c640e7049 100644
--- a/README.md
+++ b/README.md
@@ -212,12 +212,16 @@ with respect to documented and/or tested features.
 ### Unreleased
 
 - Removed: Python 3.7 support
+- Removed: `MappingMortar` and `MortarFacetBasis` in favor of
+  `skfem.supermeshing`
 - Added: Python 3.12 support
 - Added: `Mesh.load` supports new keyword arguments
   `ignore_orientation=True` and `ignore_interior_facets=True` which
   will both speed up the loading of larger three-dimensional meshes by
   ignoring facet orientation and all tags not on the boundary,
   respectively.
+- Added: `skfem.supermeshing` (requires `shapely>=2`) for creating quadrature
+  rules for interpolating between two 1D or 2D meshes.
 - Fixed: `MeshTet` uniform refine was reindexing subdomains incorrectly
 
 ## [8.1.0] - 2023-06-16
diff --git a/docs/examples/ex04.py b/docs/examples/ex04.py
index 65848d250..bb29fcf61 100644
--- a/docs/examples/ex04.py
+++ b/docs/examples/ex04.py
@@ -8,7 +8,7 @@
 
 import numpy as np
 from skfem import *
-from skfem.experimental.supermeshing import intersect, elementwise_quadrature
+from skfem.supermeshing import intersect, elementwise_quadrature
 from skfem.models.elasticity import (linear_elasticity, lame_parameters,
                                      linear_stress)
 from skfem.helpers import dot, sym_grad, jump, mul
diff --git a/docs/examples/ex47.py b/docs/examples/ex47.py
index 55bf0bb53..742a89c9e 100644
--- a/docs/examples/ex47.py
+++ b/docs/examples/ex47.py
@@ -2,7 +2,7 @@
 
 import numpy as np
 from skfem import *
-from skfem.experimental.supermeshing import intersect, elementwise_quadrature
+from skfem.supermeshing import intersect, elementwise_quadrature
 
 
 m1 = MeshLine(np.linspace(1, 10, 20))
diff --git a/docs/examples/ex49.py b/docs/examples/ex49.py
new file mode 100644
index 000000000..d5c6e74db
--- /dev/null
+++ b/docs/examples/ex49.py
@@ -0,0 +1,47 @@
+"""Projection between two meshes using supermesh."""
+
+import numpy as np
+from skfem import *
+from skfem.supermeshing import intersect, elementwise_quadrature
+
+
+m1 = MeshTri.init_tensor(np.linspace(0, 1, 4),
+                         np.linspace(0, 1, 4))
+m2 = MeshQuad().refined(2)
+e1 = ElementTriP2()
+e2 = ElementQuad2()
+
+m12, t1, t2 = intersect(m1, m2)
+
+bases = [
+    Basis(m1, e1),
+    Basis(m2, e2),
+]
+projbases = [
+    Basis(m1, e1, quadrature=elementwise_quadrature(m1, m12, t1, intorder=4), elements=t1),
+    Basis(m2, e2, quadrature=elementwise_quadrature(m2, m12, t2, intorder=4), elements=t2),
+]
+
+
+@BilinearForm
+def mass(u, v, _):
+    return u * v
+
+
+P = mass.assemble(*projbases)
+M = mass.assemble(projbases[1])
+
+y1 = bases[0].project(lambda x: x[0] ** 1.6)
+y2 = solve(M, P.dot(y1))
+
+xs = np.linspace(0, 1, 100)
+xs = np.vstack((xs, xs))
+l2 = (bases[0].interpolator(y1)(xs)
+      - bases[1].interpolator(y2)(xs)).sum()
+
+if __name__ == "__main__":
+    print('L2 error: {}'.format(l2))
+    ax = bases[0].plot(y1, colorbar=True, shading='gouraud')
+    m1.draw(ax=ax, color='ko')
+    ax = bases[1].plot(y2, colorbar=True, shading='gouraud')
+    m2.draw(ax=ax, color='ro').show()
diff --git a/requirements.txt b/requirements.txt
index 05c4e836b..3c741f1d6 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -10,3 +10,4 @@ flake8
 sphinx~=6.1
 autograd
 pep517
+shapely
diff --git a/skfem/__init__.py b/skfem/__init__.py
index cfcc3f3b5..7ed63f3c6 100644
--- a/skfem/__init__.py
+++ b/skfem/__init__.py
@@ -16,7 +16,6 @@
 __all__ = all_mesh + all_assembly + all_element + [  # noqa
     'MappingAffine',
     'MappingIsoparametric',
-    'MappingMortar',  # TODO remove due to deprecation
     'adaptive_theta',
     'build_pc_ilu',
     'build_pc_diag',
diff --git a/skfem/assembly/__init__.py b/skfem/assembly/__init__.py
index b21dd5072..c76b3f896 100644
--- a/skfem/assembly/__init__.py
+++ b/skfem/assembly/__init__.py
@@ -52,7 +52,7 @@
 from itertools import product
 
 from .basis import (Basis, CellBasis, FacetBasis, BoundaryFacetBasis,
-                    InteriorFacetBasis, MortarFacetBasis)
+                    InteriorFacetBasis)
 from .basis import InteriorBasis, ExteriorFacetBasis  # backwards compatibility
 from .dofs import Dofs, DofsView
 from .form import Form, TrilinearForm, BilinearForm, LinearForm, Functional
@@ -94,7 +94,6 @@ def asm(form: Form,
     "FacetBasis",
     "BoundaryFacetBasis",
     "InteriorFacetBasis",
-    "MortarFacetBasis",
     "Dofs",
     "DofsView",
     "TrilinearForm",
diff --git a/skfem/assembly/basis/__init__.py b/skfem/assembly/basis/__init__.py
index 7e6050d3a..f4ebdc081 100644
--- a/skfem/assembly/basis/__init__.py
+++ b/skfem/assembly/basis/__init__.py
@@ -2,7 +2,6 @@
 from .cell_basis import CellBasis  # noqa
 from .facet_basis import FacetBasis  # noqa
 from .interior_facet_basis import InteriorFacetBasis  # noqa
-from .mortar_facet_basis import MortarFacetBasis  # noqa
 
 # aliases
 Basis = CellBasis
diff --git a/skfem/assembly/basis/mortar_facet_basis.py b/skfem/assembly/basis/mortar_facet_basis.py
deleted file mode 100644
index 1e431a237..000000000
--- a/skfem/assembly/basis/mortar_facet_basis.py
+++ /dev/null
@@ -1,57 +0,0 @@
-from typing import Optional, Tuple, Any
-
-from numpy import ndarray
-from skfem.element import Element
-from skfem.mapping import MappingMortar
-from skfem.mesh import Mesh
-from skfem.generic_utils import deprecated
-
-from .facet_basis import FacetBasis
-from ..dofs import Dofs
-
-
-class MortarFacetBasis(FacetBasis):
-
-    @deprecated("FacetBasis + skfem.supermeshing")
-    def __init__(self,
-                 mesh: Mesh,
-                 elem: Element,
-                 mapping: MappingMortar,
-                 intorder: Optional[int] = None,
-                 quadrature: Optional[Tuple[ndarray, ndarray]] = None,
-                 facets: Optional[Any] = None,
-                 dofs: Optional[Dofs] = None,
-                 side: int = 0):
-        """Precomputed global basis on the mortar mesh."""
-
-        if side not in (0, 1):
-            raise Exception("'side' must be 0 or 1.")
-
-        if facets is None:
-            from copy import deepcopy
-            mapping = deepcopy(mapping)
-            mapping.side = side
-            facets = mapping.helper_to_orig[side]
-
-        facets = mesh.normalize_facets(facets)
-
-        super(MortarFacetBasis, self).__init__(
-            mesh,
-            elem,
-            mapping=mapping,
-            intorder=intorder,
-            quadrature=quadrature,
-            facets=facets,
-            dofs=dofs,
-        )
-
-    def default_parameters(self):
-        """Return default parameters for `~skfem.assembly.asm`."""
-        return {
-            **super(MortarFacetBasis, self).default_parameters(),
-            # TODO following is not valid in 3D
-            'h1': self.mapping.maps[0].detDG(self.X,
-                                             self.mapping.helper_to_orig[0]),
-            'h2': self.mapping.maps[1].detDG(self.X,
-                                             self.mapping.helper_to_orig[1]),
-        }
diff --git a/skfem/assembly/form/form.py b/skfem/assembly/form/form.py
index abdad7b16..7f45fae64 100644
--- a/skfem/assembly/form/form.py
+++ b/skfem/assembly/form/form.py
@@ -30,7 +30,7 @@ class Form:
     form: Optional[Callable] = None
 
     def __init__(self,
-                 form: Optional[Callable] = None,
+                 form: Optional[Union[Callable, 'Form']] = None,
                  dtype: Union[Type[np.float64],
                               Type[np.complex64]] = np.float64,
                  nthreads: int = 0,
diff --git a/skfem/experimental/__init__.py b/skfem/experimental/__init__.py
index 62f0b04fe..3245d8674 100644
--- a/skfem/experimental/__init__.py
+++ b/skfem/experimental/__init__.py
@@ -2,5 +2,6 @@
 
 
 logger = logging.getLogger(__name__)
-logger.warning("You are using experimental features under skfem.experimental "
-               "that might change at any time without prior notice.")
+logger.warning("Warning: You are using experimental features "
+               "under skfem.experimental that might change at "
+               "any time without prior notice.")
diff --git a/skfem/experimental/supermeshing/__init__.py b/skfem/experimental/supermeshing/__init__.py
index 04ab5f487..8599d07cf 100644
--- a/skfem/experimental/supermeshing/__init__.py
+++ b/skfem/experimental/supermeshing/__init__.py
@@ -1,80 +1,8 @@
-import numpy as np
+# flake8: noqa
+import logging
+from skfem.supermeshing import *
 
-from skfem.mesh import MeshLine, MeshTri
-from skfem.quadrature import get_quadrature
 
-
-def intersect(m1, m2):
-    p1, t1 = m1
-    p2, t2 = m2
-    if p1.shape[0] == 1 and p2.shape[0] == 1:
-        return _intersect1d(p1, t1, p2, t2)
-    elif p1.shape[0] == 2 and p2.shape[0] == 2:
-        return _intersect2d(p1, t1, p2, t2)
-    raise NotImplementedError("The given mesh types not supported.")
-
-
-def _intersect2d(p1, t1, p2, t2):
-    """Two-dimensional supermesh using shapely and bruteforce."""
-    try:
-        from shapely.geometry import Polygon
-        from shapely.strtree import STRtree
-        from shapely.ops import triangulate
-    except Exception:
-        raise Exception("2D supermeshing requires the package 'shapely'.")
-    t = np.empty((3, 0))
-    p = np.empty((2, 0))
-    polys = [Polygon(p1[:, t1[:, itr]].T) for itr in range(t1.shape[1])]
-    ixmap = {id(polys[itr]): itr for itr in range(t1.shape[1])}
-    s = STRtree(polys)
-    ix1, ix2 = [], []
-    for jtr in range(t2.shape[1]):
-        poly1 = Polygon(p2[:, t2[:, jtr]].T)
-        result = s.query(Polygon(p2[:, t2[:, jtr]].T))
-        if len(result) == 0:
-            continue
-        for poly2 in result:
-            tris = triangulate(poly1.intersection(poly2))
-            for tri in tris:
-                p = np.hstack((p, np.vstack(tri.exterior.xy)[:, :-1]))
-                diff = np.max(t) + 1 if t.shape[1] > 0 else 0
-                t = np.hstack((t, np.array([[0], [1], [2]]) + diff))
-                ix1.append(ixmap[id(poly2)])
-                ix2.append(jtr)
-    return (
-        MeshTri(p, t),
-        np.array(ix1, dtype=np.int64),
-        np.array(ix2, dtype=np.int64),
-    )
-
-
-def _intersect1d(p1, t1, p2, t2):
-    """One-dimensional supermesh."""
-    # find unique supermesh facets by combining nodes from both sides
-    p = np.concatenate((p1.flatten().round(decimals=10),
-                        p2.flatten().round(decimals=10)))
-    p = np.unique(p)
-    t = np.array([np.arange(len(p) - 1), np.arange(1, len(p))])
-    p = np.array([p])
-
-    supermap = MeshLine(p, t)._mapping()
-    mps = supermap.F(np.array([[.5]]))
-    ix1 = MeshLine(p1, t1).element_finder()(mps[0, :, 0])
-    ix2 = MeshLine(p2, t2).element_finder()(mps[0, :, 0])
-
-    return MeshLine(p, t), ix1, ix2
-
-
-def elementwise_quadrature(mesh, supermesh=None, tind=None, order=None):
-    """For creating element-by-element quadrature rules."""
-    if order is None:
-        order = 4
-    if supermesh is None:
-        supermesh = mesh
-    X, W = get_quadrature(supermesh.elem, order)
-    mmap = mesh.mapping()
-    smap = supermesh.mapping()
-    return (
-        mmap.invF(smap.F(X), tind=tind),
-        np.abs(smap.detDF(X) / mmap.detDF(X, tind=tind)) * W,
-    )
+logger = logging.getLogger(__name__)
+logger.warning("Warning: skfem.experimental.supermeshing "
+               "has been moved to skfem.supermeshing.")
diff --git a/skfem/mapping/__init__.py b/skfem/mapping/__init__.py
index 99ce2b203..e2e794523 100644
--- a/skfem/mapping/__init__.py
+++ b/skfem/mapping/__init__.py
@@ -8,4 +8,3 @@
 from .mapping import Mapping  # noqa
 from .mapping_affine import MappingAffine  # noqa
 from .mapping_isoparametric import MappingIsoparametric  # noqa
-from .mapping_mortar import MappingMortar  # noqa
diff --git a/skfem/mapping/mapping_mortar.py b/skfem/mapping/mapping_mortar.py
deleted file mode 100644
index 520538f27..000000000
--- a/skfem/mapping/mapping_mortar.py
+++ /dev/null
@@ -1,206 +0,0 @@
-from dataclasses import replace
-from typing import Optional
-
-import numpy as np
-from numpy import ndarray
-
-from ..generic_utils import deprecated
-from ..mesh import Mesh2D
-from .mapping import Mapping
-
-
-class MappingMortar(Mapping):
-    """A mapping from d-1 dimensional reference element to mortar facets."""
-
-    side = 0
-
-    @deprecated("skfem.supermeshing")
-    def __init__(self,
-                 maps,
-                 helper_to_orig,
-                 helper_maps,
-                 super_to_helper,
-                 supermap,
-                 normals):
-        """Should not be called directly."""
-        self.maps = maps
-        self.helper_to_orig = helper_to_orig
-        self.helper_maps = helper_maps
-        self.super_to_helper = super_to_helper
-        self.supermap = supermap
-        self._normals = normals
-
-    @classmethod
-    def init_2D(cls,
-                mesh1: Mesh2D,
-                mesh2: Mesh2D,
-                boundary1: ndarray,
-                boundary2: ndarray,
-                tangent: ndarray):
-        """Create mortar mappings for two 2D meshes via projection.
-
-        Parameters
-        ----------
-        mesh1
-            An object of the type :class:`~skfem.mesh.mesh_2d.Mesh2D`.
-        mesh2
-            An object of the type :class:`~skfem.mesh.mesh_2d.Mesh2D`.
-        boundary1
-            A subset of facets to use from mesh1.
-        boundary2
-            A subset of facets to use from mesh2.
-        tangent
-            A tangent vector defining the direction of the projection.
-
-        """
-        from ..mesh import MeshLine
-        tangent /= np.linalg.norm(tangent)
-
-        # find unique nodes on the two boundaries
-        p1_ix = np.unique(mesh1.facets[:, boundary1].flatten())
-        p2_ix = np.unique(mesh2.facets[:, boundary2].flatten())
-        p1 = mesh1.p[:, p1_ix]
-        p2 = mesh2.p[:, p2_ix]
-
-        def proj(p):
-            """Project onto the line defined by 'tangent'."""
-            return np.outer(tangent, tangent) @ p
-
-        def param(p):
-            """Calculate signed distances of projected points from origin."""
-            y = proj(p)
-            return np.linalg.norm(y, axis=0) * np.sign(np.dot(tangent, y))
-
-        # find unique supermesh facets by combining nodes from both sides
-        param_p1 = param(p1)
-        param_p2 = param(p2)
-        _, ix = np.unique(np.concatenate((param_p1.round(decimals=10),
-                                          param_p2.round(decimals=10))),
-                          return_index=True)
-        ixorig = np.concatenate((p1_ix, p2_ix + mesh1.p.shape[1]))[ix]
-        p = np.array([np.hstack((param(mesh1.p), param(mesh2.p)))])
-        t = np.array([ixorig[:-1], ixorig[1:]])
-
-        # create 1-dimensional supermesh from the intersections of the
-        # projected facet elements
-        p = p[:, np.concatenate((t[0], np.array([t[1, -1]])))]
-        range_max = np.min([np.max(param_p1), np.max(param_p2)])
-        range_min = np.max([np.min(param_p1), np.min(param_p2)])
-        p = np.array([p[0, (p[0] <= range_max) * (p[0] >= range_min)]])
-        t = np.array([np.arange(p.shape[1] - 1), np.arange(1, p.shape[1])])
-        m_super = MeshLine(p, t)
-
-        # helper meshes for creating the mappings
-        m1 = MeshLine(np.sort(param_p1), np.array([np.arange(p1.shape[1] - 1),
-                                                   np.arange(1, p1.shape[1])]))
-        m2 = MeshLine(np.sort(param_p2), np.array([np.arange(p2.shape[1] - 1),
-                                                   np.arange(1, p2.shape[1])]))
-
-        # construct normals by rotating 'tangent'
-        normal = np.array([tangent[1], -tangent[0]])
-        normals = normal[:, None].repeat(t.shape[1], axis=1)
-
-        # initialize mappings (for orienting)
-        map_super = m_super._mapping()
-        map_m1 = m1._mapping()
-        map_m2 = m2._mapping()
-        map_mesh1 = mesh1._mapping()
-        map_mesh2 = mesh2._mapping()
-
-        # matching of elements in the supermesh and the helper meshes
-        mps = map_super.F(np.array([[.5]]))
-        ix1 = np.digitize(mps[0, :, 0], m1.p[0]) - 1
-        ix2 = np.digitize(mps[0, :, 0], m2.p[0]) - 1
-
-        # for each element, map two points to global coordinates, reparametrize
-        # the points, and flip corresponding helper mesh element indices if
-        # sorting is wrong
-        f1mps = .5 * (mesh1.p[:, mesh1.facets[0, boundary1]] +
-                      mesh1.p[:, mesh1.facets[1, boundary1]])
-        sort_boundary1 = np.argsort(param(f1mps))
-        z1 = map_mesh1.G(map_m1.invF(map_super.F(np.array([[.25, .75]])),
-                                     tind=ix1),
-                         find=boundary1[sort_boundary1][ix1])
-        ix1_flip = np.unique(ix1[param(z1[:, :, 1]) < param(z1[:, :, 0])])
-        m1t = m1.t.copy()
-        m1t[:, ix1_flip] = np.flipud(m1t[:, ix1_flip])
-        m1 = replace(m1, t=m1t)
-
-        f2mps = .5 * (mesh2.p[:, mesh2.facets[0, boundary2]] +
-                      mesh2.p[:, mesh2.facets[1, boundary2]])
-        sort_boundary2 = np.argsort(param(f2mps))
-        z2 = map_mesh2.G(map_m2.invF(map_super.F(np.array([[.25, .75]])),
-                                     tind=ix2),
-                         find=boundary2[sort_boundary2][ix2])
-        ix2_flip = np.unique(ix2[param(z2[:, :, 1]) < param(z2[:, :, 0])])
-        m2t = m2.t.copy()
-        m2t[:, ix2_flip] = np.flipud(m2t[:, ix2_flip])
-        m2 = replace(m2, t=m2t)
-
-        # construct normals by rotating 'tangent'
-        normal = np.array([tangent[1], -tangent[0]])
-        normals = normal[:, None].repeat(t.shape[1], axis=1)
-
-        # initialize mappings (for orienting)
-        map_super = m_super._mapping()
-        map_m1 = m1._mapping()
-        map_m2 = m2._mapping()
-        map_mesh1 = mesh1._mapping()
-        map_mesh2 = mesh2._mapping()
-
-        # matching of elements in the supermesh and the helper meshes
-        mps = map_super.F(np.array([[.5]]))
-        ix1 = np.digitize(mps[0, :, 0], m1.p[0]) - 1
-        ix2 = np.digitize(mps[0, :, 0], m2.p[0]) - 1
-
-        return cls((map_mesh1, map_mesh2),
-                   (boundary1[sort_boundary1][ix1],
-                    boundary2[sort_boundary2][ix2]),
-                   (map_m1, map_m2),
-                   (ix1, ix2),
-                   map_super,
-                   normals)
-
-    def F(self, X, tind=None):
-        return self.maps[self.side].F(X, tind=tind)
-
-    def invF(self, x, tind=None):
-        return self.maps[self.side].invF(x, tind=tind)
-
-    def detDF(self, X, tind=None):
-        return self.maps[self.side].detDF(X, tind=tind)
-
-    def DF(self, X, tind=None):
-        return self.maps[self.side].DF(X, tind=tind)
-
-    def invDF(self, X, tind=None):
-        return self.maps[self.side].invDF(X, tind=tind)
-
-    def G(self,
-          X: ndarray,
-          find: Optional[ndarray] = None,
-          side: Optional[int] = None) -> ndarray:
-        side = self.side if side is None else side
-        return self.maps[side].G(
-            self.helper_maps[side].invF(self.supermap.F(X),
-                                        tind=self.super_to_helper[side]),
-            find=self.helper_to_orig[side]
-        )
-
-    def detDG(self,
-              X: ndarray,
-              find: Optional[ndarray] = None) -> ndarray:
-        if X.shape[0] == 1:
-            segs = self.G(np.array([[0., 1.]]), side=0)
-            return np.sqrt(np.diff(segs[0], axis=1) ** 2 +
-                           np.diff(segs[1], axis=1) ** 2)
-        raise NotImplementedError
-
-    def normals(self,
-                X: ndarray,
-                tind: ndarray,
-                find: ndarray,
-                t2f: ndarray) -> ndarray:
-        return np.repeat(self._normals[:, :, None],
-                         X.shape[-1],
-                         axis=2)
diff --git a/skfem/supermeshing.py b/skfem/supermeshing.py
new file mode 100644
index 000000000..322aa6195
--- /dev/null
+++ b/skfem/supermeshing.py
@@ -0,0 +1,114 @@
+import numpy as np
+
+from skfem.mesh import MeshLine, MeshTri
+from skfem.quadrature import get_quadrature
+
+
+def intersect(m1, m2):
+    """Create a supermesh between two nonmatching (1D or 2D) meshes.
+
+    Parameters
+    ----------
+    m1
+        The first mesh.
+    m2
+        The second mesh.
+
+    Returns
+    -------
+    supermesh
+        Mesh object with line (1D) or triangle (2D) elements.
+    ix1
+        Index of each supermesh element within the first mesh.
+    ix2
+        Index of each supermesh element within the second mesh.
+
+    """
+    p1, t1 = m1
+    p2, t2 = m2
+    if p1.shape[0] == 1 and p2.shape[0] == 1:
+        return _intersect1d(p1, t1, p2, t2)
+    elif p1.shape[0] == 2 and p2.shape[0] == 2:
+        return _intersect2d(p1, t1, p2, t2)
+    raise NotImplementedError("The given mesh types not supported.")
+
+
+def _intersect2d(p1, t1, p2, t2):
+    """Two-dimensional supermesh using shapely and bruteforce."""
+    try:
+        from shapely.geometry import Polygon
+        from shapely.strtree import STRtree
+        from shapely.ops import triangulate
+    except Exception:
+        raise Exception("2D supermeshing requires the package 'shapely>=2'.")
+    t = np.empty((3, 0))
+    p = np.empty((2, 0))
+    polys = [Polygon(p1[:, t1[:, itr]].T) for itr in range(t1.shape[1])]
+    ixmap = {id(polys[itr]): itr for itr in range(t1.shape[1])}
+    s = STRtree(polys)
+    ix1, ix2 = [], []
+    for jtr in range(t2.shape[1]):
+        poly1 = Polygon(p2[:, t2[:, jtr]].T)
+        result = s.query(Polygon(p2[:, t2[:, jtr]].T))
+        if len(result) == 0:
+            continue
+        for poly2 in s.geometries.take(result):
+            tris = triangulate(poly1.intersection(poly2))
+            for tri in tris:
+                p = np.hstack((p, np.vstack(tri.exterior.xy)[:, :-1]))
+                diff = np.max(t) + 1 if t.shape[1] > 0 else 0
+                t = np.hstack((t, np.array([[0], [1], [2]]) + diff))
+                ix1.append(ixmap[id(poly2)])
+                ix2.append(jtr)
+    return (
+        MeshTri(p, t),
+        np.array(ix1, dtype=np.int64),
+        np.array(ix2, dtype=np.int64),
+    )
+
+
+def _intersect1d(p1, t1, p2, t2):
+    """One-dimensional supermesh."""
+    # find unique supermesh facets by combining nodes from both sides
+    p = np.concatenate((p1.flatten().round(decimals=10),
+                        p2.flatten().round(decimals=10)))
+    p = np.unique(p)
+    t = np.array([np.arange(len(p) - 1), np.arange(1, len(p))])
+    p = np.array([p])
+
+    supermap = MeshLine(p, t)._mapping()
+    mps = supermap.F(np.array([[.5]]))
+    ix1 = MeshLine(p1, t1).element_finder()(mps[0, :, 0])
+    ix2 = MeshLine(p2, t2).element_finder()(mps[0, :, 0])
+
+    return MeshLine(p, t), ix1, ix2
+
+
+def elementwise_quadrature(mesh, supermesh=None, tind=None, intorder=None):
+    """For creating element-by-element quadrature rules.
+
+    Parameters
+    ----------
+    mesh
+        The mesh for which to create the quadrature rules.
+    supermesh
+        Created using skfem.supermeshing.intersect.
+    tind
+        A subset of elements
+    intorder
+        Integration order, by default equal to 4.
+
+    """
+    if intorder is None:
+        intorder = 4
+    if supermesh is None:
+        raise Exception("elementwise_quadrature: User must provide "
+                        "'supermesh' keyword argument which has been "
+                        "created using skfem.supermeshing.intersect.")
+    X, W = get_quadrature(supermesh.elem, intorder)
+    mmap = mesh.mapping()
+    smap = supermesh.mapping()
+    return (
+        mmap.invF(smap.F(X), tind=tind),
+        np.abs(smap.detDF(X) / mmap.detDF(X, tind=tind)) * W,
+    )
diff --git a/tests/test_basis.py b/tests/test_basis.py
index c932730f5..06e1422d0 100644
--- a/tests/test_basis.py
+++ b/tests/test_basis.py
@@ -10,9 +10,8 @@
 from skfem import BilinearForm, LinearForm, asm, solve, condense, projection
 from skfem.mesh import (Mesh, MeshTri, MeshTet, MeshHex,
                         MeshQuad, MeshLine1, MeshWedge1)
-from skfem.assembly import (CellBasis, FacetBasis, Dofs, Functional,
-                            MortarFacetBasis)
-from skfem.mapping import MappingIsoparametric, MappingMortar
+from skfem.assembly import CellBasis, FacetBasis, Dofs, Functional
+from skfem.mapping import MappingIsoparametric
 from skfem.element import (ElementVectorH1, ElementTriP2, ElementTriP1,
                            ElementTetP2, ElementHexS2, ElementHex2,
                            ElementQuad2, ElementLineP2, ElementTriP0,
@@ -360,78 +359,6 @@ def test_pickling():
     )
 
 
-@pytest.mark.parametrize(
-    "m1, m2, lenright",
-    [
-        (
-            MeshTri().refined(3),
-            MeshTri().translated((1., 0.)).refined(2),
-            1.,
-        ),
-        (
-            MeshTri.init_refdom().refined(3).with_boundaries(
-                {'right': lambda x: x[0] + x[1] == 1}
-            ),
-            MeshTri().translated((1., 0.)).refined(2),
-            np.sqrt(2),
-        )
-    ]
-)
-def test_mortar_basis(m1, m2, lenright):
-    # some sanity checks for MortarBasis
-    e = ElementTriP1()
-
-    mort = MappingMortar.init_2D(m1,
-                                 m2,
-                                 m1.boundaries['right'],
-                                 m2.boundaries['left'],
-                                 [0., 1.])
-
-    mb = [
-        MortarFacetBasis(m1, e, mapping=mort, intorder=4, side=0),
-        MortarFacetBasis(m2, e, mapping=mort, intorder=4, side=1),
-    ]
-
-    @Functional
-    def unity(w):
-        return 1.
-
-    @LinearForm
-    def load(v, w):
-        return 1. * v
-
-    @BilinearForm
-    def mass(u, v, w):
-        return u * v
-
-    assert_allclose(mb[0].default_parameters()['h1'],
-                    mb[1].default_parameters()['h1'])
-
-    assert_allclose(mb[0].default_parameters()['h2'],
-                    mb[1].default_parameters()['h2'])
-
-    assert (mb[0].default_parameters()['h1']
-            < mb[0].default_parameters()['h2']).all()
-
-    assert_almost_equal(unity.assemble(mb[0]), lenright)
-    assert_almost_equal(unity.assemble(mb[1]), lenright)
-
-    assert_almost_equal(load.assemble(mb[0]).dot(m1.p[1]), .5 * lenright)
-    assert_almost_equal(load.assemble(mb[1]).dot(m2.p[1]), .5 * lenright)
-
-    # integral is over the domain of the first argument
-    assert_almost_equal((mass.assemble(mb[0], mb[1])
-                         .dot(m1.p[0] * 0. + 1.)
-                         .dot(m2.p[0] * 0. + 1.)), lenright)
-
-    assert_almost_equal((mass.assemble(mb[1], mb[0])
-                         .dot(m2.p[0] * 0. + 1.)
-                         .dot(m1.p[0] * 0. + 1.)), lenright)
-
-    assert_allclose(mass.assemble(mb[0], mb[1]).toarray(),
-                    mass.assemble(mb[1], mb[0]).T.toarray())
-
-
 @pytest.mark.parametrize(
     "m, e, fun",
     [
diff --git a/tests/test_examples.py b/tests/test_examples.py
index 02a178f82..eb45b0b0a 100644
--- a/tests/test_examples.py
+++ b/tests/test_examples.py
@@ -396,3 +396,10 @@ class TestEx48(TestCase):
     def runTest(self):
         import docs.examples.ex48 as ex48
         self.assertAlmostEqual(np.max(ex48.u), 0.0012653264834919884)
+
+
+class TestEx49(TestCase):
+
+    def runTest(self):
+        import docs.examples.ex49 as ex49
+        self.assertLess(np.abs(ex49.l2), 0.0012)
diff --git a/tests/test_mapping.py b/tests/test_mapping.py
index e37afa377..a22de28c7 100644
--- a/tests/test_mapping.py
+++ b/tests/test_mapping.py
@@ -4,7 +4,6 @@
 from skfem.mesh import MeshHex, MeshQuad, MeshTri
 from skfem.element import ElementHex1, ElementQuad1, ElementHex2
 from skfem.assembly import FacetBasis
-from skfem.mapping import MappingMortar
 
 
 class TestIsoparamNormals(unittest.TestCase):
@@ -88,54 +87,3 @@ class TestInverseMappingHex2(TestInverseMappingHex):
     """This should be equivalent to TestInverseMappingHex."""
 
     element = ElementHex2
-
-
-class TestMortarPair(unittest.TestCase):
-    """Check that mapped points match."""
-
-    mesh1_type = MeshTri
-    mesh2_type = MeshTri
-    nrefs1 = 2
-    nrefs2 = 3
-    translate_y = 0.0
-
-    def init_meshes(self):
-        m1 = self.mesh1_type().refined(self.nrefs1)
-        m2 = self.mesh2_type().refined(self.nrefs2).translated((1.0, self.translate_y))
-        return m1, m2
-
-    def runTest(self):
-        m1, m2 = self.init_meshes()
-        mp = MappingMortar.init_2D(m1, m2,
-                                   m1.facets_satisfying(lambda x: x[0] == 1.),
-                                   m2.facets_satisfying(lambda x: x[0] == 1.),
-                                   np.array([0., 1.]))
-        test_points = np.array([np.linspace(0., 1., 7)])
-        self.assertTrue((mp.G(test_points) -
-                         mp.G(test_points) < 1e-10).all())
-
-
-class TestMortarPairTriQuad(TestMortarPair):
-
-    mesh1_type = MeshTri
-    mesh2_type = MeshQuad
-
-
-class TestMortarPairQuadQuad(TestMortarPair):
-
-    mesh1_type = MeshQuad
-    mesh2_type = MeshQuad
-
-
-class TestMortarPairNoMatch1(TestMortarPair):
-
-    mesh1_type = MeshQuad
-    mesh2_type = MeshTri
-    translate_y = 0.1
-
-
-class TestMortarPairNoMatch2(TestMortarPair):
-
-    mesh1_type = MeshQuad
-    mesh2_type = MeshTri
-    translate_y = -np.pi / 10.