A working contact detecting method for 2d case

[ZAARAOUI999]

Hello Mr. Andreas, I'm trying to give some ideas that may help you to implement Contact between solidbodies, so I have got some of good results with detecting contact points, and as I know this should be the first step to implement a contact behavior. The work below may need to be enhanced as a for-loop is taking a place there.

The code:

import felupe as fem
import numpy as np

import matplotlib.pyplot as plt
from scipy.spatial import distance_matrix
from matplotlib import animation, rc
rc('animation', html='html5')

def get_bounding_points(*args): # taked from @adtzlr solution in discussion #331
    """Get the points that lie only on the boundary surface of a solidbody"""
    rbs = list(args)
    bps = [rb.mesh.points[np.isin(np.arange(rb.mesh.npoints), rb.mesh.cells_faces)] for rb in rbs]
    return bps

def get_contact_points(rb1, rb2, tol=1):
    """Get the contact points between two regions"""
    
    #getting bounding points
    bp1, bp2 = get_bounding_points(rb1, rb2)
    #calculating the distance between the two boundaries points bp1 and bp2
    #scipy.spatial.distance_matrix: gives the distance between two ndarrays (M,K) (N,K) ==> (M,N)
    d = distance_matrix(bp1, bp2)
    d = np.round(d,4)
    #catching min distances (tolerance: tol)
    min_mask = np.where(d<=tol)
    minimal_d = d[min_mask]

    #removing duplicates
    minimal_d = np.unique(minimal_d)

    #getting neighbors points
    mask = np.any(np.isin(d,minimal_d), axis=1)
    points = d[mask]
    
    #getting contact points for first solid
    index_1 = np.arange(bp2.shape[0]).reshape(-1,1)
    index_1 = np.repeat(index_1, bp1.shape[0], axis=1).T + np.arange(bp1.shape[0]).reshape(-1,1)
    bp1_c = bp1[index_1[mask][:,:1]].reshape(-1,2)
    
    #getting contact points for second solid
    index_ = np.arange(bp2.shape[0])
    index_2 = set()
    for i, point in enumerate(points):
      for i, dmin in enumerate(minimal_d):
           mask = np.where(point==dmin)
           indice = index_[mask]
           if len(indice)>0:
               index_2.add(indice[0])
    bp2_c = bp2[list(index_2)]
    
    return bp1, bp2, bp1_c, bp2_c


from scipy.spatial import ConvexHull as conv
def plot_mesh(ax, points, c='k', name=''):
    """Plot mesh boundary"""
    hull = conv(points)
    x, y = points[hull.vertices, 0], points[hull.vertices,1]
    x = np.append(x, x[0])
    y = np.append(y, y[0])
    ax.plot(x, y, c='k')
    ax.text(x.mean(), y.mean(), name, c=c)


m1 = fem.Rectangle(b=(10,10), n=10)
angle = np.linspace(0,90,90)
m2 = [fem.mesh.rotate(fem.Rectangle(a=(0,10), b=(5,15), n=5), angle_deg=-a,
                      axis=-1, center=(5,10)) for a in angle]

rb1 = fem.RegionQuadBoundary(m1)
rb2 = [fem.RegionQuadBoundary(m) for m in m2]


fig, ax = plt.subplots()

def animate(i):
    bp1, bp2, bp1_c, bp2_c = get_contact_points(rb1, rb2[i], .6)
    ax.clear()
    ax.set_title('Contact Points')
    ax.set_xlabel(r'x$\longrightarrow$')
    ax.set_ylabel(r'y$\longrightarrow$')
    sc1 = ax.scatter(bp1_c[...,0], bp1_c[...,1], s=30, marker='o', c='b',
                label=r'Contact points $(S_1)$')
    sc2 = ax.scatter(bp2_c[...,0], bp2_c[...,1], s=80, marker='o', facecolors='none',
                edgecolors='r', label=r'Contact points $(S_2)$')

    plot_mesh(ax, bp1, c= 'b', name= r'$S_1$')
    plot_mesh(ax, bp2, c='r', name= r'$S_2$')

    return sc1, sc2
    

anim = animation.FuncAnimation(fig, animate, frames=90, interval=90, blit=True)
anim.save('contact.gif')

The result:

I think that this wasn't a waste of time, right?

[adtzlr]

Looks great, don't worry, nothing is waste of time if you learn something and it's fun! Unfortunately contact algos are a field of which I do not have much of theoretical expertice.

[ZAARAOUI999]

It's my case too, I just taked a look at some of published articles that talked about the solid contact formulations. This is should be a new challenge and it will be an amazing thing. Thanks for the positive vibes.