constraint that $|u(x)| = 1$ pointwise, rather than directly solving the stationary problem with the same constraint. This key difference allows one to avoid the necessity of using a non-linear solver because the pointwise constraint can be enforced simply by normalizing the solution in between discrete time steps (see [@ruuth_diffusion-generated_2001; @laux_analysis_2019]). We proceed by defining a discrete Laplacian operator, $\Delta_Q$ corresponding to the connection $\nabla^Q$. $$\label{eq:discrete-laplacian}
  \Delta_Q(u)|_i = \frac{1}{|\mathcal{N}(n_i)|} \sum_{n_j \in \mathcal{N}(n_i)}{\left(u_j - Q_{ij}(u_i)\right)}$$

where $\mathcal{N}(n_i)$ is the one-ring neighborhood of $n_i$ and $|\mathcal{N}(n_i)|$ is the area of that neighborhood. This discrete Laplacian operator allows us to assemble a discrete diffusion equation using a backward Euler time discretization; $$(I - \tau A)u = u_0 + \tau b$$

where $A$ is the matrix form of $\Delta_Q$. We iteratively solve this equation for a small time step $\tau$ and pointwise renormalize the solution between each iteration. The algorithm is described in detail in algorithm \[alg:MBO\].

Let $u^0$ be the solution to $Au = b$. Fix $\tau$, $\delta$, and set $k = 0$. Solve the discrete diffusion equation, $$(I - \tau A)u^{k+1} = u^k + \tau b$$ Set $u^{k+1}_j = \frac{u^{k+1}_j}{|u^{k+1}_j|}$\
$k++$

Generation of a Quad Layout with T-junctions {#sec:quad-layout}
============================================

After designing 