Author: singhishree

The Log Diffusion Solver

Post author By singhishree
Post date October 20, 2025

Teaser Figure — **Figure 1:** Comparison of distance estimation by heat methods. Standard diffuse-then-log (middle) suffers from numerical issues when the heat is small (large geodesic distances), leading to noisy and inaccurate isolines. Our proposed log-space diffusion (right) performs diffusion directly in the log-heat domain, producing smooth and accurate estimates, closely matching the ground-truth distances (left). It is worth noting that our solver outperforms others primarily for smaller time step values.

Introduction

The heat equation is one of the canonical diffusion PDEs: it models how “heat” (or any diffusive quantity) spreads out over space as time progresses. Mathematically, let $u(x,t) \in \mathbb{R}_{\geq 0}$ denote the temperature at spatial location $x \in \mathbb{R}^{d}$ and time $t \in \mathbb{R}$. The continuous heat equation is given by

\[\frac{\partial u}{\partial t} = \Delta u,\]

where $\Delta u (x, t) = \sum_{i=1}^{d} \frac{\partial^2 u}{\partial x_i^2}$ is the Laplace operator representing how local temperature gradients drive heat from hotter to colder regions. On a surface mesh, the spatial domain is typically discretized so that the PDE reduces to a system of ODEs in time, which can be integrated using stable and efficient schemes such as backward Euler. Interestingly, in many applications—including geodesic distance computation ¹ and optimal transport ²—it is often more useful to work with the logarithm of the heat rather than the heat itself. For example, Varadhan’s formula states that the geodesic distance $\phi$ between two points $x$ and $x^*$ on a curved domain can be recovered from the short-time asymptotics of the heat kernel:

\[\phi(x, x^{*}) = \lim_{t \to 0} \sqrt{-4t \,\log k_{t,x^{*}}(x)},\]

where $k_{t,x^{*}}(x) = u(x,t)$ is the heat distribution at time $t$ resulting from an initial Dirac delta at $x^*$, i.e., $u(x,0) = \delta(x-x^{*})$. In practice, this leads to a two-step workflow: (1) simulate the heat equation to obtain $u_t$, and (2) take the logarithm to recover the desired quantity. However, this approach can be numerically unstable in regions where $u$ is very small: tiny relative errors in $u$ can be greatly magnified in the log domain. For example, an absolute error of $10^{-11}$ in $u$ becomes an error of roughly $25.33$ in $\log u$.

A natural idea is therefore to evolve the logarithm of the heat directly. Let $u>0$ satisfy the heat equation and define $w = \log u$. One can show that $w$ obeys the nonlinear PDE ³

\[\frac{\partial w}{\partial t} = \Delta w + \|\nabla w\|_2^2.\]

This formulation has the appealing property that it avoids taking the logarithm after simulation, thereby reducing the amplification of small relative errors in $u$ when $u$ is tiny. However, the price of this reformulation is the additional nonlinear term $\|\nabla w\|_2^2$, which complicates time integration as each step requires iterative optimization, and introduces a discrepancy with the discretized heat equation. In practice, evaluating $\|\nabla w\|_2^2$ requires a discrete gradient operator, and the resulting ODE is generally not equivalent to the one obtained by first discretizing the heat equation for $u$ (which involves only a discrete Laplacian) up to a log transform. This non-equivalence stems from (1) there is a new spatially discretized gradient operator, and (2) the fact that the logarithm is nonlinear and does not commute with discrete differential operators.

Motivated by this discrepancy, we take a different approach that preserves the exact relationship to the standard discretization. We first discretize the heat equation in space, yielding the familiar linear ODE system for $u$. We then apply the transformation $w = \log u$ to this discrete system, producing a new ODE for $w$. Now, by construction, the new ODE is equivalent to the heat ODE up to a log transform of variable. This still allows us to integrate directly in log-space, which allows more accurate and stable computation. We show that this formulation offers improved numerical robustness in regions where $u$ is very small, and demonstrate its effectiveness in applications such as geodesic distance computation.

**Table 1:** Side-by-side procedures for computing $\log u$ in heat diffusion, note the highlighted steps are the same for standard heat diffusion and ours.
Standard (post-log)	Continuous log-PDE	Proposed discrete-log
1. Start with the heat PDE for $u$.	1. Start with the log-heat PDE for $w=\log u$.	1. Start with the heat PDE for $u$.
2. Discretize in space $\Rightarrow$ ODE for $u$.	2. Discretize in space for $\nabla$ and $\Delta$ $\Rightarrow$ ODE for $w$.	2. Discretize in space $\Rightarrow$ ODE for $u$.
3. Time-integrate in $u$-space.	3. Time-integrate in $w$-space and output $w$.	3. Apply change of variables $w=\log u$ $\Rightarrow$ ODE for $w$.
4. Compute $w=\log u$ as a post-process.		4. Time-integrate in $w$-space and output $w$.

Background

We begin by reviewing how previous approaches are carried out, focusing on the original heat method ¹ and the log-heat PDE studied in ³. In particular, we discuss two key stages of the discretization process: (1) spatial discretization, where continuous differential operators such as the Laplacian and gradient are replaced by their discrete counterparts, turning the PDE into an ODE; and (2) time discretization, where a numerical time integrator is applied to solve the resulting ODE. Along the way, we introduce notation that will be used throughout.

Heat Method

Spatial discretization. In the heat method ¹, the heat equation is discretized on triangle mesh surfaces by representing the heat values at mesh vertices as a vector $u \in \mathbb{R}^{|V|}$, where $|V|$ is the number of vertices. A corresponding discrete Laplace operator is then required to evolve $u$ in time. Many Laplacian discretizations are possible, each with its own trade-offs ⁴. The widely used cotangent Laplacian is written in matrix form as $M^{-1}L$, where $M, L \in \mathbb{R}^{|V| \times |V|}$. The matrix $L$ is defined via

\[L_{ij} = \begin{cases} \frac{1}{2}(\cot \alpha_{ij} + \cot \beta_{ij}) \quad & i \neq j \\ -\sum_{i \neq j}L_{ij} & i = j \end{cases},\]

where $i,j$ are vertex indices and $\alpha_{ij}, \beta_{ij}$ are the two angles opposite the edge $(i,j)$ in the mesh. The mass matrix $M = \mathrm{diag}(m_1, \dots, m_{|V|})$ is a lumped (diagonal) area matrix, with $m_i$ denoting the area associated with vertex $i$. Note that $L$ has several useful properties: it is sparse, symmetric, and each row sums to zero. This spatial discretization converts the continuous PDE into an ODE for the vector-valued heat function $u(t)$:

