armadillo.txt

numPoints 117
0.365071 0.0626429
0.489804 0.0447059
0.595851 0.0482979
0.705 0.0581818
1.20242 0.0495
1.3474 0.0344667
1.49 0.0335058
1.63627 0.0534804
0.403448 0.148879
0.491482 0.220926
0.725119 0.232738
0.738461 0.138333
1.17868 0.142075
1.18687 0.226161
1.44703 0.221139
1.55032 0.149435
0.466786 0.356875
0.605415 0.382114
0.750454 0.329091
1.17 0.314375
1.3172 0.331585
1.45265 0.349167
0.453529 0.494779
0.726864 0.529915
0.761097 0.43378
0.925569 0.560633
1.1625 0.426413
1.1875 0.523636
1.47825 0.487302
0.451154 0.644423
0.7529 0.6351
0.859196 0.646071
0.967213 0.675082
1.10859 0.655128
1.18045 0.591136
1.48483 0.639489
0.454333 0.7975
0.62754 0.753659
0.820471 0.72153
0.983883 0.776861
1.04658 0.718816
1.26177 0.746573
1.44942 0.770814
0.518971 0.925
0.611719 1.01461
0.753619 0.926785
0.926779 0.925082
1.11516 0.900788
1.28965 0.970139
1.38754 0.878261
0.66625 1.13742
0.85003 1.0886
1.03465 1.09401
1.19906 1.17576
1.21595 1.05182
0.53575 1.35017
0.635737 1.28016
0.783853 1.24314
0.928477 1.26568
1.08683 1.3065
1.25554 1.2723
1.42645 1.35261
0.216884 1.51594
0.327563 1.44937
0.416615 1.37992
0.635131 1.48694
0.792381 1.41122
0.962164 1.44965
1.17129 1.47516
1.31612 1.33712
1.52507 1.38718
1.61787 1.44448
1.72295 1.49098
0.124944 1.62736
0.295738 1.66236
0.477847 1.63569
0.589643 1.6442
0.804853 1.61017
0.988894 1.59443
1.12152 1.65621
1.31786 1.64111
1.42638 1.64181
1.60178 1.59812
1.80711 1.561
1.86587 1.65262
0.112658 1.74304
0.363425 1.83904
0.431786 1.73893
0.673726 1.73451
0.736154 1.8076
0.945217 1.77546
1.17545 1.8017
1.25552 1.71396
1.50422 1.71882
1.58638 1.79297
1.74036 1.73077
1.89868 1.76535
0.0832031 1.85531
0.115672 1.98075
0.323824 1.97088
0.696111 1.98651
0.749649 1.88149
1.00455 1.95513
1.18425 1.87917
1.24842 1.9773
1.67562 1.86877
1.80462 1.96415
1.86156 1.88125
0.199789 2.07894
0.302011 2.11529
0.694941 2.12241
0.807469 2.08877
0.889414 1.9673
1.11571 2.06635
1.26174 2.11452
1.68422 1.97789
1.72221 2.03558
numTriangle 414
0 1 8
1 2 9
8 1 9
9 2 10
2 3 11
10 2 11
4 5 12
12 5 13
5 6 14
6 7 15
14 6 15
9 10 17
16 9 17
17 10 18
13 5 20
5 14 20
19 13 20
20 14 21
16 17 22
22 17 23
17 18 24
23 17 24
19 20 26
26 20 27
20 21 28
27 20 28
22 23 29
31 25 32
27 28 35
34 27 35
29 23 37
23 30 37
36 29 37
37 30 38
31 32 38
38 32 39
33 34 41
34 35 41
40 33 41
41 35 42
36 37 43
43 37 44
37 38 45
44 37 45
38 39 46
45 38 46
39 40 47
40 41 47
46 39 47
47 41 48
41 42 49
48 41 49
44 45 50
45 46 51
50 45 51
46 47 52
51 46 52
47 48 54
52 47 54
53 52 54
50 51 57
56 50 57
51 52 58
57 51 58
52 53 59
58 52 59
59 53 60
55 56 65
56 57 66
65 56 66
57 58 66
58 59 67
66 58 67
59 60 68
67 59 68
68 60 69
62 63 74
73 62 74
63 64 75
64 55 75
55 65 75
74 63 75
75 65 76
65 66 77
76 65 77
66 67 77
67 68 78
77 67 78
78 68 79
68 69 80
69 61 80
79 68 80
61 70 81
80 61 81
70 71 82
81 70 82
71 72 82
82 72 83
73 74 85
85 74 86
74 75 87
86 74 87
76 77 88
88 77 89
77 78 90
89 77 90
78 79 90
90 79 91
79 80 92
91 79 92
81 82 93
93 82 94
82 83 95
83 84 95
94 82 95
95 84 96
85 86 97
97 86 99
98 97 99
89 90 101
90 91 102
102 91 103
94 95 105
95 96 107
105 95 107
106 105 107
98 99 108
108 99 109
110 100 111
101 90 112
90 102 112
100 101 112
111 100 112
102 103 113
103 104 113
113 104 114
105 106 115
115 106 116

