Starting with an established engine is a good idea.

I’ve written something akin to what you’re looking for from scratch and it’s a project unto itself.

Honestly I would start 2d. 3d more than doubles the complexity.

I’ve been out of the loop for a good long while, but Box2d would be my first port of call. (They also recently launched Box3d, so maybe the interaction would be familiar once you’ve finished with 2d)

Yeah … I’d just go with Blender, honestly.

What you’re describing is actually already extremely close to how it processes cloth/jiggle physics internally. All you’d need to do is create your array of points, link them with lines, and then set physics parameters for each one. Plus, it already comes with ways to simulate external forces on the system, if you want.

It’s actually quite simple to simulate with python or octave, even with a simple update scheme. Take for instance a look at https://interactivecomputergraphics.github.io/physics-simulation/, especially its “mass spring system” as a starting point.

Here there are examples on how simulate various systems from scratch.

Oh, I read it wrong before, but basically the same principle of my previous answer, but the differential equation between dots, so there is a direction unit vector associated with the two connected points, the differential equation would need the spring length too, but that would be easy to put, as well as the constants.

The old wrong answer:

You can model this grid pattern using partial differential equations, the springs + masses setup reduces to a 2D wave equation, your varying masses and spring constants reduce to a same order and size matrix of constants, instead a global value, the values are on the grid. With finite differences, it’s very easy to solve.

It will be just a grid updating, you can build this in less than 100 lines of code:

import numpy as np


def free_wave_equation_2d(lenx, leny, sol_len, c, ew, ew_len, dt):
    """
    Parameters:
    lenx      : int          - Number of points in x direction
    leny      : int          - Number of points in y direction
    sol_len   : int          - Number of time steps
    c         : ndarray (2D): Wave speed coefficient array
    ew        : ndarray (1D): Source/initial condition array
    ew_len    : int          - Length of the ew array
    dt        : float        - Time step

    Returns:
    array     : ndarray (3D) - Resulting solution (lenx, leny, sol_len)
    """

    array = np.zeros((lenx, leny, sol_len))
    acceleration = np.zeros((lenx, leny))
    half_lenx = lenx // 2
    half_leny = leny // 2

    # Set initial condition at t=0
    array[half_lenx, half_leny, 0] = ew[0]

    for i in range(1, sol_len):
        # u_curr represents u(t)
        u_curr = array[:, :, i - 1]

        acceleration[:, :] = 0

        # --- Finite Difference Laplacian Calculation ---
        # X-direction
        acceleration[0, :] = u_curr[1, :] - u_curr[0, :]
        acceleration[-1, :] = u_curr[-2, :] - u_curr[-1, :]
        acceleration[1:-1, :] = u_curr[2:, :] + u_curr[:-2, :] - 2 * u_curr[1:-1, :]

        # Y-direction
        acceleration[:, 0] += u_curr[:, 1] - u_curr[:, 0]
        acceleration[:, -1] += u_curr[:, -2] - u_curr[:, -1]
        acceleration[:, 1:-1] += (
            u_curr[:, 2:, :] + u_curr[:, :-2, :] - 2 * u_curr[:, 1:-1, :]
        )

        # Apply wave speed coefficient
        acceleration = acceleration * (c**2)

        # --- Finite Difference Integration Step ---
        if i == 1:
            # For the very first step
            array[:, :, i] = u_curr + 0.5 * acceleration * (dt**2)
        else:
            # Standard Leapfrog / Central Difference Integration:
            u_prev = array[:, :, i - 2]
            array[:, :, i] = 2 * u_curr - u_prev + acceleration * (dt**2)

        if i < ew_len:
            array[half_lenx, half_leny, i] = ew[i]

    return array

Edit: The code was ported from my college time fortran code.