The source code can be cloned or downloaded from the Github repository.

The Book and the Code


The first edition of the book corresponds to the code version 1.0.0. After the release of v1.0.0, the code will evolve, and the book-version of the code will be maintained as a separate branch book-1st-edition in the GitHub repository. It will take bug fixes only, but no feature addition or API changes to keep in sync with the book.


Although I tried not to diverge from the code in the book, several important enhancement and bug fixes were applied to the code after finishing the manuscript for the book. Below is the list of the notable differences between the book and v1.0.0.

Page 63~67

The Surface3 interface has changed. The main idea of this update is to separate the world and local coordinates, and let the deriving classes to override the local version only. For example, the book shows closestPoint is virtual. The new code adds closestPointLocal as a new virtual function to be overriden and make closestPoint to be non-virtual.

Page 69

The ImplicitSurface3 interface has changed. Similar to the update above, this modification also let the sub-classes to override the function in local frame (signedDistanceLocal), but keep the world frame version as a non-virtual function.

Page 120~

The member variable _mass in SphSystemData3 is removed. The new code uses mass() instead.

Page 184

There was a bug in the code. The fixed version is shown below:

auto u = velocity->uAccessor();
auto uPos = velocity->uPosition();

velocity->parallelForEachUIndex([&](size_t i, size_t j, size_t k) {
    Vector3D pt = uPos(i, j, k);
    if (isInsideSdf(_colliderSdf.sample(pt))) {
        Vector3D colliderVel = collider()->velocityAt(pt);
        Vector3D vel = velocity->sample(pt);
        Vector3D g = _colliderSdf.gradient(pt);
        if (g.lengthSquared() > 0.0) {
            Vector3D n = g.normalized();
            Vector3D velr = vel - colliderVel;
            Vector3D velt = projectAndApplyFriction(
                velr, n, collider()->frictionCoefficient());

            Vector3D velp = velt + colliderVel;
            uTemp(i, j, k) = velp.x;
        } else {
            uTemp(i, j, k) = colliderVel.x;
    } else {
        uTemp(i, j, k) = u(i, j, k);

u.parallelForEachIndex([&](size_t i, size_t j, size_t k) {
    u(i, j, k) = uTemp(i, j, k);

Page 213~

The interface GridPressureSolver3::solve has changed to:

virtual void solve(
    const FaceCenteredGrid3& input,
    double timeIntervalInSeconds,
    FaceCenteredGrid3* output,
    const ScalarField3& boundarySdf
        = ConstantScalarField3(kMaxD),
    const VectorField3& boundaryVelocity
        = ConstantVectorField3({0, 0, 0}),
    const ScalarField3& fluidSdf
        = ConstantScalarField3(-kMaxD)) = 0;