Implementation

import numpy as np
import scipy
from scipy.sparse import bmat, csr_matrix, coo_matrix, lil_matrix

import skfem
from skfem.helpers import dot, ddot, sym_grad
from skfem import InteriorBasis, FacetBasis, ElementTetP1, ElementVectorH1

from skfem.supermeshing import intersect, elementwise_quadrature


mesh_top = skfem.MeshTet.init_tensor(
    np.linspace(0, 2, 20),
    np.linspace(0, 2, 20),
    np.linspace(0, 1, 10)
)
mesh_bot = skfem.MeshTet.init_tensor(
    np.linspace(0.5, 1.5, 20),
    np.linspace(0.5, 1.5, 20),
    np.linspace(-0.5, 0.0, 5)
)


def is_contact_surface(x):
    return np.isclose(x[2], 0.0)


def is_force_surface(x):
    return np.isclose(x[2], 1.0)


def is_dirichlet_surface(x):
    return np.isclose(x[2], -0.5)


mesh_top = mesh_top.with_boundaries({
    "contact": is_contact_surface,
    "force": is_force_surface
})
mesh_bot = mesh_bot.with_boundaries({
    "contact": is_contact_surface,
    "dirichlet": is_dirichlet_surface
})


elem = skfem.ElementTetP1()
u_elem = skfem.ElementVectorH1(elem)
basis_top = skfem.InteriorBasis(mesh_top, u_elem)
basis_bot = skfem.InteriorBasis(mesh_bot, u_elem)

E, nu = 210e9, 0.3
mu = E/(2*(1+nu))
lam = E*nu/((1+nu)*(1-2*nu))
@skfem.BilinearForm
def a_elastic(u, v, w):
    eps_u = sym_grad(u)      # symmetric gradient of u
    eps_v = sym_grad(v)
    return lam * np.trace(eps_u) * np.trace(eps_v) + 2*mu * ddot(eps_u, eps_v)


K1 = a_elastic.assemble(basis_top)
K2 = a_elastic.assemble(basis_bot)

interface_facets = mesh_bot.facets_satisfying(lambda pts: np.isclose(pts[2], 0.0))
mortar_basis = FacetBasis(mesh_bot, elem, facets=interface_facets) 


basis_top = skfem.Basis(mesh_top, ElementVectorH1(skfem.ElementTetP1()))
basis_bot = skfem.Basis(mesh_bot, ElementVectorH1(skfem.ElementTetP1()))

facets_if = mesh_bot.facets_satisfying(is_contact_surface)
lambda_basis = skfem.FacetBasis(
    mesh_bot,
    skfem.ElementTetP1(),
    facets=facets_if
)


elem_s = ElementTetP1()
elem_v = ElementVectorH1(elem_s)
m1t, orig1 = mesh_top.trace('contact', mtype=skfem.MeshTri, project=lambda p: p[[0, 1]])
m2t, orig2 = mesh_bot.trace('contact', mtype=skfem.MeshTri, project=lambda p: p[[0, 1]])
m12, t1, t2 = intersect(m1t, m2t)
quad1 = elementwise_quadrature(m1t, m12, t1)
quad2 = elementwise_quadrature(m2t, m12, t2)

fb_u_top = FacetBasis(
    mesh_top, elem_v, facets=orig1[t1], quadrature=quad1
)
fb_u_bot = FacetBasis(
    mesh_bot, elem_v, facets=orig2[t2], quadrature=quad2
)
fb_lam_v = FacetBasis(mesh_bot, elem_v, facets=orig2[t2], quadrature=quad2)


@skfem.BilinearForm
def b_top(u1, u2, v1, v2, w):
    if v1.ndim == 3:
        return dot(v1, u2)
    return 0.0 * (w.n[0] * 0 + 1)

@skfem.BilinearForm
def b_bot(u1, u2, v1, v2, w):
    if v1.ndim == 3:
        return - dot(v1, u2)
    return 0.0 * (w.n[0] * 0 + 1)


Bfull_top = skfem.asm(b_top, fb_u_top*fb_lam_v)
Bfull_bot = skfem.asm(b_bot, fb_u_bot*fb_lam_v)

Nu_top = int(fb_u_top.N)
Nu_bot = int(fb_u_bot.N)
Nlam_v = int(fb_lam_v.N)

B1 = Bfull_top[:Nu_top, Nu_top:Nu_top + Nlam_v]
B2 = Bfull_bot[:Nu_bot, Nu_bot:Nu_top + Nu_bot]


print("K1:", K1.shape)
print("K2:", K1.shape)
print("B1:", B1.shape)
print("B2:", B1.shape)
print("basis_top.N:", basis_top.N)

B1_full = lil_matrix((B1.shape[0], basis_top.N))
B2_full = lil_matrix((B2.shape[0], basis_bot.N))

K = bmat([
    [K1,        None,     B1],
    [None,      K2,      -B2],
    [B1.T,     -B2.T,     None]
], format='csr')

# Assemble RHS vector
F_facet_ids = mesh_top.facets_satisfying(is_force_surface)
# F_facet_ids = mesh.facets_satisfying(lambda x: np.isclose(x[0], 1.0))

fbasis = skfem.FacetBasis(mesh_top, elem_v, facets=F_facet_ids)
total_force = 100.0  # [N]


@skfem.Functional
def surface_measure(w):
    return 1.0


A_proj_z = surface_measure.assemble(fbasis)
pressure = total_force / A_proj_z


# ---- traction: t = p * n （l(v) = ∫_Γ (p n) · v ds）----
@skfem.LinearForm
def l_comp(v, w):
    return pressure * v[2]

# @skfem.LinearForm
# def l_comp(v, w):
#     return pressure * dot(w.n, v)


top_F = l_comp.assemble(fbasis)
bot_F = np.zeros(mesh_bot.p.shape[1] * mesh_bot.dim())
F = np.hstack([top_F, bot_F, np.zeros(B1.T.shape[0])])


D_facets = mesh_bot.facets_satisfying(is_dirichlet_surface)
D_dofs = basis_bot.get_dofs(facets=D_facets).all() + K1.shape[0]

AII, bI, xI, I = skfem.condense(K, F, D=D_dofs)