corotated_fem_pbd.py

"""
Then we use XPBD-FEM (gpu version, Jacobian solver) to simulate the deformation of 2D object
"""
import taichi as ti
from taichi.lang.ops import sqrt
from readObj import Objfile
import time

ti.init(arch=ti.gpu)

h = 0.003  # timestep size
staticPoint = 8
omega = 0.5  # SOR factor
compliance = 1.0e-4  # Fat Tissuse compliance, for more specific material,please see: http://blog.mmacklin.com/2016/10/12/xpbd-slides-and-stiffness/
alpha = compliance * (1.0 / h / h
                      )  # timestep related compliance, see XPBD paper
obj = Objfile()
# obj.read("./2dMesh.obj")
obj.readTxt("./armadillo.txt")
vertices = obj.getVertice()
triangles = obj.getFaces()

NV = obj.getNumVertice()
NF = obj.getNumFaces()  # number of faces
pos = ti.Vector.field(2, float, NV)
oldPos = ti.Vector.field(2, float, NV)
vel = ti.Vector.field(2, float, NV)  # velocity of particles
invMass = ti.field(float, NV)  #inverse mass of particles


f2v = ti.Vector.field(3, int, NF)  # ids of three vertices of each face
B = ti.Matrix.field(2, 2, float, NF)  # D_m^{-1}
F = ti.Matrix.field(2, 2, float, NF)  # deformation gradient
lagrangian = ti.field(float, NF)  # lagrangian multipliers

gradient = ti.Vector.field(2,float, 3 * NF)
dLambda = ti.field(float, NF)

gravity = ti.Vector([0, -1.2])
MaxIte = 20
NumSteps = 10

# For validation
dualResidual = ti.field(float, ())
primalResidual = ti.field(float, ())

attractor_pos = ti.Vector.field(2, float, ())
attractor_strength = ti.field(float, ())


@ti.kernel
def init_pos():
    for i in range(NV):
        pos[i] *= ti.Vector([0.4, 0.4])
        pos[i] += ti.Vector([0.1, 0.01])
        oldPos[i] = pos[i]
        vel[i] = ti.Vector([0, 0])
        invMass[i] = 1.0
    for i in range(staticPoint):
        invMass[i] = 0.0
    for i in range(NF):
        ia, ib, ic = f2v[i]
        a, b, c = pos[ia], pos[ib], pos[ic]
        B_i_inv = ti.Matrix.cols([b - a, c - a])
        B[i] = B_i_inv.inverse()


@ti.kernel
def resetLagrangian():
    for i in range(NF):
        lagrangian[i] = 0.0

@ti.func
def vec(x):
    return ti.Vector([x[0,0], x[1,0], x[0,1], x[1, 1]])