\[\frac{\mathrm{d}u}{\mathrm{d}t} = M^{-1}L\,u.\]

Time discretization and integration. The ODE is typically advanced in time using a fully implicit backward Euler scheme, which is unconditionally stable. For a time step $\Delta t = t_{n+1} – t_n$, backward Euler yields:

\[(M – \Delta t\,L)\,u_{n+1} = M\,u_n,\]

where $u_n$ is given and we solve for $u_{n+1}$. Thus, each time stepping requires solving a sparse linear system with system matrix $M – \Delta t\,L$. Since this matrix remains constant throughout the simulation, it can be factored once (e.g., by Cholesky decomposition) and reused at every iteration, greatly improving efficiency.

Log Heat PDE

Spatial discretization. In prior work ³ on the log-heat PDE, solving this requires discretizing not only the Laplacian term but also the additional nonlinear term $\|\nabla w\|_2^{2}$. After spatial discretization, the squared gradient norm of a vertex-based vector $w \in \mathbb{R}^{|V|}$ is computed as

\[g(w) = W \left( \sum_{k=1}^{3} (G_k w) \odot (G_k w) \right) \in \mathbb{R}^{|V|}\]

where $G_{k} \in \mathbb{R}^{|F| \times |V|}$ is the discrete gradient operator on mesh faces for the $k$-th coordinate direction, and $W \in \mathbb{R}^{|V| \times |F|}$ is an averaging operator that maps face values to vertex values. With this discretization, the ODE for $w(t) \in \mathbb{R}^{|V|}$ becomes

\[\frac{\mathrm{d}w}{\mathrm{d}t} = M^{-1}L\,w + g(w)\]

Time discretization and integration. If we apply fully implicit Euler, we would find $w^{n+1}$ by solving a nonlinear equation. To avoid this large nonlinear solve, prior work ³ splits the problem into two substeps:

a linear diffusion part: $\frac{\mathrm{d}w}{\mathrm{d}t} = M^{-1}L\,w$, which is stiff but easy to solve implicitly,
a nonlinear part: the full equation, which is nonlinear but convex.

This separation lets each component be treated with the method best suited for it: a direct linear solve for the diffusion, and a convex optimization problem for the nonlinear term. Using Strang-Marchuk splitting:

Half-step diffusion: $(M – \frac{\Delta t}{2} L)\,w^{(1)} = Mw^n$
Full-step nonlinear term: Solve a convex optimization problem to find $w^{(2)}$
Half-step diffusion: $(M – \frac{\Delta t}{2} L)\,w^{(3)} = Mw^{(2)}$, and set $w^{n+1} = w^{(3)}$

A New ODE in Log-Domain

Let us revisit the spatially discretized heat (or diffusion) equation. For each vertex $i \in \{1, \dots |V|\}$ on a mesh, its evolution is described by

\[\frac{\mathrm{d} u_i}{\mathrm{d} t} = \frac{1}{m_i} (L u)_i\]

where $u$ represents our heat, and $L$ is the discrete Laplacian operator. A key insight is that the Laplacian term can be rewritten as a sum of local differences between a point and its neighbors:

\[(Lu)_i = \sum_{j \sim i}L_{ij} (u_j – u_i).\]

This shows that the change in $u_i$ depends on how different it is from its neighbors, weighted by the Laplacian matrix entries. Now, before doing time discretization, we use a change of variable $w \Leftarrow\log u$, accordingly, $u = \mathrm{e}^{w}$, to obtain the ODE

\[\mathrm{e}^{w_i}\frac{\mathrm{d} w_i}{\mathrm{d}t} = \frac{1}{m_i} \sum_{j \sim i}L_{ij} (\mathrm{e}^{w_j} – \mathrm{e}^{w_i}),\]

which can be rearranged into

\[\frac{\mathrm{d} w_i}{\mathrm{d}t} = \frac{1}{m_i} \sum_{j \sim i}L_{ij} (\mathrm{e}^{w_j \, – \, w_i} – 1).\]

In vectorized notation, we construct a matrix $\tilde{L} \Leftarrow L – \mathrm{diag}(L_{11}, \dots, L_{|V| \,|V|})$, which zeros out the diagonals of $L$. Then the log heat $w(t) \in \mathbb{R}^{|V|}$ follows the ODE

\[\frac{\mathrm{d} w}{\mathrm{d}t} = M^{-1} \left( \mathrm{diag}(\mathrm{e}^{-w})\, \tilde{L}\,\mathrm{e}^{w} – \tilde{L} 1 \right),\]

and another equivalent way to put it is

\[\frac{\mathrm{d} w}{\mathrm{d}t} = M^{-1} \left( \left(\tilde{L} \odot \exp(1 w^{\top} – w 1^{\top} ) \right) 1 – \tilde{L} 1 \right).\]

This ODE is equivalent to the standard heat ODE up to a log transform of the integrated results, but crucially it avoids the amplification of small relative errors in $u$ when $u$ is tiny by integrating directly in log-space.

Numerical Solvers for Our New ODE

Forward (Explicit Euler)

The most straightforward approach uses explicit Euler to update the solution at each time step. The update rule is

\[w_{n+1} = w_{n} + \Delta t\, M^{-1} \left( \mathrm{diag}(\mathrm{e}^{-w_n})\, \tilde{L}\,\mathrm{e}^{w_n} – \tilde{L} 1 \right)\]

While this is efficient and can be solved linearly, it causes instability.

Semi-Implicit

We can discretize with both explicit and implicit parts by decomposing the exponential function as a linear and nonlinear part:

\[\mathrm{e}^{x} = \underbrace{(1 + x)}_{\text{implicit}} + \underbrace{(\mathrm{e}^{x} – x – 1)}_{\text{explicit}}\]

Applying this decomposition yields a linear system for solving $w_{n+1}$:

\[\left( M + \Delta t \left(\tilde{L} – \mathrm{diag}(\tilde{L} 1) \right) \right) w_{n+1} = M w_{n} – \Delta t \left( \mathrm{diag}(\mathrm{e}^{-w_n}) \tilde{L} \mathrm{e}^{w_{n}} – \tilde{L} w_{n} + \mathrm{diag}(w_{n}) \tilde{L} 1 – \tilde{L} 1 \right)\]

This retains efficiency but does not fully address instability issues.

Backward (Convex Optimization)

Solving the backward Euler equation involves finding the vector $w^{t+1}$ that satisfies

\[\frac{w^{t+1}_i – w^{t}_i}{\Delta t} = \frac{1}{m_i} \sum_{j \sim i}L_{ij} (\mathrm{e}^{w^{t+1}_j \, – \, w^{t+1}_i} – 1)\]

