<p>Disclaimer: I haven't read most of the papers that the book refers to, so what you will find below is my own understanding based on very little data (i.e. I may have severely misunderstood some part of the book).</p>
<p><strong>Laplace-Beltrami Operator</strong></p>
<p>In a finite element formulation you typically have a stiffness matrix <span class="math-container">$W$</span> and a mass matrix <span class="math-container">$M$</span>. For a 2-manifold one possible discretisation is the cotangent one:
<span class="math-container">$$W_{ij} = 
\begin{cases}
-\frac{1}{2}(\cot\alpha_{ij}+\cot\beta_{ij}) &amp;\text{for } j\in\mathcal{N}(i),\\
\sum_{k\in\mathcal{N}(i)}C_{ik} &amp;\text{for } j=i, \\
0 &amp;\text{for } j\ne i \land j\notin \mathcal{N}(i),
\end{cases}
$$</span>
where <span class="math-container">$\alpha_{ij}$</span> and <span class="math-container">$\beta_{ij}$</span> are the angles opposite to the edge <span class="math-container">$ij$</span>. You can also find the same up to a factor <span class="math-container">$\frac{1}{2}$</span> in the appendix of the book you linked (note that in said book they use instead of <span class="math-container">$W$</span> and <span class="math-container">$M$</span>, <span class="math-container">$M$</span> and <span class="math-container">$D$</span>), or a few pages prior to the one you cite. In the appendix of the book that you have the mass matrix <span class="math-container">$M$</span> is a lumped (diagonal) matrix, where <span class="math-container">$M_{ii}$</span> is the Voronoi area for vertex <span class="math-container">$i$</span>. Then finally they use <span class="math-container">$\Delta_S\approx L = -M^{-1}W$</span> (here <span class="math-container">$\Delta_S$</span> describes the Laplace-Beltrami operator on the manifold <span class="math-container">$S$</span>).</p>
<p><strong>Continuous Harmonic Formulation</strong></p>
<p>As far as I can tell the continuous formulation of the fairing approach with a membrane surface would be something like this:
<span class="math-container">\begin{align}
\Delta_S x &amp;= 0, \quad x\in S, \\
x &amp;= p, \quad x\in \partial S,
\end{align}</span>
where <span class="math-container">$S$</span> is the blue part in the images, <span class="math-container">$\partial\Omega$</span> is the boundary (the interface between white and blue), and <span class="math-container">$p$</span> are the (fixed) positions at said boundary. You may note that the above is actually a non-linear equation. This is the case because <span class="math-container">$\Delta_{S}$</span> depends on the surface itself, and <span class="math-container">$x$</span> also depends on the surface itself. Thus I would actually suggest (for the sake of devising a numerical algorithm) to instead consider the time evolution equation:
<span class="math-container">\begin{alignat}{3}
\partial_t x(t) &amp;= \Delta_{S(t)} x(t), &amp;\quad&amp; x\in S(t), &amp;\quad&amp; t\in(0,\infty) \\
x(t) &amp;= p, &amp;\quad&amp; x\in \partial S(t), &amp;\quad&amp; t\in[0,\infty), \\
S(0) &amp;= S_0.
\end{alignat}</span>
This should be equivalent to the first formulation at the steady state <span class="math-container">$\partial_t x = 0$</span>.</p>
<p><strong>Discrete Harmonic Formulation</strong></p>
<p>The first continuous formulation can be discretised just as:
<span class="math-container">$$(C-(I-C)L)X = CP,$$</span>
where <span class="math-container">$\Delta_{S} \approx L = -M^{-1}W\in\mathbb{R}^{n\times n}$</span>, <span class="math-container">$X\in\mathbb{R}^{n\times 3}$</span> are the <span class="math-container">$n$</span> vertices' positions that are inside the blue region and also the ones from the boundary (vertices that are from the white region but have a blue neighbour), <span class="math-container">$P\in\mathbb{R}^{n\times 3}$</span> are the initial positions. Finally <span class="math-container">$C = \operatorname{diag}(c)$</span>, where <span class="math-container">$c \approx 1_{\partial S}$</span>:
<span class="math-container">$$c_i=\begin{cases} 1, \quad \text{$i$ is a boundary vertex}, \\ 0, \quad \text{$i$ is not a boundary vertex}. \end{cases}$$</span>
<span class="math-container">$(C-(I-C)L)$</span> essentially negated <span class="math-container">$L$</span> and replaces its rows corresponding to Dirichlet vertices with <span class="math-container">$e_k$</span> (where <span class="math-container">$k$</span> is the index of the row being replaced). You can fully eliminate the Dirichlet vertices by considering the reduced formulation as in FEM also. Note however that in either case <span class="math-container">$L$</span> is actually a function of <span class="math-container">$X$</span> (since these cotan angles and areas depend on the position of the vertices), this means that the system is non-linear. A standard approach is to try and use a fixed point iteration, i.e. have an initial guess <span class="math-container">$X^0$</span>, then compute <span class="math-container">$L^0$</span> from it, and solve the system <span class="math-container">$(C-(I-C)L^0)X^1 = CP$</span>. Repeating this yields the iteration <span class="math-container">$(C-(I-C)L^k)X^{k+1} = CP$</span>, and you stop when e.g. <span class="math-container">$\|X^{k+1}-X^k\|&lt;\epsilon$</span> in some iteration (note that the superscripts here are iteration indices and not powers). I have no idea how stable this is however. You can of course try other approaches for solving nonlinear systems of equations.</p>
<p>For the time-evolution equation we can discretise the time derivative, e.g. using explicit or semi-implicit Euler. In the explicit variant you have:
<span class="math-container">$$\partial_t x(t) = \Delta_{S(t)} x(t) \implies \frac{x^{k+1}-x^k}{\tau} = \Delta_{S^k} x^k,$$</span>
while in the semi-implicit variant you have:
<span class="math-container">$$\partial_t x(t) = \Delta_{S(t)} x(t) \implies \frac{x^{k+1}-x^k}{\tau} = \Delta_{S^k} x^{k+1}.$$</span>
You can of course also try a fully-implicit variant, but that will result in a nonlinear system of equations that you must solve. Once you discretise also the spatial domain and incorporate the boundary conditions you get:
<span class="math-container">\begin{alignat}{2}
x^{k+1} &amp;= (I-C)(I+\tau L^k)x^k + CP &amp;\quad&amp;&amp; \text{(explicit Euler)}\\
(C+(I-C)(I-\tau L^k))x^{k+1} &amp;= CP &amp;\quad&amp;&amp; \text{(semi-implicit Euler)}.
\end{alignat}</span>
The semi-implicit variant should be stable for all <span class="math-container">$\tau&gt;0$</span> I believe, but larger step sizes may give poor results.</p>
<p><strong>Polyharmonic Formulation</strong></p>
<p>For the polyharmonic formulation you need additional boundary conditions. What I can suggest are Navier conditions:
<span class="math-container">\begin{align}
(-\Delta_S)^q x &amp;= 0, \quad x\in S, \\
(-\Delta)^l x &amp;= p^l, \quad x\in \partial S, \, 0\leq l \leq q-1.
\end{align}</span>
Then this gets discretised as
<span class="math-container">$$\left(C\sum_{l=0}^{q-1}(-L)^{l}+(I-C)(-L)^q\right)X = C\sum_{l=0}^{q-1}P^l.$$</span>
Note that the superscripts here refer to powers, except for the superscript of <span class="math-container">$P$</span> which refers to just different matrix indices.</p>
<p><strong>Implementation Details</strong></p>
<p>When you solve the linear systems you do not have to solve simultaneously for each coordinate, you can solve them one at a time. The system matrix remains the same, only the right-hand side changes. You may also want symmetrize the <span class="math-container">$L$</span> as described in the book, then you could use a conjugate gradient solver since the matrices should be symmetric positive definite (provided you have at least 1 Dirichlet vertex).</p>
<p>When you have to deal with the <span class="math-container">$(-L)^j$</span> powers you can just implement them initially as iterated matrix-vector products. Later on you can compute the explicit form of the matrices in order to reuse cache more.</p>
<p>The choice of the Navier conditions for the polyharmonic case lets you directly use <span class="math-container">$(-L)^l$</span> (as seen in the last discretisation) instead of having to derive special modification for the boundary terms.</p>
<p>Also, please post some screenshots if you make this work. It's nice to know if what I mentioned actually works out and see the result of it.</p>