@ti.func
def computeGradient_f(idx, U, S, V):
    isSuccess = True
    sumSigma = sqrt((S[0, 0] - 1)**2 + (S[1, 1] - 1)**2)
    if sumSigma < 1.0e-6:
        isSuccess = False

    dcdS = 1.0 / sumSigma * ti.Vector([S[0, 0] - 1, 0, 0, S[1, 1] - 1
                                       ])  # (dcdS11, dcdS12, dcdS21, dcdS22)
    # Compute (dFdx)^T
    dFdxT = ti.Vector([[-B[idx][0,0]-B[idx][1,0],                         0, -B[idx][0,1]-B[idx][1,1],            0],
                      [                       0,  -B[idx][0,0]-B[idx][1,0],                        0, -B[idx][0,1]-B[idx][1,1]],
                      [             B[idx][0,0],                         0,              B[idx][0,1],            0],
                      [                       0,               B[idx][0,0],                        0,  B[idx][0,1]],
                      [             B[idx][1,0],                         0,              B[idx][1,1],            0],
                      [                       0,               B[idx][1,0],                        0,  B[idx][1,1]]])
    # Compute (dsdF)
    dF00 = ti.Vector([[1,0], [0, 0]])
    dF10 = ti.Vector([[0,0], [1, 0]])
    dF01 = ti.Vector([[0,1], [0, 0]])
    dF11 = ti.Vector([[0,0], [0, 1]])
    dsdF00 = U.transpose() @ dF00 @ V
    dsdF10 = U.transpose() @ dF10 @ V
    dsdF01 = U.transpose() @ dF01 @ V
    dsdF11 = U.transpose() @ dF11 @ V

    dsdF = ti.Matrix.cols([vec(dsdF00), vec(dsdF10), vec(dsdF01), vec(dsdF11)])

    g = dFdxT @ dsdF.transpose() @ dcdS
    g0 = ti.Vector([g[0], g[1]])

    g1 = ti.Vector([g[2], g[3]])
    g2 = ti.Vector([g[4], g[5]])
    return g0, g1, g2, isSuccess


@ti.kernel
def semiEuler():
    # semi-Euler update pos & vel
    for i in range(NV):
        if (invMass[i] != 0.0):
            vel[i] = vel[i] + h * gravity + attractor_strength[None] * (
                attractor_pos[None] - pos[i])
            oldPos[i] = pos[i]
            pos[i] = pos[i] + h * vel[i]

@ti.kernel
def updteVelocity():
    # update velocity
    for i in range(NV):
        if (invMass[i] != 0.0):
            vel[i] = (pos[i] - oldPos[i]) / h


@ti.kernel
def computeGradientVector():
    for i in range(NF):
        ia, ib, ic = f2v[i]
        a, b, c = pos[ia], pos[ib], pos[ic]
        invM0, invM1, invM2 = invMass[ia], invMass[ib], invMass[ic]
        sumInvMass = invM0 + invM1 + invM2
        if sumInvMass < 1.0e-6:
            print("wrong invMass function")
        D_s = ti.Matrix.cols([b - a, c - a])
        F[i] = D_s @ B[i]
        U, S, V = ti.svd(F[i])
        constraint = sqrt((S[0, 0] - 1)**2 + (S[1, 1] - 1)**2)
        g0, g1, g2, isSuccess = computeGradient_f(i, U, S, V)
        if isSuccess:
            l = invM0 * g0.norm_sqr() + invM1 * g1.norm_sqr(
            ) + invM2 * g2.norm_sqr()
            dLambda[i] = -(constraint + alpha * lagrangian[i]) / (l + alpha)
            lagrangian[i] = lagrangian[i] + dLambda[i]
            gradient[3 * i + 0]  = g0
            gradient[3 * i + 1]  = g1
            gradient[3 * i + 2]  = g2

@ti.kernel
def updatePos():
    for i in range(NF):
        ia, ib, ic = f2v[i]
        invM0, invM1, invM2 = invMass[ia], invMass[ib], invMass[ic]
        if (invM0 != 0.0):
            pos[ia] += omega * invM0 * dLambda[i] * gradient[3 * i + 0]
        if (invM1 != 0.0):
            pos[ib] += omega * invM1 * dLambda[i] * gradient[3 * i + 1]
        if (invM2 != 0.0):
            pos[ic] += omega * invM2 * dLambda[i] * gradient[3 * i + 2]