Taking inspiration from prior work ³, we find that relaxing this equality into an inequality leads to a convex constraint on $w^{t+1}_i$:

\[\frac{w^{t+1}_i – w^{t}_i}{\Delta t} \geq \frac{1}{m_i} \sum_{j \sim i}L_{ij} (\mathrm{e}^{w^{t+1}_j \, – \, w^{t+1}_i} – 1)\]

It can be shown that minimizing the objective function

\[f(w) = \sum_{i}M_iw_i\]

with the inequality constraint defined above yields the exact solution to the original equation. We have written a complete proof for this finding here. We implement this constrained optimization problem using the Clarabel solver ⁶ in CVXPY ⁷, a Python library for convex optimization problems.

Distance computations by different solvers — **Figure 2:** Visualization of distance computations on the mesh by different solvers. Ground Truth, HeatBackward, LogHeatBackward, NewLogHeatBackward (from left to right)

Experiments

As a common application, we demonstrate the effectiveness of log heat methods for computing fast approximations of geodesic distances. As ground truth, we used libigl’s ⁵ implementation of exact discrete geodesic distances. We used the mesh of a kitten with 10 diffusion time-steps equally spaced between $t=0$ and $t=10^{-3}$.

**Figure 3:** Correlation plot for log heat methods against an exact geodesic computation.

**Table 2:** Correlation analysis results.
Mesh	# of verts	Method	Pearson $r$ ↑	RMSE ↓	MAE ↓
Cat	9447	HeatBackward	0.986	0.285	0.241
		LogHeatBackward	-0.183	0.577	0.522
		NewLogHeatBackward (ours)	0.879	0.163	0.112

The flat region for the NewLogHeatBackward Solver probably indicates regions of numerical overflow/underflow, since the log of tiny heat values is going to be incredibly negative. Through our experimentation we observed that the solver for this our formulated PDE outperforms others mostly in cases where the diffusion time-step is extremely small.

References

[1] Keenan Crane, Clarisse Weischedel, and Max Wardetzky. “The Heat Method for Distance Computation.” Communications of the ACM, vol. 60, no. 11, October 2017, pp. 90–99.

[2] Justin Solomon, Fernando de Goes, Gabriel Peyré, Marco Cuturi, Adrian Butscher, Andy Nguyen, Tao Du, and Leonidas Guibas. “Convolutional Wasserstein Distances: Efficient Optimal Transportation on Geometric Domains.” ACM Transactions on Graphics, vol. 34, no. 4, August 2015.

[3] Leticia Mattos Da Silva, Oded Stein, and Justin Solomon. “A Framework for Solving Parabolic Partial Differential Equations on Discrete Domains.” ACM Transactions on Graphics, vol. 43, no. 5, October 2024.

[4] Max Wardetzky, Saurabh Mathur, Felix Kälberer, and Eitan Grinspun. “Discrete Laplace Operators: No Free Lunch.” In Proceedings of the Fifth Eurographics Symposium on Geometry Processing, pages 33–37, Barcelona, Spain, 2007.

[5] Alec Jacobson et al. “libigl: A Simple C++ Geometry Processing Library.” Available at https://github.com/libigl/libigl, 2025.

[6] Paul J. Goulart and Yuwen Chen. “Clarabel: An Interior-Point Solver for Conic Programs with Quadratic Objectives.” arXiv preprint arXiv:2405.12762, 2024.

[7] Steven Diamond and Stephen Boyd. “CVXPY: A Python-Embedded Modeling Language for Convex Optimization.” Journal of Machine Learning Research, vol. 17, no. 83, pages 1–5, 2016.

Post author By singhishree
Post date October 20, 2025
Tags diffusion, pde, solver

Research

Gaussian Fluids

Post author By singhishree
Post date September 26, 2025

SGI Fellows: Haojun Qiu, Nhu Tran, Renan Bomtempo, Shree Singhi and Tiago Trindade
Mentors: Sina Nabizadeh and Ishit Mehta
SGI Volunteer: Sara Samy

Introduction

The field of Computational Fluid Dynamics (CFD) has long been a cornerstone of both stunning visual effects in movies and video games, and also a critical tool in modern engineering, for designing more efficient cars and aircrafts for example.
At the heart of any CFD application there exists a solver that essentially predicts how a fluid moves in space at discrete time intervals. The dynamic behavior of a fluid is governed by a famous set of partial differential equations, called the Navier-Stokes equations, here presented in vector form:

$$
\begin{align}
\frac{\partial \mathbf{u}}{\partial t}
+ (\mathbf{u} \cdot \nabla) \mathbf{u} &=
-\frac{1}{\rho} \nabla p
+ \nu \nabla^{2} \mathbf{u}
+ \mathbf{f},
\quad
\text{subject to}\quad \nabla \cdot \mathbf{u}
&=
0,
\end{align}
$$

where $t$ is time, $\mathbf{u}$ is the varying velocity vector field, $\rho$ is the fluids density, $p$ is the pressure field, $\nu$ is the kinematic viscosity coefficient and $\mathbf{f}$ represents any external body forces (e.g. gravity, electromagnetic, etc.). Let’s quickly break down what each term means in the first equation, which describes the conservation of momentum:

$\dfrac{\partial \mathbf{u}}{\partial t}$ is the local acceleration, describing how the velocity of the fluid changes at a fixed point in space.
$(\mathbf{u} \cdot \nabla) \mathbf{u}$ is the advection (or convection) term. It describes how the fluid’s momentum is carried along by its own flow. This non-linear term is the source of most of the complexity and chaotic behavior in fluids, like turbulence.
$\frac{1}{\rho}\nabla p$ is the pressure gradient. It’s the force that causes fluid to move from areas of high pressure to areas of low pressure.
$\nu \nabla^{2} \mathbf{u}$ is the viscosity or diffusion term. It accounts for the frictional forces within the fluid that resist flow and tend to smooth out sharp velocity changes. For the “Gaussian Fluids” paper, this term is ignored ($\nu=0$) to model an idealized inviscid fluid.

The second equation, $\nabla \cdot \mathbf{u} = 0$, is the incompressibility constraint. It mathematically states that the fluid is divergence-free, meaning it cannot be compressed or expanded.

The first step to numerically solving these equations on some spatial domain $\Omega$ using classical solvers is to discretize this domain. To do that, two main approaches have dominated the field from its very beginning:

Eulerian methods: which use fixed grids that partition the spatial domain into a (un)structured collection of cells on which the solution is approximated. However, these methods often suffer from “numerical viscosity,” an artificial damping that smooths away fine details;
Lagrangian methods: which use particles to sample the solution at certain points in space and that move with the fluid flow. However, these methods can struggle with accuracy and capturing delicate solution structures.