CC BY-SA 4.0
2023-09-13T22:32:34.467Z
2023-09-13T22:42:07.803Z
2023-09-13T22:42:07.803Z
lightxbulb
<p>In the book <a href="http://staff.ustc.edu.cn/%7Elgliu/Courses/DGP_2014_autumn-winter/References/Book_Polygon%20Mesh%20Processing.pdf" rel="nofollow noreferrer">Polygon Mesh Processing</a>, there is a section dedicated to fairing. In particular page 59 has this lovely figure:</p>
<p><a href="https://i.sstatic.net/nTlAb.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/nTlAb.png" alt="enter image description here" /></a></p>
<p>I am specially interested in the rightmost image. According to the text this is done by fin ding the set of values of x such that of <span class="math-container">$\Delta ^ 3 x = 0$</span></p>
<p>The text claims that all that is needed in these cases is to solve</p>
<p><span class="math-container">$$L^3x = 0$$</span></p>
<p>Where <span class="math-container">$L$</span> is a Laplacian matrix.</p>
<p>What I am not following is how to construct this Laplacian matrix, in particular, the text mentions that it is simpler to fix 2 rings of vertices for the case where you want the middle figure, so I assume you want to fix 3 rings for the rightmost figure.</p>
<p>I was hoping someone could explain to me how you actually feed in the coefficients of the fixed vertices into the system.</p>

