Modeling a fixed joint

Hello everyone,

I would like to model a fixed joint between two parts and I use the 3D box with holes example as inspiration.

So, this part and a second one are modeled and I would like the joint to be modeled through the holes as if the two pieces were linked through bolts or screws. I fix one face of a part and apply a force on a face of the other one and the aim is to transfer loads through linked holes (and not to fix them all).

Here is my code missing the link part which I hope you help me with:

from netgen.occ import *
from ngsolve import *
import netgen.gui

box1 = Box((0,0,0), (3,0.6,1))
cyl1 = sum( [Cylinder((0.5+i,0,0.5), Y, 0.25,0.8) for i in range(3)] )
cyl1.faces.name="cyl1"
geo1 = box1-cyl1

box2 = Box((0,0.6,0), (3,1.2,1))
cyl2 = sum( [Cylinder((0.5+i,0.6,0.5), Y, 0.25,0.8) for i in range(3)] )
cyl2.faces.name="cyl2"
geo2 = box2-cyl2

geo = geo1 + geo2

geo1.faces.Min(X).name = "fix"
geo2.faces.Max(X).name = "force"

mesh = geo.GenerateMesh(maxh=0.1).Curve(3)

E, nu = 210, 0.2
mu  = E / 2 / (1+nu)
lam = E * nu / ((1+nu)*(1-2*nu))

def Stress(strain):
    return 2*mu*strain + lam*Trace(strain)*Id(3)
fes = VectorH1(mesh, order=3, dirichlet="fix")
u,v = fes.TnT()
gfu = GridFunction(fes)

with TaskManager():
    a = BilinearForm(InnerProduct(Stress(Sym(Grad(u))), Sym(Grad(v))).Compile()*dx)
    pre = preconditioners.BDDC(a)
    a.Assemble()

force = CF((1e-3,0,0))
f = LinearForm(force*v*ds("force")).Assemble()
inv = solvers.CGSolver(a.mat, pre, plotrates=True, tol=1e-8)
gfu.vec.data = inv * f.vec

with TaskManager():
    fesstress = MatrixValued(H1(mesh,order=3), symmetric=True)
    gfstress = GridFunction(fesstress)
    gfstress.Interpolate (Stress(Sym(Grad(gfu))))
Draw(gfu, mesh, deformation=True)

Thank you in advance.
Sabri