The core problem has always been a trade-off. Achieving high detail with traditional solvers often requires a massive increase in memory and computational power, hitting the “curse of dimensionality”. Hybrid Lagrangian-Eulerian methods that try to combine the best of both worlds can introduce their own errors when transferring data between grids and particles. The quest has been for a representation that is expressive enough to capture the rich, chaotic nature of fluids without being computationally prohibitive.

A Novel Solution: Gaussian Spatial Representation

By now you have probably come across an emerging technology in the field of 3D rendering called 3D Gaussian Splatting (3DGS). It attempts to represent a 3D scene using, not traditional polygonal meshes or volumetric voxels, but in a manner similar to an artist using paint brush strokes to approximate a colored picture, Gaussian Splatting uses 3D gaussians to “paint” a 3D scene. These gaussians are nothing more than a ellipsoid with a color value and opacity associated.

Think of each Gaussian not as a solid, hard-edged ellipsoid, but as a soft, semi-transparent cloud (Fig.2). It’s most opaque and dense at its center and it smoothly fades to become more transparent as you move away from the center. The specific shape, size, and orientation of this “cloud” are controlled by its parameters. By blending thousands or even millions of these colored, transparent splats together, Gaussian Splatting can represent incredibly detailed and complex 3D scenes.

Fig.1 Example 3DGS of a Lego build (left) and a zoomed in view of the gaussians that make it up (right).

Drawing from this incredible effectiveness of using gaussians to encode complex geometry and lighting of 3D scenes, a new paper presented at SIGGRAPH 2025 titled “Gaussian Fluids: A Grid-Free Fluid Solver based on Gaussian Spatial Representation” by Jingrui Xing et al. introduces a fascinating new approach to computational domain representation for CFD applications.

The new method, named Gaussian Spatial Representation (GSR), uses gaussians to essentially “paint” vector fields. While 3DGS encodes local color information into each gaussian, GSR encodes local vector fields. In this way each gaussian stores a vector that defines a local direction of the velocity field, where the magnitude of the vectors is greater closer the gaussian’s center, and quickly gets smaller when moving farther from it. In Fig.2 we can see an example of 2D gaussian viewed with color data (left), as well as with a vector field (right).

Fig.2 Example plot of a 2D Gaussian with scalar data (left) and vector data (right)

By splatting a lot of these gaussians and adding the vectors we can then define a continuously differentiable vector field around a given domain, which allows us to solve the complex Navier-Stokes equations not through traditional discretization, but as a physics-based optimization problem at each time step. A custom loss function ensures the simulation adheres to physical laws like incompressibility (zero-divergence) and, crucially, preserves the swirling motion (vorticity) that gives fluids their character.

Each point represents a Gaussian splat, with its color mapped to the logarithm of its anisotropy ratio. This visualization highlights how individual Gaussians are stretched or compressed in different regions of the domain.

Enforcing Physics Through Optimization: The Loss Functions

The core of the “Gaussian Fluids” method is its departure from traditional time-stepping schemes. Instead of discretizing the Navier-Stokes equations directly, the authors reframe the temporal evolution of the fluid as an optimization problem. At each time step, the solver seeks to find the set of Gaussian parameters that best satisfies the governing physical laws. This is achieved by minimizing a composite loss function, which is a weighted sum of several terms, each designed to enforce a specific physical principle or ensure numerical stability.

Two of the main loss terms are:

Vorticity Loss ($L_{vor}$): For a two dimensional ideal (inviscid) fluid, vorticity, defined as the curl of the velocity field ($\omega=\nabla \times u$), is a conserved quantity that is transported with the flow. This loss function is designed to uphold this principle. It quantifies the error between the vorticity of the current velocity field and the target vorticity field, which has been advected from the previous time step. By minimizing this term, the solver actively preserves the rotational and swirling structures within the fluid. This is critical for preventing the numerical dissipation (an artificial smoothing of details) that often plagues grid-based methods and for capturing the fine, filament-like features of turbulent flows.
Divergence Loss ($L_{div}$): This term enforces the incompressibility constraint, $\nabla \cdot u=0$. From a physical standpoint, this condition ensures that the fluid’s density remains constant; it cannot be locally created, destroyed, or compressed. The loss function achieves this by evaluating the divergence of the Gaussian-represented velocity field at numerous sample points and penalizing any non-zero values. Minimizing this term is essential for ensuring the conservation of mass and producing realistic fluid motion.

Beyond these two, to function as a practical solver, the system must also handle interactions with boundaries and maintain stability over long simulations. The following loss terms address these requirements:

Boundary Losses ($L_{b1}$, $L_{b2}$): The behavior of a fluid is critically dependent on its interaction with solid boundaries. These loss terms enforce such conditions. By sampling points on the surfaces of geometries within the scene, the loss function penalizes any fluid velocity that violates the specified boundary conditions. This can include “no-slip” conditions (where the fluid velocity matches the surface velocity, typically zero) or “free-slip” conditions (where the fluid can move tangentially to the surface but not penetrate it). This sampling-based strategy provides an elegant way to handle complex geometries without the need for intricate grid-cutting or meshing procedures.
Regularization Losses ($L_{pos}$, $L_{aniso}$, $L_{vol}$): These terms are included to ensure the numerical stability and well-posedness of the optimization problem.
- The Position Penalty ($L_{pos}$) constrains the centers of the Gaussians, preventing them from deviating excessively from the positions predicted by the initial advection step. This regularizes the optimization process, improves temporal coherence, and helps prevent particle clustering.
- The Anisotropy ($L_{aniso}$) and Volume ($L_{vol}$) losses act as geometric constraints on the Gaussian basis functions themselves. They penalize Gaussians that become overly elongated, stretched, or distorted. This is crucial for maintaining a high-quality spatial representation and preventing the numerical instabilities that could arise from ill-conditioned Gaussian kernels.

However, when formulating a fluid simulation as a Physics-Based optimization problem, one challenge that soon becomes apparent is that choosing good loss functions can be really tricky. The reliance on “soft constraints”, such as the divergence loss function, in the optimization process means that small errors can accumulate over time. Also, the handling of complex domains with holes (non-simply connected domains) is also a big challenge.
When running the Leapfrog simulation we can see that the divergence of the solution, although relatively small, isn’t as close to zero as it should be (ideally it should be exactly zero).

This “Leapfrog” simulation highlights the challenge of “soft constraints.” The “Divergence” plot (bottom-left) shows the small, non-zero errors that can accumulate, while the “Vorticity” plot (bottom-right) shows the swirling structures the solver aims to preserve.

Divergence-free Gaussian Representation

