Add supermeshing, remove MortarFacetBasis and MortarMapping

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

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

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))

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()

requirements.txt

@@ -10,3 +10,4 @@ flake8
 sphinx~=6.1
 autograd
 pep517
+shapely

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',

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",

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

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,

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.")

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.")

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