Material Point Method Implementation Notes in Houdini

Deformation Gradient

The deformation gradient F serves as a fundamental tensor in continuum mechanics, capturing the local deformation of material points. It provides a complete description of stretch, shear, and rigid body rotation within an infinitesimal neighborhood around each material point.

Mathematical Definition:

The mapping relationship follows dx = F · dX, where dX represents an infinitesimal line segment in the reference configuration and dx represents its deformed counterpart in the current configuration.

Key Properties:


Volume Change Analysis: J = det(F)

The determinant of the deformation gradient, denoted as J, quantifies the volumetric strain. This can be decomposed into elastic and plastic components: J =Jp × Je.

The elastic component Je accounts for recoverable volume changes, while the plastic component Jp represents irreversible deformations. In typical metal plasticity scenarios, volume conservation is assumed, resulting in Jp = 1. However, for porous materials or materials undergoing significant plastic compression, Jp may deviate from unity, requiring explicit treatment of plastic volumetric strain.


Particle-Based POP Emission Source Selection

This approach derives from the official MasterClass workflow. The technique involves extracting particle sources from an MPM simulation for secondary particle emission.

A particularly useful criterion relies on monitoring the incremental change in plastic volume ΔJp. When ΔJp > 0, the material experiences volumetric expansion, indicating active plastic flow suitable for emission sources.

// Retrieve source particle data
int srcParticle = idtopoint(1, @id);

// Calculate velocity from position delta
@v = (@P - vector(point(1, "P", srcParticle))) / @TimeInc;

// Track plastic volume increment
f@dJp = @Jp - point(1, "Jp", srcParticle);

// Filter by plastic volume expansion
if(@dJp < 0.01)
    removepoint(0, @ptnum);

// Filter by velocity threshold
if(length(@v) < 10)
    removepoint(0, @ptnum);

// Filter by proximity to collision surface
float surfaceDist = 0.15;
float sdf = volumesample(2, "surface", @P);
if(sdf < -surfaceDist || sdf > surfaceDist)
    removepoint(0, @ptnum);

SDF-Based Collision Handling

Input port 1 receives the combined SDF and velocity field for collision detection.

float signedDist;
vector gradient, fieldVelocity;
    
signedDist = volumesample(1, 0, @P) - @pscale;
    
if (signedDist < 0)
{
    gradient = normalize(volumegradient(1, 0, @P));
    fieldVelocity = volumesamplev(1, "vel", @P);
    
    @P -= gradient * signedDist;

    vector particleVelocity = @v - fieldVelocity;
    float frictionCoeff = 1.0;
    float stickThreshold = 0.0;

    float normalVelocity = dot(particleVelocity, gradient);
    if (normalVelocity < stickThreshold)
    {
        vector tangentialComponent = particleVelocity - gradient * normalVelocity;
        normalVelocity = abs(normalVelocity);
        float tangentialMagnitude = length(tangentialComponent);
        
        if (tangentialMagnitude <= frictionCoeff * normalVelocity)
        {
            // Static friction mode
            particleVelocity = 0;
        }
        else
        {
            // Dynamic friction mode
            particleVelocity = tangentialComponent - (frictionCoeff * normalVelocity * tangentialComponent) / tangentialMagnitude;
        }
    
        @v = particleVelocity + fieldVelocity;
    }
}

Temporal SDF Interpolation for Multi-Substep Simulations

When executing particle-based simulations with substeps exceeding one, accurate SDF evaluation requires temporal interpolation between frames.

float timeBias = frac(@Frame);
    
vector velPrev = volumesamplev(0, "vel", @P);
vector velCurr = volumesamplev(1, "vel", @P);
vel = lerp(velPrev, velCurr, timeBias);

// Compute sample positions for both temporal samples
vector samplePosPrev = @P - timeBias * @TimeInc * vel;
vector samplePosCurr = @P + (1.0 - timeBias) * @TimeInc * vel;

// Sample SDF at interpolated positions
float sdfPrev = volumesample(0, "surface", samplePosPrev);
float sdfCurr = volumesample(1, "surface", samplePosCurr);


sdf = lerp(sdfPrev, sdfCurr, timeBias) - @pscale;

Input ports 0 and 1 should recieve velocity and SDF fields from consecutive time samples.


Computing Aggregate Orientation

For point clouds where individual particles carry orientation attributes, determining the mean orientation requires aggregation of rotational data.

First, collect all orientation values into a detail array:

p[]@orients = detail(1, "orient");

Then perform the averaging computation:

vector zAccumulator = 0;
vector yAccumulator = 0;

int sampleCount = len(p[]@orients);

for (int i = 0; i < sampleCount; i++)
{
    zAccumulator += qrotate(@orients[i], {0, 0, 1}) / float(sampleCount);
    yAccumulator += qrotate(@orients[i], {0, 1, 0}) / float(sampleCount);
}

@orient = quaternion(matrix3(maketransform(normalize(zAccumulator), normalize(yAccumulator))));

This method projects each orientation onto the Z and Y axes, averages the resulting vectors, and reconstructs a valid orientation quaternion from the normalized averages.

Tags: Houdini MPM VEX Material Point Method Physics Simulation

Posted on Sat, 09 May 2026 13:50:33 +0000 by NathanLedet