Now we want to extend on this paper a bit. Let’s focus on the 2D case for now. The problem of regularizing divergence is that the velocity vector field won’t be exactly divergence free. I can perhaps help this with an extra step every iteration. Define $\mathbf{J} \in \mathbb{R}^{2 \times 2}$ as a rotation matrix

$$
\mathbf{J} = \begin{bmatrix}
0 & -1\\
1 & 0
\end{bmatrix}.
$$

It is known that for any potential function $\phi: \mathbb{R}^{2} \to \mathbb{R}$, we have $\mathbf{J} \ \nabla\phi(\mathbf{x})$ divergence free:

$$
\nabla \cdot (\mathbf{J} \ \nabla \phi) = –
\frac{\partial }{\partial x} \frac{\partial \phi}{ \partial y} + \frac{\partial}{\partial y} \frac{\partial \phi}{\partial x} = 0.
$$

This potential function can be constructed with the same gaussian mixture by carrying a $\phi_{i} \in \mathbb{R}$ at every gaussian, and thus a continuous function having full support

$$
\phi(\mathbf{x}) = \sum_{i=1}^{N} \phi_i G_i(\mathbf{x}).
$$

similarly we can replace $G_i(\mathbf{x})$ with the clamped gaussians $\tilde{G}_i(\mathbf{x})$ for efficiency. Now we can construct a vector field that is unconditionally divergence-free

$$
\mathbf{u}_{\phi}(\mathbf{x}) = \mathbf{J} \ \nabla \phi(\mathbf{x})
$$

and the gradient of the potential can be computed

$$
\begin{align}
\nabla \phi(\mathbf{x}) &= \sum_{i=1}^{N} \phi_i \nabla G_i(\mathbf{x}) \\
&= – \sum_{i=1}^{N} \phi_i G_i(\mathbf{x})
\mathbf{\Sigma}_{i}^{-1}
(\mathbf{x} – \mathbf{\mu}_i).
\end{align}
$$

so the equation can be written as

$$
\mathbf{u}_{\phi}(\mathbf{x})
= \sum_{i=1}^{N}
\underbrace{\left(
– \phi_i \mathbf{J} \mathbf{\Sigma}_{i}^{-1}(\mathbf{x} – \mathbf{\mu}_i)
\right)}_{\text{per-gaussian vector}}
G_i(\mathbf{x})
$$

We want to somehow fit this $\mathbf{u}_{\phi}(\mathbf{x})$ to the vector field $\mathbf{u}(x)$ directly from the gaussian velocities

$$
\mathbf{u}(\mathbf{x}) = \sum_{i=1}^{N} \mathbf{u}_i G_i(\mathbf{x}).
$$

This can be a simple loss

$$
\underset{\{\phi_i\}_{i=1}^{N}}{\arg\min} \frac{1}{d|\mathcal{D}|}
\int_{\mathcal{D}}
\|\mathbf{u}(\mathbf{x}) – \mathbf{J} \ \nabla \phi(\mathbf{x}) \|_{2}^{2}
\ \mathrm{d}\mathbf{x}
$$

and we use the similar approach of Monte Carlo way of sampling $\mathbf{x}$ uniformly in space.

Doing this, we managed to achieve a divergence much closer to zero, as shown in the image below.

References

Jingrui Xing, Bin Wang, Mengyu Chu, and Baoquan Chen. 2025. Gaussian Fluids: A Grid-Free Fluid Solver based on Gaussian Spatial Representation. In Proceedings of the Special Interest Group on Computer Graphics and Interactive Techniques Conference Conference Papers (SIGGRAPH Conference Papers ’25). Association for Computing Machinery, New York, NY, USA, Article 9, 1–11. https://doi.org/10.1145/3721238.3730620

Post author By singhishree
Post date September 26, 2025
Tags fluid-simulation, gaussian-splatting

Research

Topology Control: Pathfinding for Genus Preservation (2/2)

Post author By singhishree
Post date August 10, 2025

SGI Fellows:
Stephanie Atherton, Marina Levay, Ualibyek Nurgulan, Shree Singhi, Erendiro Pedro

Recall (from post 1/2) that while DeepSDFs are a powerful shape representation, they are not inherently topology-preserving. By “preserving topology,” we mean maintaining a specific topological invariant throughout shape deformation. In low-dimensional latent spaces $z_{\text{dim}} = 2$, linear interpolation between two topologically similar shapes can easily stray into regions with a different genus, breaking that invariant.

Visualizing genera count in 2D latent space

Higher-dimensional latent spaces may offer more flexibility—potentially containing topology-preserving paths—but these are rarely straight lines. This points to a broader goal: designing interpolation strategies that remain within the manifold of valid shapes. Instead of naive linear paths, we focus on discovering homotopic trajectories—routes that stay on the learned shape manifold while preserving both geometry and topology. In this view, interpolation becomes a trajectory-planning problem: navigating from one latent point to another without crossing into regions of invalid or unstable topology.

Pathfinding in Latent Space

Suppose you are a hiker hoping to a reach a point $A$ from a point $B$ amidst a chaotic mountain range. Now your goal is to plan your journey so that there will be minimal height change in your trajectory, i.e., you are a hiker that hates going up and down much! Fortunately, an oracle gives us a magical device, let us call it $f$, that can give us the exact height of any point we choose. In other words, we can query $f(x)$ for any position $x$, and this device is differentiable, $f'(x)$ exists!

Metaphors aside, the problem of planning a path from a particular latent vector $A$ to another $B$ in the learned latent space would greatly benefit from another auxiliary network that learns the mapping from the latent vectors to the desired topological or geometric features. We will introduce a few ways to use this magical device – a simple neural network.

Gradient as the Compass

Now the best case scenario would be to stay at the same height for the whole of the journey, no? This greedy approach puts a hard constraint on the problem, but it also greatly reduces the possible space of paths to take. For our case, we would love to move toward the direction that does not change our height, and this set of directions precisely forms the nullspace of the gradient.

No height change in mathematical terms means that we look for directions where the derivatives equal zero, as in

$$D_vf(x) = \nabla f(x)\cdot v = 0,$$

where the first equality is a fact of the calculus and the second shows that any desirable direction $v\in \text{null}(\nabla f(x))$.

This does not give any definitive directions for a position $x$ but a set of possible good directions. Then a greedy approach is to take the direction that aligns most with our relative position against the goal – the point $B$.

Almost done, but the problem is that if we always take steps with the assumption that we are on the correct isopath, we would end up accumulating errors. So, we actually want to nudge ourselves a tiny bit along the gradient to land back on the isopath. Toward this, let $x$ be the current position and $\Delta x=\alpha \nabla f(x) + n$, where $n$ is the projection of $B-x$ to $\text{null}(\nabla f(x))$. Then we want $f(x+\Delta x) = f(B)$. Approximate