CC BY-SA 4.0
2023-09-12T09:24:36.457Z
$L^2$ is just matrix multiplication of $L$ with itself: $L^2=L\cdot L$, similarly for $L^3$. In practice you can also just apply the matrix-vector products one at a time: $L^3v = L(L(Lv))$. Note that the above formulation is actually slightly wrong - the mistake dates back to Kobbelt's paper, there should be a mass matrix involved in there, i.e. $LM^{-1}LM^{-1}L$ instead of $L^3$. I will write a more detailed answer once I have the time.
lightxbulb
CC BY-SA 4.0
2023-09-12T09:41:14.790Z
Btw it's a 2-ring neighbourhood, because iterating two times of a connectivity matrix of a 1-ring neghbourhood yields a 2-ring neighbourhood matrix.
lightxbulb
CC BY-SA 4.0
2023-09-12T19:12:06.793Z
@lightxbulb thank you for the answer. I am still confused. I don't fully understand how to fill in L. In particular, in this case the Laplacian is being applied to the coordinates, so they are interdependent. So I don't fully understand how to linearized it.
Makogan
CC BY-SA 4.0
2023-09-12T21:19:03.653Z
You can use the cotan formula on a mesh to compute $L$: https://en.wikipedia.org/wiki/Discrete_Laplace_operator#Mesh_Laplacians The mass matrix is not in general the identity, you can build a lumped mass matrix by summing the triangle areas around a vertex. As far as $x$ being position and $L(x)$ also being a function of position - this is correct. One way to solve it is to freeze $L(x^k)$ at iteration $k$ and solve for $x^{k+1}$: $L(x^k)x^{k+1} = 0$. I don't know whether something additional needs to be done here for stability. Note that the Dirichlet conditions matter.
lightxbulb
CC BY-SA 4.0
2023-09-12T21:43:48.487Z
@lightxbulb sorry but this makes no sense to me. The blue area in the picture is going to vary widely depending on what one wants the final shape to be. But if you use the standard mesh laplacian, the coefficients of your matrix are set by the original guess of the mesh geometry. That sounds wrong, they should be updated over time, no?
Makogan
CC BY-SA 4.0
2023-09-12T22:57:22.107Z
@lightxbulb I think the best way to formulate my question is as follows, I only know how to construct a Laplaacian matrix for scalar functions, in this case we don't have a scalar function, that is throwing me off completely.
Makogan
<mathematics><computational-geometry><geometry>
Faring a region of a mesh connecting two curves
Makogan
# Faring a region of a mesh connecting two curves

(Asked by: username_0 on 2023-09-12)