if __name__ == "__main__":
    # pos.from_numpy(0.2 * vertices[:,0:3:2]) # for .obj file
    pos.from_numpy(vertices)  # for txt file
    f2v.from_numpy(triangles)
    init_pos()
    pause = False
    gui = ti.GUI('Corotated FEM XPBD')
    programStart = time.time()
    sumDrawTime = 0
    drawPoint = True
    while gui.running:
        realStart = time.time()
        for e in gui.get_events():
            if e.key == gui.ESCAPE:
                gui.running = False
            elif e.key == gui.SPACE:
                pause = not pause
            elif e.key == 'p':
                drawPoint = not drawPoint
        mouse_pos = gui.get_cursor_pos()
        attractor_pos[None] = mouse_pos
        attractor_strength[None] = gui.is_pressed(gui.LMB) - gui.is_pressed(
            gui.RMB)
        gui.circle(mouse_pos, radius=15, color=0x336699)
        if not pause:
            for i in range(NumSteps):
                semiEuler()
                resetLagrangian()
                for ite in range(MaxIte):
                    computeGradientVector()
                    updatePos()
                updteVelocity()
        ti.sync()
        start = time.time()
        if not drawPoint:
            faces = f2v.to_numpy()
            for i in range(NF):
                ia, ib, ic = faces[i]
                a, b, c = pos[ia], pos[ib], pos[ic]
                gui.triangle(a, b, c, color=0x00FF00)

        positions = pos.to_numpy()
        if drawPoint:
            gui.circles(positions, radius=4, color=0xFF00)
        for i in range(staticPoint):
            gui.circle(positions[i], radius=5, color=0xFF0000)
        gui.show()
        end = time.time()
        sumDrawTime += end - start
    programEnd = time.time()
    print("Average Draw Time Ratio: ",
          sumDrawTime / (programEnd - programStart))

readObj.py

"""
Read ``.obj`` files and store vertices and triangles. 
Quad mesh is not supported right now.
Normals, texture coordinates are not extracted right now. 
"""

import numpy as np
import sys


class Objfile:
    m_numVertices = 0
    m_numFaces = 0
    m_vertices = []
    m_indices = []
    m_filename = ""
    def read(self, filename=""):
        """
            Read ``.obj`` file created by Blender
        """
        self.m_filename = filename
        if filename == "":
            print("please give me the obj file name")
            sys.exit(1)

        with open(filename,"r") as file:
            lines = [line.strip("\n") for line in file.readlines()]
            for line in lines:
                if '#' in line or 'mtl' in line or 'o' in line:
                    continue
                x = line.split(" ")
                if x[0] == 'v': # vertex
                    pos = x[1:]
                    pos = [float(p) for p in pos]
                    self.m_vertices.append(pos)
                    self.m_numVertices += 1
                elif x[0] == 'vt':
                    continue
                elif x[0] == 'vn':
                    continue
                elif x[0] == 'f': #face
                    if len(x) == 5:
                        print("please use triangle mesh")
                        sys.exit(1)
                    tri =[int(xx) for xx in x[1:]]
                    self.m_indices.append(tri)
                    self.m_numFaces += 1
            file.close()
    def readTxt(self, filename=""):
        """
            Read 2D Mesh generate by Miles Macklin's program
        """
        self.m_filename = filename
        if filename == "":
            print("please give me the obj file name")
            sys.exit(1)

        with open(filename,"r") as file:
            point = False
            triangle = False
            lines = [line.strip("\n") for line in file.readlines()]
            for line in lines:
                if 'numPoints' in line:
                    point = True
                    triangle = False
                    continue
                elif 'numTriangle' in line:
                    triangle = True
                    point = False
                    continue

                if point:
                    pos = [float(p) for p in line.split(" ")]
                    self.m_vertices.append(pos)
                    self.m_numVertices += 1
                    continue
                elif triangle:
                    tri =[int(xx) for xx in line.split(" ")]
                    self.m_indices.append(tri)
                    self.m_numFaces += 1
                    continue
            file.close()

    def ouputObjfile(self):
        print(f"{self.m_numVertices } vertices, {self.m_numFaces} faces")
        for i in range(self.m_numVertices):
            print(i, " ", self.m_vertices[i])
        for i in range(self.m_numFaces):
            print(i, " ", self.m_indices[i])

    def getVertice(self):
        return np.asarray(self.m_vertices)
    def getFaces(self):
        return np.asarray(self.m_indices)
    def getNumVertice(self):
        return self.m_numVertices
    def getNumFaces(self):
        return self.m_numFaces

if __name__ == "__main__":
    objFile = Objfile()
    # objFile.read("2dMesh.obj")
    objFile.readTxt("bunny.txt")
    vertices = objFile.getVertice()
    print(vertices)
    faces = objFile.getFaces()
    print(faces)