$$f(B) = f(x+\Delta x) \approx f(x) + \nabla f(x)\cdot \Delta x,$$

and substituting for $\Delta x$, see that

$$f(x) + \nabla f(x) \cdot (\alpha\nabla f(x) + n) \
\Longleftrightarrow
\alpha \approx \frac{f(B) – f(x)}{|\nabla f(x)|^2}.$$

Now we just take $\Delta x = \frac{f(B) – f(x)}{|\nabla f(x)|^2}\nabla f(x) + \text{proj}_{\text{null}(\nabla f(x))} (B- x)$ for each position $x$ on the path.

Optimizing Vertical Laziness

Sometimes it is impossible to hope for the best case scenario. Even when $f(A)=f(B)$, it might happen that the latent space is structured such that $A$ and $B$ are in different components of the level set of the function $f$. Then there is no hope for a smooth hike without ups and downs! But the situation is not completely hopeless if we are fine with taking detours as long as such undesirables are minimized. We frame the problem as

$$\text{argmin}_{x_i \in L, \forall i} \sum_{i=1}^n |f(x_i) – f(B)|^2 + \lambda\sum_{i=1}^{n-2} |(x_{i+2}-x_{i+1}) – (x_{i+1} – x_i)|^2,
$$

where $L$ is the latent space and ${x_i}_{i=1}^n$ is the sequence defining the path such that $x_1 = A$ and $x_n=B$. The first term is to encourage consistency in the function values, while the second term discourages sudden changes in curvature, thereby ensuring smoothness. Once defined, various gradient-based optimization algorithms can be used on this problem, which is now imposed with a relatively soft constraint.

Results of pathfinding via optimization in different contexts.

A Geodesic Path

One alternative attempt of note to optimize pathfinding used the concept of geodesics, or the shortest, smoothest distance between two points on a Riemannian manifold. In the latent space, one is really working over a latent data manifold, so thinking geometrically, we experimented with a geodesic path algorithm. In writing this algorithm, we set our two latent endpoints and optimized the points in between by minimizing the energy of our path on the output manifold. This alternative method worked similarly well to the optimized pathfinding algorithm posed above!

Conclusion and Future Work

Our journey into DeepSDFs, gradient methods, and geodesic paths has shown both the promise and the pitfalls of blending multiple specialized components. The decoupled design gave us flexibility, but also revealed a fragility—every part must share the same initialized latent space, and retraining one element (like the DeepSDF) forces a retraining of its companions, the regressor and classifier. With each component carrying its own error, the final solution sometimes felt like navigating with a slightly crooked compass.

Yet, these challenges are also opportunities. Regularization strategies such as Lipschitz constraints, eikonal terms, and Laplacian smoothing may help “tame” topological features, while alternative frameworks like diffeomorphic flows and persistent homology could provide more robust topological guarantees. The full integration of these models remains a work in progress, but we hope that this exploration sparks new ideas and gives you, the reader, a sharper sense of how topology can be preserved—and perhaps even mastered—in complex geometric learning systems.

Genus is a commonly known and widely applied shape feature in 3D computer graphics, but there is also a whole mathematical menu of topological invariants to explore along with their preservation techniques! Throughout our research, we also considered:

Would topological accuracy benefit from a non-differentiable algorithm (e.g. RRT) that can directly encode genus as a discrete property? How would we go about certifying such an algorithm?
How can homotopy continuation be used for latent space pathfinding? How will this method complement continuous deformations of homotopic and even non-homotopic shapes?
What are some use cases to preserve topology, and what choice of topological invariant should we pair with those cases?

We invite all readers and researchers to take into account our struggles and successes in considering the questions posed.

Acknowledgements. We thank Professor Kry for his mentorship, Daria Nogina and Yuanyuan Tao for many helpful huddles, Professor Solomon for his organization of great mentors and projects, and the sponsors of SGI for making this summer possible.

References

DeepSDF

SDFConnect

Post author By singhishree
Post date August 10, 2025
Tags Deep Learning, Topology

Tutorial week

Ironing Out Wrinkles with Mesh Fairing

Post author By singhishree
Post date August 5, 2025

Picture your triangle mesh as that poor, crumpled T-shirt you just pulled out of the dryer. It’s full of sharp folds, uneven patches, and you wouldn’t be caught dead wearing it in public—or running a physics simulation! In geometry processing, these “wrinkles” manifest as skinny triangles, janky angles, and noisy curvature that spoil rendering, break solvers, or just look… ugly. Sometimes, we might want to fix these problems using something geometry nerds call surface fairing.

Figure 1. (Left) Clean mesh, (Right) Mesh damaged by adding Gaussian noise to each vertex (Middle) Result after applying surface fairing to the damaged mesh.

In essence, surface fairing is an algorithm one uses to smoothen their mesh. Surface fairing can be solved using energy minimization. We define an “energy” that measures how wrinkly our mesh is, then use optimization algorithms to slightly nudge each vertex so that the total energy drops.

Spring Energy
Think of each mesh edge (i,j) as a spring with rest length 0. If a spring is too long, it costs energy. Formally:
E_mem = ½ Σ_(i,j)∈E w_ij‖x_i − x_j‖² with weights w_ij uniform or derived from entries of the discrete Laplacian matrix (using cotangents of opposite angles). We’ll call this term the Membrane Term
Bending Energy
Sometimes, we might want to additionally smooth the mesh even if it looks less wrinkly. In this case, we penalize the curvature of the mesh, i.e., the rate at which normals change at the vertices:
E_bend = ½ Σ_i A_i‖Δx_i‖² where Δx_i is the discrete Laplacian at vertex i and A_i is its “vertex area”. For the uninitiated, the idea of a vertex having an area might seem a little weird. But it represents how much area a vertex controls and is often calculated as one-third of the summed areas of triangles incident to vertex i. We’ll call this energy term the Laplacian Term.
Angle Energy
As Nicholas has warned the fellows, sometimes, having long skinny triangles can cause numerical instability. So we might like our triangles to be more equilateral. We can additionally punish deviations from 60° for each triangle, by adding another energy term: E_angle = Σ_k=1..3(θ_k−60°)². We’ll call this the Angle Term.

Figure 2. (Left) Damaged mesh. (In order) Surface fairing using only 1) Spring Energy 2) Bending Energy 3) Angle Energy. Notice how the last two energy minimizations fail to smooth the mesh properly alone.