In the book [Polygon Mesh Processing](http://staff.ustc.edu.cn/%7Elgliu/Courses/DGP_2014_autumn-winter/References/Book_Polygon%20Mesh%20Processing.pdf), there is a section dedicated to fairing. In particular page 59 has this lovely figure:

[![enter image description here](https://i.sstatic.net/nTlAb.png)](https://i.sstatic.net/nTlAb.png)

I am specially interested in the rightmost image. According to the text this is done by fin ding the set of values of x such that of $\Delta ^ 3 x = 0$

The text claims that all that is needed in these cases is to solve

$$L^3x = 0$$

Where $L$ is a Laplacian matrix.

What I am not following is how to construct this Laplacian matrix, in particular, the text mentions that it is simpler to fix 2 rings of vertices for the case where you want the middle figure, so I assume you want to fix 3 rings for the rightmost figure.

I was hoping someone could explain to me how you actually feed in the coefficients of the fixed vertices into the system.

**Tags:** mathematics, computational-geometry, geometry

> **Comments:**
> **username_1**: $L^2$ is just matrix multiplication of $L$ with itself: $L^2=L\cdot L$, similarly for $L^3$. In practice you can also just apply the matrix-vector products one at a time: $L^3v = L(L(Lv))$. Note that the above formulation is actually slightly wrong - the mistake dates back to Kobbelt's paper, there should be a mass matrix involved in there, i.e. $LM^{-1}LM^{-1}L$ instead of $L^3$. I will write a more detailed answer once I have the time.  
> **username_1**: Btw it's a 2-ring neighbourhood, because iterating two times of a connectivity matrix of a 1-ring neghbourhood yields a 2-ring neighbourhood matrix.  
> **username_0**: @username_1 thank you for the answer. I am still confused. I don't fully understand how to fill in L. In particular, in this case the Laplacian is being applied to the coordinates, so they are interdependent. So I don't fully understand how to linearized it.  
> **username_1**: You can use the cotan formula on a mesh to compute $L$: https://en.wikipedia.org/wiki/Discrete_Laplace_operator#Mesh_Laplacians The mass matrix is not in general the identity, you can build a lumped mass matrix by summing the triangle areas around a vertex. As far as $x$ being position and $L(x)$ also being a function of position - this is correct. One way to solve it is to freeze $L(x^k)$ at iteration $k$ and solve for $x^{k+1}$: $L(x^k)x^{k+1} = 0$. I don't know whether something additional needs to be done here for stability. Note that the Dirichlet conditions matter.  
> **username_0**: @username_1 sorry but this makes no sense to me. The blue area in the picture is going to vary widely depending on what one wants the final shape to be. But if you use the standard mesh laplacian, the coefficients of your matrix are set by the original guess of the mesh geometry. That sounds wrong, they should be updated over time, no?  
> **username_0**: @username_1 I think the best way to formulate my question is as follows, I only know how to construct a Laplaacian matrix for scalar functions, in this case we don't have a scalar function, that is throwing me off completely.  

---

## 1 answers

### Answer by username_1 on 2023-09-13. Score: 1 | Accepted answer

Disclaimer: I haven't read most of the papers that the book refers to, so what you will find below is my own understanding based on very little data (i.e. I may have severely misunderstood some part of the book).

**Laplace-Beltrami Operator**

In a finite element formulation you typically have a stiffness matrix $W$ and a mass matrix $M$. For a 2-manifold one possible discretisation is the cotangent one:
$$W\_{ij} =
\begin{cases}
-\frac{1}{2}(\cot\alpha\_{ij}+\cot\beta\_{ij}) &\text{for } j\in\mathcal{N}(i),\\
\sum\_{k\in\mathcal{N}(i)}C\_{ik} &\text{for } j=i, \\
0 &\text{for } j\ne i \land j\notin \mathcal{N}(i),
\end{cases}
$$
where $\alpha\_{ij}$ and $\beta\_{ij}$ are the angles opposite to the edge $ij$. You can also find the same up to a factor $\frac{1}{2}$ in the appendix of the book you linked (note that in said book they use instead of $W$ and $M$, $M$ and $D$), or a few pages prior to the one you cite. In the appendix of the book that you have the mass matrix $M$ is a lumped (diagonal) matrix, where $M\_{ii}$ is the Voronoi area for vertex $i$. Then finally they use $\Delta\_S\approx L = -M^{-1}W$ (here $\Delta\_S$ describes the Laplace-Beltrami operator on the manifold $S$).

**Continuous Harmonic Formulation**

As far as I can tell the continuous formulation of the fairing approach with a membrane surface would be something like this:
\begin{align}
\Delta\_S x &= 0, \quad x\in S, \\
x &= p, \quad x\in \partial S,
\end{align}
where $S$ is the blue part in the images, $\partial\Omega$ is the boundary (the interface between white and blue), and $p$ are the (fixed) positions at said boundary. You may note that the above is actually a non-linear equation. This is the case because $\Delta\_{S}$ depends on the surface itself, and $x$ also depends on the surface itself. Thus I would actually suggest (for the sake of devising a numerical algorithm) to instead consider the time evolution equation:
\begin{alignat}{3}
\partial\_t x(t) &= \Delta\_{S(t)} x(t), &\quad& x\in S(t), &\quad& t\in(0,\infty) \\
x(t) &= p, &\quad& x\in \partial S(t), &\quad& t\in[0,\infty), \\
S(0) &= S\_0.
\end{alignat}
This should be equivalent to the first formulation at the steady state $\partial\_t x = 0$.

**Discrete Harmonic Formulation**

The first continuous formulation can be discretised just as:
$$(C-(I-C)L)X = CP,$$
where $\Delta\_{S} \approx L = -M^{-1}W\in\mathbb{R}^{n\times n}$, $X\in\mathbb{R}^{n\times 3}$ are the $n$ vertices' positions that are inside the blue region and also the ones from the boundary (vertices that are from the white region but have a blue neighbour), $P\in\mathbb{R}^{n\times 3}$ are the initial positions. Finally $C = \operatorname{diag}(c)$, where $c \approx 1\_{\partial S}$:
$$c\_i=\begin{cases} 1, \quad \text{$i$ is a boundary vertex}, \\ 0, \quad \text{$i$ is not a boundary vertex}. \end{cases}$$
$(C-(I-C)L)$ essentially negated $L$ and replaces its rows corresponding to Dirichlet vertices with $e\_k$ (where $k$ is the index of the row being replaced). You can fully eliminate the Dirichlet vertices by considering the reduced formulation as in FEM also. Note however that in either case $L$ is actually a function of $X$ (since these cotan angles and areas depend on the position of the vertices), this means that the system is non-linear. A standard approach is to try and use a fixed point iteration, i.e. have an initial guess $X^0$, then compute $L^0$ from it, and solve the system $(C-(I-C)L^0)X^1 = CP$. Repeating this yields the iteration $(C-(I-C)L^k)X^{k+1} = CP$, and you stop when e.g. $\|X^{k+1}-X^k\|<\epsilon$ in some iteration (note that the superscripts here are iteration indices and not powers). I have no idea how stable this is however. You can of course try other approaches for solving nonlinear systems of equations.

For the time-evolution equation we can discretise the time derivative, e.g. using explicit or semi-implicit Euler. In the explicit variant you have:
$$\partial\_t x(t) = \Delta\_{S(t)} x(t) \implies \frac{x^{k+1}-x^k}{\tau} = \Delta\_{S^k} x^k,$$
while in the semi-implicit variant you have:
$$\partial\_t x(t) = \Delta\_{S(t)} x(t) \implies \frac{x^{k+1}-x^k}{\tau} = \Delta\_{S^k} x^{k+1}.$$
You can of course also try a fully-implicit variant, but that will result in a nonlinear system of equations that you must solve. Once you discretise also the spatial domain and incorporate the boundary conditions you get:
\begin{alignat}{2}
x^{k+1} &= (I-C)(I+\tau L^k)x^k + CP &\quad&& \text{(explicit Euler)}\\
(C+(I-C)(I-\tau L^k))x^{k+1} &= CP &\quad&& \text{(semi-implicit Euler)}.
\end{alignat}
The semi-implicit variant should be stable for all $\tau>0$ I believe, but larger step sizes may give poor results.

**Polyharmonic Formulation**

For the polyharmonic formulation you need additional boundary conditions. What I can suggest are Navier conditions:
\begin{align}
(-\Delta\_S)^q x &= 0, \quad x\in S, \\
(-\Delta)^l x &= p^l, \quad x\in \partial S, \, 0\leq l \leq q-1.
\end{align}
Then this gets discretised as
$$\left(C\sum\_{l=0}^{q-1}(-L)^{l}+(I-C)(-L)^q\right)X = C\sum\_{l=0}^{q-1}P^l.$$
Note that the superscripts here refer to powers, except for the superscript of $P$ which refers to just different matrix indices.

**Implementation Details**

When you solve the linear systems you do not have to solve simultaneously for each coordinate, you can solve them one at a time. The system matrix remains the same, only the right-hand side changes. You may also want symmetrize the $L$ as described in the book, then you could use a conjugate gradient solver since the matrices should be symmetric positive definite (provided you have at least 1 Dirichlet vertex).

When you have to deal with the $(-L)^j$ powers you can just implement them initially as iterated matrix-vector products. Later on you can compute the explicit form of the matrices in order to reuse cache more.

The choice of the Navier conditions for the polyharmonic case lets you directly use $(-L)^l$ (as seen in the last discretisation) instead of having to derive special modification for the boundary terms.

Also, please post some screenshots if you make this work. It's nice to know if what I mentioned actually works out and see the result of it.


---

