Barrier Aware Projected Newton in Python

[rc]

Hello,

I would like to have the Barrier Aware Projected Newton algorithm as described in your article Li2020IPC available in Python. Does the ipc-toolkit contain all the necessary ingredients or should I look elsewhere? Specifically:

a) How do I compute/update the constraint set Ĉ ←ComputeConstraintSet(x, d) - or is this done implicitly when computing the barrier function?

b) How do I update BCs and equality constraints when scripting an object motion?

(Context: about a year ago I started using ipctk in sfepy. Now I have some time to try to improve/fix the problems described in sfepy/sfepy#1132 (comment). I think the problems with no convergence might be related to not following the Algorithm 1 well enough.)

Best regards,
Robert Cimrman

[zfergus]

Hello Robert,

Thank you for your interest in the IPC Toolkit and for your work integrating it with sfepy.

While the toolkit does not provide a complete, out-of-the-box implementation of the Barrier-Aware Projected Newton method as described in the original IPC paper, it does offer the fundamental building blocks needed to implement this method within your own simulation framework.

Specifically, to address your questions:

a) Computing/Updating the Constraint Set $\hat{C}$: There are two core classes in the ipc-toolkit that facilitate this: NormalCollisions and BarrierPotential. The NormalCollisions class is responsible for detecting potential collisions based on the current configuration of your mesh and a specified barrier distance dhat. You can use its build method to update the set of active constraints $\hat{C}$ at each iteration. The BarrierPotential class computes the barrier energy and its derivatives given the current set of collisions. In your simulation loop, you would typically call C.build(...) to update the collision set and then use B(...) to compute the barrier energy and its gradient/Hessian.

b) Updating BCs and Equality Constraints: The ipc-toolkit does not directly manage boundary conditions (BCs) or equality constraints. What you want to do is create a penalty or Augmented lagrangian term for these constraints and add them to your total energy functional.

To illustrate how you might implement a simple solver loop using these components, I've included a sample code snippet below. This example demonstrates how to set up the collision detection and barrier potential, and how to perform a basic iterative update of the vertex positions while respecting the barrier constraints.

I hope this helps you better understand how to utilize the IPC Toolkit for your needs. If you have further questions or need additional assistance, please feel free to reach out.

Best,
Zach

---

import numpy as np
import scipy
import meshio
import ipctk

mesh = meshio.read("../tests/data/two-cubes-close.ply")
vertices = mesh.points
faces = mesh.cells_dict["triangle"]


class Inertia:
    """Simple inertia with prescribed target positions."""

    def __init__(self):
        import igl
        mid_x = (vertices[:, 0].min() + vertices[:, 0].max()) / 2
        self.xhat = vertices.copy()
        self.xhat[vertices[:, 0] < mid_x, 0] += 0.2
        self.xhat[vertices[:, 0] > mid_x, 0] -= 0.2
        self.mass = igl.massmatrix(
            vertices, faces, igl.MASSMATRIX_TYPE_VORONOI)
        self.mass = np.repeat(self.mass.diagonal(), 3)
        self.mass = scipy.sparse.diags(self.mass)

    def __call__(self, x):
        dx = (x - self.xhat).flatten()
        return dx.transpose() @ self.mass @ dx

    def gradient(self, x):
        return self.mass @ (x - self.xhat).flatten()

    def hessian(self, x):
        return self.mass


inertia = Inertia()

# -----------------------------------------------------------------------------

# Contact parameters
dhat = 0.001
kappa = 100 * inertia.mass.toarray().mean()

# Build the collision mesh
collision_mesh = ipctk.CollisionMesh(vertices, ipctk.edges(faces), faces)

# Build the collision set
C = ipctk.NormalCollisions()
C.use_area_weighting = True
C.use_improved_max_approximator = True
C.build(collision_mesh, vertices, dhat)

# Create a barrier potential
B = ipctk.BarrierPotential(dhat, use_physical_barrier=True)

# Initial energy and gradient
prev_energy = (
    inertia(vertices)
    + kappa * B(C, collision_mesh, vertices)
    # + ...
)
grad = (
    inertia.gradient(vertices)
    + kappa * B.gradient(C, collision_mesh, vertices)
    # + ...
)

iter = 0
while np.linalg.norm(grad) > 1e-3 and iter < 1000:
    meshio.Mesh(vertices, {"triangle": faces}).write(
        f"solver-step-{iter:03d}.ply")

    # Compute the Hessian
    hess = (
        inertia.hessian(vertices)
        + kappa * B.hessian(
            C, collision_mesh, vertices,
            project_hessian_to_psd=ipctk.PSDProjectionMethod.CLAMP)
        # + ...
    )

    # Solve for the step direction dx = -hess⁻¹ grad
    dx = scipy.sparse.linalg.spsolve(hess, -grad)
    assert grad.dot(dx) < 0

    next_vertices = vertices + dx.reshape(-1, 3, order="C")

    # ------------------------------------------------------------------------------
    # Perform a line search to find a step size that decreases the barrier potential

    # Candidates for continuous collision detection
    candidates = ipctk.Candidates()
    candidates.build(collision_mesh, vertices, next_vertices, dhat)

    # Initial step size for line search (ensure no collisions)
    alpha = candidates.compute_collision_free_stepsize(
        collision_mesh, vertices, next_vertices)

    # Backtrack until the barrier potential decreases
    condition = True
    ls_iter = 0
    while condition:
        next_vertices = vertices + alpha * dx.reshape(-1, 3, order="C")
        C.build(candidates, collision_mesh, next_vertices, dhat)
        curr_energy = (
            inertia(next_vertices)
            + kappa * B(C, collision_mesh, next_vertices)
            # + ...
        )
        condition = curr_energy > prev_energy and ls_iter < 100
        if condition:
            alpha *= 0.5
            ls_iter += 1

    vertices = next_vertices

    prev_energy = curr_energy
    grad = (
        inertia.gradient(vertices)
        + kappa * B.gradient(C, collision_mesh, vertices)
        # + ...
    )

    print("||grad||:", np.linalg.norm(grad))

    iter += 1

meshio.Mesh(vertices, {"triangle": faces}).write(f"solver-step-{iter:03d}.ply")

output.mp4

[rc]

Thank you for your fast and detailed reply! So it seems I was doing the a) OK (besides not using the area weighing and improved max approximator), but not using the Candidates at all. I will try the code snippet in my setting and think about the b) implementation.