Note that in most cases using angle energy and bending energy terms are optional, but using the spring energy term is important! (however, you may encounter a few special cases where you can make do without the spring energy?? i’m not too sure, don’t take my word for it). Once we have computed these energies, we weight them appropriately with scaling factors and sum them up. The next step is to minimize them. I love Machine Learning and hate numerical solvers; and so, it brings me immense joy to inform you that since the problem is non-convex, we ~~can~~ should use gradient descent! At each vertex, compute the gradient ∇_xE and take a tiny step towards the minima: $x_{i} \leftarrow x_i – \lambda \cdot (\nabla_{x} E)_{i}$

Now, have a look at our final surface fairing algorithm in its fullest glory. Isn’t it beautiful? Maybe you and I are nerds after all 🙂

Figure 1. (Left) Clean mesh, (Right) Damaged mesh (Middle) Resulting mesh after each step of gradient descent on the three energy terms combined. Note how the triangles in this case are *more equilateral* than the red mesh in Figure 2, because we’ve combined all energy terms.

So, the next time your mesh looks like you’ve tossed a paper airplane into a blender, remember: with a bit of math and a few iterations, you can make it runway-ready; and your mesh-processing algorithm might just thank you for sparing its life. You can find the code for the blog at github.com/ShreeSinghi/sgi_mesh_fairing_blog.

Happy geometry processing!

References

Desbrun, M., Meyer, M., Schröder, P., & Barr, A. H. (1999). Implicit fairing of irregular meshes using diffusion and curvature flow.
Pinkall, U., & Polthier, K. (1993). Computing discrete minimal surfaces and their conjugates.

Post author By singhishree
Post date August 5, 2025
Tags mesh fairing, optimization

Research

Topology Control: Training a DeepSDF (1/2)

Post author By singhishree
Post date August 5, 2025

SGI Fellows:
Stephanie Atherton, Marina Levay, Ualibyek Nurgulan, Shree Singhi, Erendiro Pedro

In the Topology Control project mentored by Professor Paul Kry and project assistants Daria Nogina and Yuanyuan Tao, we sought to explore preserving topological invariants of meshes within the framework of DeepSDFs. Deep Signed Distance Functions are a neural implicit representation used for shapes in geometry processing, but they don’t come with the promise of respecting topology. After finishing our ML pipeline, we explored various topology-preserving techniques through our simple, initial case of deforming a “donut” (torus) into a mug.

DeepSDFs

Signed Distance Field (SDF) representation of a 3D bunny. The network predicts the signed distance from each spatial point to the surface. Source: (Park et al., 2019).

Signed Distance Functions (SDFs) return the shortest distance from any point in 3D space to the surface of an object. Their sign indicates spatial relation: negative if the point lies inside, positive if outside. The surface itself is defined implicitly as the zero-level set: the locus where $\text{SDF}(x) = 0 $.

In 2019, Park et al. introduced DeepSDF, the first method to learn a continuous SDF directly using a deep neural network (Park et al., 2019). Given a shape-specific latent code $ z \in \mathbb{R}^d $ and a 3D point $ x \in \mathbb{R}^3 $, the network learns a continuous mapping:

$$
f_\theta(z_i, x) \approx \text{SDF}^i(x),
$$

where $ f_\theta $ takes a latent code $ z_i $ and a 3D query point $ x $ and returns an approximate signed distance.

The training set is defined as:

$$
X := {(x, s) : \text{SDF}(x) = s}.
$$

Training minimizes the clamped L1 loss between predicted and true distances:

$$
\mathcal{L}\bigl(f_\theta(x), s\bigr)
= \bigl|\text{clamp}\bigl(f_\theta(x), \delta\bigr) – \text{clamp}(s, \delta)\bigr|
$$

with the clamping function:

$$
\text{clamp}(x, \delta) = \min\bigl(\delta, \max(-\delta, x)\bigr).
$$

Clamping focuses the loss near the surface, where accuracy matters most. The parameter $ \delta $ sets the active range.

This is trained on a dataset of 3D point samples and corresponding signed distances. Each shape in the training set is assigned a unique latent vector $ z_i $, allowing the model to generalize across multiple shapes.

Once trained, the network defines an implicit surface through its decision boundary, precisely where $ f_\theta(z, x) = 0 $. This continuous representation allows smooth shape interpolation, high-resolution reconstruction, and editing directly in latent space.

Training Field Notes

We sampled training data from two meshes, torus.obj and mug.obj using a mix of blue-noise points near the surface and uniform samples within a unit cube. All shapes were volume-normalized to ensure consistent interpolation.

DeepSDF is designed to intentionally overfit. Validation is typically skipped. Effective training depends on a few factors: point sample density, network size, shape complexity, and sufficient epochs.

After training, the implicit surface can be extracted using Marching Cubes or Marching Tetrahedra to obtain a polygonal mesh from the zero-level set.

Training Parameters
SDF Delta	1.0
Latent Mean	0.0
Latent SD	0.01
Loss Function	Clamped L1
Optimizer	Adam
Network Learning Rate	0.001
Latent Learning Rate	0.01
Batch Size	2
Epochs	5000
Max Points per Shape	3000
Network Architecture
Latent Dimension	16
Hidden Layer Size	124
Number of Layers	8
Input Coordinate Dim	3
Dropout	0.0
Point Cloud Sampling
Radius	0.02
Sigma	0.02
Mu	0.0
Number of Gaussians	10
Uniform Samples	5000

For higher shape complexity, increasing the latent dimension or training duration improves reconstruction fidelity.

Latent Space Interpolation

One compelling application is interpolation in latent space. By linearly blending between two shape codes $ z_a $ and $ z_b $, we generate new shapes along the path

$$
z(t) = (1 – t) \cdot z_a + t \cdot z_b,\quad t \in [0,1].
$$

Latent space interpolation between mug and torus.

While DeepSDF enables smooth morphing between shapes, it exposes a core limitation: a lack of topological consistency. Even when the source and target shapes share the same number of genus, interpolated shapes can exhibit unintended holes, handles, or disconnected components. These are not artifacts, they reveal that the model has no built-in notion of topology.

However, this limitation also opens the door for deeper exploration. If neural fields like DeepSDF lack an inherent understanding of topology, how can we guide them toward preserving it? In the next post, we explore a fundamental topological property—genus—and how maintaining it during shape transitions could lead us toward more structurally meaningful interpolations.

References

https://arxiv.org/pdf/1901.05103

Post author By singhishree
Post date August 5, 2025
Tags Deep Learning, Topology

Tutorial week

Tutorial Week Day 5: Debugging and Robustness

Post author By singhishree
Post date July 23, 2025

