Support composite Basis (e.g., for mortaring)

[kinnala]

Right now the use of "composite forms" such as

def form(u1, u2, v1, v2, w):
    return ...

requires that we create ElementComposite, e.g., multiplying ElementTriP2() * ElementTriP1() and then u1 and u2 would correspond to different unknowns but be defined over the same mesh.

Now there are situations where you would like to have u1 and u2 be defined using different Basis objects instead. One simple example is when using InteriorFacetBasis(..., side=0) and InteriorFacetBasis(..., side=1) to implement DG methods. Another example is when using mortaring (see ex04) so that the FacetBasis objects are coming from the boundary of two different meshes.

This PR includes works toward something like BasisComposite which allows multiplying two compatible Basis objects for the same effect, i.e. using both of them within forms.

[kinnala]

Here is what I've been testing:

from skfem import *
from skfem.experimental.autodiff import *
from skfem.experimental.autodiff.helpers import *
import jax.numpy as np
import numpy as onp
import matplotlib.pyplot as plt

m1 = MeshTri.init_tensor(np.linspace(0, 5, 60), np.linspace(0, 0.25, 5))
m2 = MeshTri.init_tensor(np.linspace(0, 5, 51), np.linspace(-0.5, -0.25, 4))

from skfem.supermeshing import intersect, elementwise_quadrature

m1t, orig1 = m1.trace('bottom', mtype=MeshLine, project=lambda p: np.array(p[0]))
m2t, orig2 = m2.trace('top', mtype=MeshLine, project=lambda p: np.array(p[0]))
m12, t1, t2 = intersect(m1t, m2t)
e1 = ElementVector(ElementTriP1())
e2 = ElementVector(ElementTriP1())

basis1 = Basis(m1, e1)
basis2 = Basis(m2, e2)
fbasis1 = FacetBasis(m1, e1, quadrature=elementwise_quadrature(m1t, m12, t1), facets=orig1[t1])
fbasis2 = FacetBasis(m2, e2, quadrature=elementwise_quadrature(m2t, m12, t2), facets=orig2[t2])

fbasis = fbasis1 * fbasis2

@NonlinearForm(hessian=True)
def g(u1, u2, w):
    eps = 1e-2
    return 1 / eps * np.minimum(0.25 + u1[1] - u2[1], 0) ** 2

@NonlinearForm(hessian=True)
def J1(u, w):
    eps = 1/2 * (grad(u) + transpose(grad(u)) + mul(transpose(grad(u)), grad(u)))
    sig = 2 * eps + eye(trace(eps), 2)
    return 1/2 * ddot(sig, eps)


@NonlinearForm(hessian=True)
def J2(u, w):
    eps = 1/2 * (grad(u) + transpose(grad(u)) + mul(transpose(grad(u)), grad(u)))
    sig = 2 * eps + eye(trace(eps), 2)
    return 1/2 * ddot(sig, eps) - 1e-2 * u[1]

x = fbasis.zeros()
dx = onp.array(fbasis.zeros())

for itr in range(40):
    jac1, rhs1 = J1.assemble(basis1, x=x[:basis1.N])
    jac2, rhs2 = J2.assemble(basis2, x=x[basis1.N:])
    jac = bmat([[jac1, None], [None, jac2]])
    rhs = np.concatenate((rhs1, rhs2))
    jacg, rhsg = g.assemble(fbasis, x=x)
    jac = jac + jacg
    rhs = rhs + rhsg
    dx = solve(*enforce(jac, onp.array(rhs), D=onp.array(onp.concatenate((
        basis1.get_dofs({'left', 'right'}).all(),
        basis2.get_dofs({'left', 'right'}).all() + basis1.N,
    )))))
    x += dx
    print(np.linalg.norm(dx))
    if np.linalg.norm(dx) < 1e-8:
        break
    mdefo1 = m1.translated(x[:basis1.N][basis1.nodal_dofs])
    mdefo2 = m2.translated(x[basis1.N:][basis2.nodal_dofs])
    ax = mdefo1.draw()
    mdefo2.draw(ax=ax)

[kinnala]

Now here is something to consider:

• How should basis.N be picked? (Mortaring vs. DG is different?)
• Make sure mesh types can be different (for mortaring)
• How about projection from one mesh to another?

[kinnala]

Fancy pics from the code above. One of the earlier iterations:

The final iteration: