[Blender] Blender intergration

[archibate]

Concisely describe the proposed feature
I'd like to intergrate Taichi into Blender's Python scripting module so that people could play with their customized Taichi physics engine in Blender.

What we want is to create animation in Blender. For example, a ball falls down on to the surface of a cloth, the cloth breaks violently.
The cloth will encounter deformations and even mesh breaks, that means we must treat each frame as a different mesh.

Describe the solution you'd like (if any)

Solution 1

What I know about how to treat each frame as a different mesh is the Stop-motion-OBJ plugin by @neverhood311. Currently it is able to load a mesh sequence into Blender and play, with each frame a different mesh.
However, it seems to be not yet able to create a mesh without reading from disk. I asked @neverhood311 if we could create a mesh sequence from Python in neverhood311/Stop-motion-OBJ#104, inputs are welcome!
Before neverhood311/Stop-motion-OBJ#104 is resolved, I used some tricks to make it possible to create a mesh sequence from Python:

First, create a empty mesh sequence of length 250:

import taichi as ti
import numpy as np
import sys

n_frames = 250
output = sys.argv[1]

x = np.ones(1)
for frame in range(n_frames):
    print('generating frame', frame)
    writer = ti.PLYWriter(num_vertices=1)
    writer.add_vertex_pos(x, x, x)
    writer.export_frame(frame, output)

Then, import it into Blender via Stop-motion-OBJ.

Then, iterate through the mesh sequence, and change the current frame mesh name one by one:

# reblend.py

import bpy
import numpy as np


def _array(x):
    if x is not None and not isinstance(x, np.ndarray):
        x = np.array(x)
    return x


## Blender-NumPy interface

def to_numpy(b, key, dim=None):
    dim = dim or len(getattr(b[0], key))
    dim = len(getattr(b[0], key))
    seq = [0] * (len(b) * dim)
    b.foreach_get(key, seq)
    return np.array(seq).reshape(len(b), dim)


def from_numpy(b, key, a, dim=None):
    a = _array(a)
    dim = dim or len(getattr(b[0], key))
    assert len(a.shape) == 2
    assert a.shape[1] == dim
    if len(b) < a.shape[0]:
        b.add(a.shape[0] - len(b))
    seq = a.reshape(a.shape[0] * dim).tolist()
    seq = seq + [0] * (len(b) * dim - len(seq))
    b.foreach_set(key, seq)


## Blender wrappers

def object_frames(obj):
    if isinstance(obj, str):
        obj = bpy.data.objects[obj]
    return (bpy.data.meshes[x.key] for x in obj.mesh_sequence_settings.meshNameArray)


def set_object_frame(obj, i, mesh):
    if isinstance(obj, str):
        obj = bpy.data.objects[obj]
    obj.mesh_sequence_settings.meshNameArray[i].key = mesh.name


def new_mesh(name, pos, edges=[], faces=[], uv=None):
    pos = _array(pos)
    edges = _array(edges)
    faces = _array(faces)
    uv = _array(uv)
    mesh = bpy.data.meshes.new(name)
    mesh.from_pydata(pos.tolist(), edges.tolist(), faces.tolist())
    if uv is not None:
        mesh.uv_layers.new()
        from_numpy(mesh.uv_layers.active.data, 'uv', uv)
    return mesh


def new_object(name, mesh, col=None):
    obj = bpy.data.objects.new(name, mesh)
    if col is not None:
        col.objects.link(obj)
    bpy.context.view_layer.objects.active = obj


## NumPy helpers

def meshgrid(n):
    def _face(x, y):
        return np.array([(x, y), (x + 1, y), (x + 1, y + 1), (x, y + 1)])

    def _edge1(x, y):
        return np.array([(x, y), (x + 1, y)])

    def _edge2(x, y):
        return np.array([(x, y), (x, y + 1)])

    n_particles = n**2
    n_edges1 = (n - 1) * n
    n_edges2 = (n - 1) * n
    n_faces = (n - 1)**2
    xi = np.arange(n)
    yi = np.arange(n)
    xs = np.linspace(0, 1, n)
    ys = np.linspace(0, 1, n)
    pos = np.array(np.meshgrid(xs, ys)).swapaxes(0, 2).reshape(n_particles, 2)
    faces = _face(*np.meshgrid(xi[:-1], yi[:-1])).swapaxes(0, 1).swapaxes(1, 2).swapaxes(2, 3)
    faces = (faces[1] * n + faces[0]).reshape(n_faces, 4)
    edges1 = _edge1(*np.meshgrid(xi[:-1], yi)).swapaxes(0, 1).swapaxes(1, 2).swapaxes(2, 3)
    edges2 = _edge2(*np.meshgrid(xi, yi[:-1])).swapaxes(0, 1).swapaxes(1, 2).swapaxes(2, 3)
    edges1 = (edges1[1] * n + edges1[0]).reshape(n_edges1, 2)
    edges2 = (edges2[1] * n + edges2[0]).reshape(n_edges2, 2)
    edges = np.concatenate([edges1, edges2], axis=0)
    pos = np.concatenate([pos, np.zeros((n_particles, 1))], axis=1)
    uv = pos[faces, :2].reshape(faces.shape[1] * faces.shape[0], 2)
    return pos, edges, faces, uv


__all__ = '''
np
bpy
from_numpy
to_numpy
object_frames
new_object
new_mesh
set_object_frame
meshgrid
'''.strip().splitlines()

---

from reblend import *
import taichi as ti
import taichi_glsl as tl

ti.init(arch=ti.cpu)

dt = 2e-4
steps = 90
stiffness = 2e3
damping = 1e-2
resistance = 0.3
gravity = 2
tolerance = 0.69
tol_randomness = 0.3
mass_scale = 0.27
cloth_height = 1.5
ball_radius = 0.5
ground_height = -0.5

cloth = bpy.data.collections['Cloth'].objects[0].data
_pos = to_numpy(cloth.vertices, 'co').astype(np.float32)
_edges = to_numpy(cloth.edges, 'vertices').astype(np.int32)
_faces = to_numpy(cloth.polygons, 'vertices').astype(np.int32)

ball = bpy.data.collections['Ball'].objects[0].data
ball_vertices = to_numpy(ball.vertices, 'co').astype(np.float32)
ball_edges = to_numpy(ball.edges, 'vertices').astype(np.int32)
ball_faces = to_numpy(ball.polygons, 'vertices').astype(np.int32)

_pos[:, 2] += cloth_height
pos = ti.Vector.field(3, float, _pos.shape[0])
vel = ti.Vector.field(3, float, _pos.shape[0])
edges = ti.Vector.field(2, int, _edges.shape[0])
rest = ti.field(float, _edges.shape[0])
ball_pos = ti.Vector.field(3, float, ())
ball_vel = ti.Vector.field(3, float, ())

ball_pos[None] = [0.2, -0.06, 2.8]
pos.from_numpy(_pos)
edges.from_numpy(_edges)


@ti.kernel
def init():
    for e in edges:
        p, q = edges[e]
        disp = pos[q] - pos[p]
        rest[e] = disp.norm()

@ti.kernel
def substep():
    for e in edges:
        if rest[e] >= 1e3:
            continue
        p, q = edges[e]
        disp = pos[q] - pos[p]
        disv = vel[q] - vel[p]
        k = disp.norm() - rest[e]
        if k > rest[e] * tolerance * (1 + (ti.random() * 2 - 1) * tol_randomness):
            rest[e] = 1e4
        acc = disp * k / rest[e]**2
        acc += disv * damping
        vel[p] += stiffness * acc * dt
        vel[q] -= stiffness * acc * dt
    for p in pos:
        vel[p].z -= gravity * dt
        new_vel = tl.ballBoundReflect(pos[p], vel[p], ball_pos[None], ball_radius + 0.037, 6)
        dv = vel[p] - new_vel
        ball_vel[None] += dv * mass_scale / pos.shape[0]
        vel[p] = new_vel
        if pos[p].z <= ground_height and vel[p].z < 0:
            vel[p].z = 0
        if pos[p].z == cloth_height:
            if (pos[p].x == -1 or pos[p].x == 1) or (pos[p].y == -1 or pos[p].y == 1):
                vel[p] *= 0
        vel[p] *= ti.exp(-dt * resistance)
        pos[p] += vel[p] * dt
    ball_vel[None].z -= gravity * dt
    ball_pos[None] += ball_vel[None] * dt
    if ball_pos[None].z <= ball_radius + ground_height and ball_vel[None].z < 0:
        ball_vel[None].z *= -0.6


def get_edges():
    cols = np.extract(_rest >= 1e3, np.arange(_rest.shape[0]))
    return np.delete(_edges, cols, axis=0)

def get_faces():
    endpoint0 = _pos[_faces[:, 0]]
    endpoint1 = _pos[_faces[:, 1]]
    endpoint2 = _pos[_faces[:, 2]]
    endpoint3 = _pos[_faces[:, 3]]
    disp0 = np.sum((endpoint0 - endpoint1)**2, axis=1)
    disp1 = np.sum((endpoint1 - endpoint2)**2, axis=1)
    disp2 = np.sum((endpoint2 - endpoint3)**2, axis=1)
    disp3 = np.sum((endpoint3 - endpoint0)**2, axis=1)
    disp = np.maximum(disp0, np.maximum(disp1, disp2, disp3))
    cols = np.extract(disp >= 1.04 * ((1 + tolerance) * _rest_average)**2, np.arange(_faces.shape[0]))
    return np.delete(_faces, cols, axis=0)
    
init()
_rest = rest.to_numpy()
_rest_average = np.average(_rest)
for i, (cloth, ball) in enumerate(zip(object_frames('p_sequence'), object_frames('q_sequence'))):
    print('rendering frame', i)
    for s in range(steps):
        substep()
    _pos = pos.to_numpy()
    _rest = rest.to_numpy()
    cloth = new_mesh('p_seq', _pos, get_edges(), get_faces())
    set_object_frame('p_sequence', i, cloth)
    ball = new_mesh('q_seq', ball_vertices + ball_pos.to_numpy().reshape(1, 3), ball_edges, ball_faces)
    set_object_frame('q_sequence', i, ball)

However, this solution seems to have 3 downside:

1. Serious memory leakage, it seems to be filling up all my 16GB memory with ~20 replay.
2. Cannot reload the generated meshes after reopening .blend file, it is stationary at the current frame.
3. Can't do other things during the rendering process, the Blender is stucked until the Python script is done.

Solution 2

It seems to me that Blender already have several built-in physics solver.
Is that possible to make taichi a blender plugin and become one of these solvers, like taichi_elements did?
@PavelBlend Would you share some insights on how to add a user-written solver (instead of our specific MPMSolver) to Blender?
It should be not only able to do the particle simulations, but also mesh cloth simulations like mentioned above.

Additional context
@PavelBlend @neverhood311 @yuanming-hu What's your opinion on these solutions? Inputs are welcome!

[PavelBlend]

Hello, @archibate

I think that you can create one mesh object and update its geometry when the frame changes.
To do this, just add a function to update the geometry in handlers:

bpy.app.handlers.frame_change_pre.append(my_function)