# Deforming a mesh using MPM particles?

Hello!
I’m using a triangular mesh to seed a shape in a MPM simulation. Is there a way of moving the mesh together with the particles? (I mostly need the mesh for rendering, but also for computing stuff like its curvature)
Thank you in advance!
roberto

I imagine @yuanming used something like that for generating the armadillo-man images?

1 Like

Hi! I’m not super sure how @yuanming made this cool image at that moment, but if you want to render 3D meshes, hope my extension library Taichi THREE, despite it’s still under-development, could helps

2 Likes

thank you @archibate! will try it

In my case, the starting point is a mesh. I used code from taichi_elements to voxelise the mesh and generate the MPM points. After that, my algorithm deforms the mesh, quite massively . At the end of the simulation, the points had move well far away from their original position. I would like to deform the mesh accordingly.
I was thinking to represent the position of each vertex in the original mesh as a linear combination of the k>=3 nearest MPM points. At the end of the simulation, I would compute the new position of the mesh vertices using the weights of the initial linear combination. Does this sound sensible to you? I’m afraid that computing k nearest neighbours with so many points and vertices may be super slow

Yes, that’s indeed super slow, but you may use a neighbour lookup table or octree for quicker lookup, e.g., pseudo-code:

``````for i in pos:  # scan neighbours
x, y = int(pos[i])
ti.append(neighbour, [x, y], i)
``````

Hello! I found a way of deforming the mesh, and it works like a charm . I thought I would share it here in case it can help someone else.
Here’s my setup: I have a mesh with vertices v and faces f, which I use to seed particles p. The particles get deformed, and turn into q. And what I want is to move the vertices v following the deformation of p into q.
There’s two steps: first, i compute k neighbouring points p for each mesh vertex v (using a kdtree). Second, I use moving least squares to compute a deformation using the k p neighbours of each vertex as control points.
Here’s the complete algorithm:

``````import numpy as np
from sklearn.neighbors import KDTree

# you need to get the following numpy arrays
# v: mesh vertices
# p: initial particle positions
# q: final particle positions
# The new vertices will be in V
k=5
tree = KDTree(p)
dist, ind = tree.query(v, k=k)
alpha = 1
w = np.array(1/dist**alpha)
w[w==np.inf]=2**31-1
V = v.copy()
for index in range(len(v)):
vi=v[index]
pi=p[ind[index]]
qi=q[ind[index]]
wi=w[index]
pcentroid = pi.T.dot(wi)/np.sum(wi)
qcentroid = qi.T.dot(wi)/np.sum(wi)
ptr = pi - pcentroid
qtr = qi - qcentroid
M = ptr.T.dot(np.diag(wi).dot(ptr))
detM = np.linalg.det(M)
if detM < 1e-8:
V[index] = vi - pcentroid + qcentroid
else:
invM = np.linalg.inv(M)
for j in range(5): # k = 5
res = (vi-pcentroid).dot(invM.dot(wi[j]*ptr[j]))
A[index,j] = res
V[index] = np.sum(np.diag(A[index,:]).dot(qi),axis=0) + qcentroid
``````

Finally, here’s pictures of one of my brain meshes before and after deformations.

3 Likes

I realised that most of the times detM<1e-8, and the code was simply doing an average of the k nearest neighbours.
Here’s an improved version: Instead of Moving Least Squares, I compute a polar transformation of the original k-nearest particles to their final positions, plus a scaling factor. I then apply that rotation and scaling to the vertex. Here’s the new code:

``````import numpy as np
from sklearn.neighbors import KDTree
from numpy.linalg import svd

# you need to get the following numpy arrays
# v: mesh vertices
# p: initial particle positions
# q: final particle positions
# The new vertices will be in V

def polar(P, Q):
'''
Input: moving and reference Nx3 arrays of coordinates
Returns 3x3 rotation matrix and linear scaling
'''
# scale
TP=P.T.dot(P)
TQ=Q.T.dot(Q)
s = (np.linalg.det(TQ)/np.linalg.det(TP))**(1/2*1/3)

# rotation
T = Q.T.dot(P)
U,_,V=svd(T)
R = V.T.dot(U.T)
return R, s

k=10
tree = KDTree(p)
dist, ind = tree.query(v, k=k)
alpha = 1
w = np.array(1/dist**alpha)
w[w==np.inf]=2**31-1
V = v.copy()

for index in range(len(v)):
vi=v[index]
pi=p[ind[index]]
qi=q[ind[index]]
wi=w[index]
pcentroid = pi.T.dot(wi)/np.sum(wi)
qcentroid = qi.T.dot(wi)/np.sum(wi)
ptr = pi - pcentroid
qtr = qi - qcentroid
vtr = vi - pcentroid
R, scale = polar(ptr, qtr)
V[index] = scale*vtr.dot(R) + qcentroid
``````