1 Introduction

The Laplace equation and Possion equation are partial differential equations (PDEs) which arise naturally from analyzing heat flow, Maxwell’s equations of electromagnetism, and shape recognition in computer graphics [?]. We studied a new way to solve these equations on overlapping triangular mesh domains, in order to smoothen the computer animation of objects that are represented using those meshes. We designed several variants of a divide-and-conquer solution method, and tested them against each other for accuracy.

2 Background

The Laplace equation models the steady-state distribution of heat on a domain $A \subset \mathbb{R}^2$ with each point $(p_1,p_2)$ on its boundary $\partial A$ held at a fixed temperature $g(p_1,p_2)$. The steady-state temperature distribution $u(x,y) : A \to \mathbb{R}$ must satisfy

where (2.1) is the Laplace equation, $ \Delta = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}$ is the Laplace operator and (2.2) is the boundary condition. (2.1) can be generalized to the Poisson equation.

\begin{equation} \Delta u = h \text{ for constant } h \in \mathbb{R}. (2.3) \end{equation}

The usual method for solving these PDEs numerically, known as finite element discretization, consists in approximating $A$ using a triangle mesh with vertices $v_1,\dotsc,v_n$ and substituting every function $f$ defined on $A$ by the vector of its finite evaluations $f_i=f(v_i)$. From this discrete perspective, solving either of the equations above becomes an optimization problem over a finite number of variables. The solutions to these optimization problems approximate the smooth solutions of the original PDE.

In the case of the Laplace and Poisson equations, this optimization problem amounts to minimizing the discrete Dirichlet energy function, which has the form

\begin{equation}\label{Eq:DirichletEnergy} E(u) = \frac{1}{2}{u^T}_A Qu_A + {c^T}u_A \text{, where } u_A = (2.4) \end{equation}

for some constant vector $c \in \mathbb{R}^n$. The $n \times n$ cotangent matrix $Q$ is ubiquitous in geometry processing and serves as the discrete equivalent of the continuous Laplace operator seen above; its derivation can be found in [?. ]. This energy is to be minimized over all $u_A \in \mathbb{R}^n$, subject to the boundary condition that fixes the entries of $u_A$ that correspond to vertices on $\partial A$. Expanding the matrix products in (2.4) shows that the right-hand-side is quadratic in the entries of $u_A$, therefore standard techniques of quadratic programming that minimize quadratic functions subject to constraints yield the energy minimizer, or the vector $u_A$ that minimizes (2.4). However, the precise method by which these optimization problems are solved will not be a concern throughout the paper.

3 Domain Decomposition

How well the energy minimizers approximate solutions to the smooth PDEs depends on the quality of the mesh. Two important parameters of a mesh that determine its quality are its coarseness, as measured by the average lengths of edges (segments connecting vertices), and its regularity, as measured by the aspect ratios of its triangles. Very regular meshes have triangles that are roughly equilateral, while long and thin triangles make an irregular mesh. As mesh quality improves, the energy minimizers should converge to solutions to the PDEs.

Generating quality meshes that represent complex shapes is an active field of research. In 3D, it is often much easier to decompose complex 3D shapes into the union of simpler domains, and then generate quality meshes for those domains. Thus, we study a divide-and-conquer approach to minimize the Dirichlet energy over the union of many domains by minimizing parts of that energy over each individual domain, while demanding that the minimizers agree wherever the domains overlap, to ensure a compatible combined minimizer. We will refer to this as the consistency requirement (formalized below). We implemented this approach to handle general pairs of overlapping 2D domains $A$ and $B$, with corresponding meshes.

However, functions over $A$ and $B$ cannot be compared because they are represented, respectively, by vectors $u_A$ and $u_B$, whose entries are values on different sets of vertices. In order to formulate that consistency requirement, we linearly interpolate from $u_A$ over the triangles of $A$’s mesh to reconstruct a function $u_A : A \to \mathbb{R}$; this is to say, we take $u_A$ to be the only function which is linear on every triangle of $A$’s mesh while taking the values of $u_A$ at its vertices. Similarly, we reconstruct $u_B : B \to \mathbb{R}$. This reconstruction can be plotted out, as in Figure 1. Our contribution is to analyze the following two formulations of the consistency requirement as the linear constraints of quadratic programs, in the hopes that they can be generalized to 3D.