On July 11th, Nicholas Sharp, the creator of the amazing Polyscope library, gave the SGI Community an insightful talk on debugging and robustness in geometry processing—insights that would later save several fellows hours of head-scratching and mental gymnastics. The talk was broadly organized into five parts, each corresponding to a paradigm where bugs commonly arise:

Representing Numbers
Representing Meshes
Optimization
Simulation and PDEs
Geometric Machine Learning

Representing Numbers

The algorithms developed for dealing with geometry are primarily built to work with clean and pretty real numbers. However, computers deal with floats, and sometimes they behave ugly. Floats are computers’ approximations of real numbers. Sometimes, it may require an infinite amount of storage to correctly represent a single number of arbitrary precision (unless you can get away with $π=3.14$). Each floating number can either be represented using 32 bits (single precision) or 64 bits (double precision).

Figure 1. An example of the float32 representation

In the floating realm, we have quirks like:

$(x+y)+z \neq x+(y+z)$
$ a>0, b>0$ but $a+b=a$

It is important to note that floating-point arithmetic is not random; it is simply misaligned with real arithmetic. That is, the same operation will consistently yield the same result, even if it deviates from the mathematically expected one. Possible alternatives to floats are integers and binary fraction representations; however, they come with their own obvious limitations.

Who’s your best friend? The NaN! The NaN is a special “floating-point” computers spit out when they’re asked to perform invalid operations like

$\frac{0}{0}\rightarrow $ NaN
$\sqrt{-1} \rightarrow$ NaN (not $i$ haha)

Every operation against a NaN results in… NaN. Hence, one small slip-up in the code can result in the entire algorithm being thrown off its course. Well then… why should I love NaNs? Because a screaming alarm is better than a silent error. If your algorithm is dealing with positive numbers and it computes $\sqrt{-42}$ somewhere, you would probably want to know if you went wrong.

Well, what are some good practices to minimize numerical error in your code?

Don’t use equality for comparison, use thresholds like: $\left| x – x^* \right| < \epsilon$ or $\frac{\left| x – x^* \right|}{\left| x^* \right|} < \epsilon$
Avoid transcendental functions wherever possible. A really cool example of this is to avoid $\theta = \arccos \left( \frac{ \mathbf{u} \cdot \mathbf{v} }{ \left| \mathbf{u} \right| \left| \mathbf{v} \right| } \right)$ and use $\cot \theta = \frac{ \cos \theta }{ \sin \theta } = \frac{ \mathbf{u} \cdot \mathbf{v} }{ \left| \mathbf{u} \times \mathbf{v} \right| }$ instead.
Clamp inputs to safe bounds, e.g. $\sqrt{x} \rightarrow \sqrt{max(x, 0)}$. (Note that while this keeps your code running smoothly, it might convert NaNs into silent errors!)
Perturb inputs to singular functions to ensure numerical stability, e.g. $\frac{1}{\left| x \right|} \to \frac{1}{\left| x \right| + \epsilon}$

Some easy solutions include:

Using a rational arithmetic system, as mentioned before. Caveats include: 1) no transcendental operations 2) some numbers might require very large integers to represent them, leading to performance issues, in terms of memory and/or speed.
Use robust predicates for specific functional implementations, the design of which is aware of the common floating-point problems

Representing Surface Meshes

Sometimes, the issue lies not with the algorithm but with the mesh itself. However, we still need to ensure that our algorithm works seamlessly on all meshes. Common problems (and solutions) include:

Unreferenced vertices and repeated vertices: throw them out
A mixture of quad and triangular meshes: subdivide, retriangulate, or delete
Degenerate faces and spurious topology: either repair these corner cases or adjust the algorithm to handle these situations
Non-manifold and non-orientable meshes: split the mesh into multiple manifolds or orientable patches
Foldover faces, poorly tessellated meshes, and disconnected components: use remeshing algorithms like Instant Meshes or use Generalized Winding Numbers instead of meshes

Optimization

Several geometry processing algorithms involve optimization of an objective function. Generally speaking, linear and sparse linear solvers are well-behaved, whereas more advanced methods like gradient descent or non-linear solvers may fail. A few good practices include:

Performing sanity checks at each stage of your code, e.g. before applying an algorithm that expects SPD matrices, check if the matrix you’re passing is actually SPD
When working with gradient-descent-like solvers, check if the gradient magnitudes are too large; that may cause instability for convergence

Simulations and PDEs

Generally, input meshes and algorithms are co-designed for engineering and scientific computing applications, so there aren’t too many problems with their simulations. However, visual computing algorithms need to be robust to arbitrary inputs, as they are typically just one component of a larger pipeline. Algorithms often fail when presented with “bad” meshes, even if it is perfect in terms of connectivity (like being a manifold, being oriented, etc.). Well then, what qualifies as a bad mesh? Meshes with skinny or obtuse triangles are particularly problematic. The solution is to remesh them using more equilateral triangles

Figure 2. (Left) An example of a “bad” mesh with skinny triangles. (Middle) A high-fidelity re-mesh that might be super-expensive to process. (Right) A low-fidelity re-mesh that trades fidelity for efficiency.

Geometric Machine Learning

Most geometric machine learning stands very directly atop geometry processing; hence, it’s important to get it right. The most common problems encountered in geometric machine learning are not so different from those encountered in standard machine learning. These problems include:

Array shape errors: You can use the Python aargh library to monitor and validate tensor shapes
NaNs and infs: Maybe your learning rate is too big? Maybe you’re passing bad inputs into singular functions? Use torch.autograd.set_detect_anomaly(mode, check_nan=True) to track these problematic numbers at inception.
“My trained model works on shape A, but not shape B.”: Is your normalization, orientation, and resolution consistent?

It is a good idea to overfit your ML model on a single shape and ensure that it works on simple objects like cubes and spheres before moving on to more complex examples.

And above all, the golden rule when your algorithm fails:
Visualize everything.

Post author By singhishree
Post date July 23, 2025

Standard (post-log)	Continuous log-PDE	Proposed discrete-log
1. Start with the heat PDE for \(u\).	1. Start with the log-heat PDE for \(w=\log u\).	1. Start with the heat PDE for \(u\).
2. Discretize in space \(\Rightarrow\) ODE for \(u\).	2. Discretize in space for \(\nabla\) and \(\Delta\) \(\Rightarrow\) ODE for \(w\).	2. Discretize in space \(\Rightarrow\) ODE for \(u\).
3. Time-integrate in \(u\)-space.	3. Time-integrate in \(w\)-space and output \(w\).	3. Apply change of variables \(w=\log u\) \(\Rightarrow\) ODE for \(w\).
4. Compute \(w=\log u\) as a post-process.		4. Time-integrate in \(w\)-space and output \(w\).