Overlap sampling constraint (OSC) $u_A$ and $u_B$ must agree at the combined vertices of both meshes.

Dirichlet sampling constraint (DSC) $u_A$ and $u_B$ must agree at the vertices of both meshes that lie on $\partial(A \cap B)$ (a small subset of the overlap sampling constraint).

Figure 1: Energy minimizers for the OSC (left) and the DSC (right), subject to the same boundary condition $g(x,y) = (x+0.3)^2y^2$. Note how the two generate virtually identical solutions, even though the second one is only enforcing equality on a much smaller number of vertices.

4 Experiments and results

From each constraint we derived a quadratic program to solve the Laplace or Poisson equation over $A \cup B$, subject to a boundary condition $g : \partial(A \cup B) \to \mathbb{R}$. Namely, we minimized a weighted combination of $E(u_A)$ and $E(u_B)$, such that $u_A$ and $u_B$ satisfy that constraint, and take values on vertices in $\partial(A \cup B)$ as dictated by $g$. We subjected each type of quadratic program to a convergence test over two annular domains with annular union (see figure 1), and measured the maximum of the distance of the approximations provided by each method to the known analytical solution. We repeated this test for $\Delta u=0$, $\Delta u=1$ and for different levels of mesh regularity. $\Delta u = 1$ was chosen for a simple nonzero right-hand side; using any other nonzero constant would have simply scaled $u$.

We concurrently tested the “ideal” method: combining the meshes of $A$ and $B$ into a combined mesh for $A \cup B$, and treating the minimizer of the Dirichlet energy over that combined mesh as a ground truth (GT) proxy to the analytic solution, when the latter is unknown. This ground truth would be unavailable in 3D where combining meshes is far more difficult, but it serves as a useful benchmark of our quadratic programs in 2D.

Figure 2: A pair of annular domains, in blue and red.

Our results, presented in Figure 3, show the Dirichlet sampling constraint to be clearly superior. While the overlap constraint appears to be accurate when solving $\Delta u=0$ on regular meshes, it fails and does not even converge to the analytical solution for $\Delta u = 1$, or when the meshes become increasingly irregular. As expected, the ground truth deviates the least from the analytic solution and converges to it the fastest.

We performed similar tests on different domain shapes, which can be seen in Figure 4. Since no analytical solution is available for these pairs of meshes, we simply measured the error between our methods and the co-refined ground truth (yellow line in figure 3), safely assuming that this last one will surely converge to the smooth solution of the PDE. The results obtained from these experiments agreed completely with those observed for the pair of annular domains.

Figure 3: Convergence of our methods for annular domains. The top row shows the results for regular meshes, while the bottom one shows the results for irregular meshes. The left column contains results for the equation $\Delta u=0$, while the right one does so for $\Delta u=1$.

Figure 4: Examples of pairs of meshes over which we have tested convergences

5 Conclusion

The overlap sampling constraint performed inconsistently between $\Delta u = 0$ and $\Delta u = 1$, as well as between regular and irregular meshes, although it was surprisingly accurate when solving the Laplace equation on regular meshes.

The Dirichlet sampling constraint converged in every one of our experiments for both the Laplace and Poisson equations, regardless of regularity and shape of the meshes. While converging slower than the ground truth, its rate of convergence appears to be at least linear. It is the best candidate for future research to generalize to 3D, progressing toward our goal of solving PDEs over the union of multiple domains when computing a common mesh is infeasible.

Acknowledgement

The authors would like to thank Prof. Alec Jacobson for his patient and dedicated mentorship, as well as the Fields Institute for bringing our research group together, and supporting our research with a stipend, through the Fields Undergraduate Summer Research Program. The third author would also like to express deep gratitude to the María Cristina Masavéu Peterson Foundation for its financial support of her studies.