<p>Yes, your understanding it correct. The Laplace-Beltrami depends on the current state <span class="math-container">$x$</span> and you have to recompute <span class="math-container">$L$</span> if <span class="math-container">$x$</span> changes. Therefore you cannot write matix of <span class="math-container">$\Delta$</span> without knowing the state <span class="math-container">$x$</span>.</p>

<hr>

<p>To expand on solving <span class="math-container">$\frac{\partial x}{\partial t} = -\Delta x$</span>:</p>

<p>Discretization in space turns <span class="math-container">$\Delta$</span> to <span class="math-container">$L$</span>, discretization in time turns <span class="math-container">$\frac{\partial x}{\partial t}$</span> to <span class="math-container">$\frac{x_{n+1} - x_n}{h}$</span> and decide if we evaluate the right hand side at time <span class="math-container">$t_n$</span> or <span class="math-container">$t_{n+1}$</span></p>

<p><strong>Option 1 = Forward Euler</strong>:
<span class="math-container">$$
x_{n+1} - x_n = h L_n x_n
$$</span>
where <span class="math-container">$L_n = L(x_n)$</span>, <span class="math-container">$L$</span> is considered as a function of the state <span class="math-container">$x$</span>. As you have mentioned this method can be unstable if we take too big step <span class="math-container">$h$</span>. Therefore, we often resort to the second option.</p>

<p><strong>Option 2 = Backward Euler</strong>:
<span class="math-container">$$
x_{n+1} - x_n = h L_{n+1} x_{n+1}
$$</span>
Because <span class="math-container">$L_{n+1}$</span> depends on <span class="math-container">$x_{n+1}$</span> the above equation is non-linear system, we cannot solve it with a simple linear solve.</p>

<p>The common approach is to Taylor expand <span class="math-container">$L$</span> around the point <span class="math-container">$x_n$</span>, general Taylor expansion around point <span class="math-container">$x$</span> is
<span class="math-container">$$
L(x + \Delta x) = L(x) + \nabla L(x) \cdot \Delta x + \dots
$$</span>
therefore for <span class="math-container">$x= x_n$</span> and <span class="math-container">$\Delta x = x_{n+1} - x_n$</span> we get
<span class="math-container">$$
L_{n+1} = L(x_{n+1}) = L(x_n + \Delta x) = L(x_n) + \nabla L(x_n) \cdot \Delta x + \dots
$$</span></p>

<p>The simplest thing to do is to take just the first term of the Taylor expansion, i.e. we replace <span class="math-container">$L_{n+1}$</span> with <span class="math-container">$L_n$</span>. This way, we end up with a linear problem
<span class="math-container">$$
x_{n+1} - x_n = h L_{n} x_{n+1}
$$</span></p>

<p>You can take higher order terms or solve the Backward Euler equation with non-linear solver. However, in your case that is probably unnecessary and the headache of computing the gradient of <span class="math-container">$L$</span> is not worth it(I might be wrong here, I do not have hands on experience with solving this equation). </p>

CC BY-SA 4.0
2018-11-04T16:03:16.867Z
Ok, your answer gave me a better insight. One small observation, in your option 2 you said the equation isn't linear, actually the equation isn't linear even in the first case.
user8469759
CC BY-SA 4.0
2018-11-05T12:11:34.343Z
In the first case, it is linear in $x_{n+1}$ but not in $x_n$.
tom
CC BY-SA 4.0
2018-11-05T12:16:41.577Z
I think the linearity should be tested against the $\left\{x_n\right\}_{n \in \mathbb{N}}$ as a sequence, not against single samples.
user8469759
CC BY-SA 4.0
2018-11-04T08:16:36.350Z
2018-11-04T13:15:15.707Z
2018-11-04T13:15:15.707Z
tom
<p>I'm reading through <a href="https://books.google.co.uk/books?id=8zX-2VRqBAkC&amp;printsec=frontcover&amp;dq=polygon%20mesh%20processing&amp;hl=en&amp;sa=X&amp;ved=0ahUKEwip4tz8t5DeAhUrL8AKHZHQCcsQ6AEIKTAA#v=onepage&amp;q=polygon%20mesh%20processing&amp;f=false" rel="noreferrer">Polygon Mesh Processing</a>, section 4.3 (Surface fairing).</p>

<p>The reference I mentioned defined defines "Surface Fairing" as</p>

<blockquote>
  <p>The goal of surface fairing is to compute shapes that are as smooth as possible... in general fair surface should follow the principle of simplest shapes, the surfaces should be free by any unnecessary details or oscillation".</p>
</blockquote>

<p>It continues...</p>

<blockquote>
  <p>This can be modeled by a suitable energy that penalizes unaesthetic behavior of the surface. A minimization of this fairness energy - subject to user defined constraints - eventually yields to the desired shape.</p>
</blockquote>

<p>The energy defined is function of a surface <span class="math-container">$x : \Omega \to \mathcal{S}$</span> and the form is</p>

<p><span class="math-container">$$
E_M(x) = \int_{\Omega} \sqrt{det(\bf{I})}dudv
$$</span></p>

<p>where <span class="math-container">$\bf{I}$</span> is the first fundamental form, however instead of this the dirichlet energy</p>

<p><span class="math-container">$$
\tilde{E}_M(x) = \int_{\Omega} \lVert x_u \rVert^2 + \lVert x_v \rVert^2 dudv
$$</span></p>

<p>Is the one usually used. Higher order energies are also used, depending on how smooth is the shape we're looking for. Using the Euler-Lagrange equation (in a continuous setting) we end up with the simple PDE</p>

<p><span class="math-container">$$
\Delta x = 0,
$$</span></p>

<p>subject to suitable boundary conditions. Discretazion to meshes is done by means of the Laplace-beltrami operator which I'll define below.</p>

<p>In general it is my understanding that everything boils down to solve a problem of the form</p>

<p><span class="math-container">$$
\Delta^k x = 0
$$</span></p>

<p>where <span class="math-container">$x = x(u,v)$</span> is our surface. If we descretize to a mesh the equation above might be rewritten as</p>

<p><span class="math-container">$$
L^k x = 0
$$</span></p>

<p>Where <span class="math-container">$L^k$</span> is the <span class="math-container">$k^{th}$</span> power of the Laplacian matrix, and <span class="math-container">$L$</span> itself is defined using the Discrete Laplace Beltrami operator, which is something like</p>

<p><span class="math-container">$$
\Delta x = \frac{1}{2A_T} \sum_{j} (cotg (\alpha_{i,j}) + cotg (\beta_{i,j}))(x_j - x_i)
$$</span></p>

<p>However despite the fact it might be very simple to understand in principle at least I'm trying to figure how this would be implemented.</p>

<p>This is because to me the angles depends from the mesh vertices and therefore the system isn't actually linear. The other approach I was thinking is to use a gradient flow (diffusion) approach of the form</p>

<p><span class="math-container">$$
\frac{\partial x}{\partial t} = -\Delta^k x
$$</span></p>

<p>the above, in principle can be discretized as follows</p>

<p><span class="math-container">$$
x(t + h) = x(t) - h L^k x(t)
$$</span></p>

<p>And <span class="math-container">$L$</span> might be computed on the current known state of the surface <span class="math-container">$x(t)$</span>, however the above is an explicit integrator which might lead to instability, and solving the issue would be done by solving</p>

<p><span class="math-container">$$
(id + h L^k) x(t + h) = x(t)
$$</span></p>

<p>But in such case there's the problem of <span class="math-container">$L$</span> which would depend from the unknown state <span class="math-container">$x(t+h)$</span>, so the problem I'm having here is actually similar to the first one I exposed...</p>

<p>So the actual question is:</p>

<p>If the mesh <span class="math-container">$x$</span> is unknown so is the discretization of <span class="math-container">$\Delta^k$</span>, how is this operator constructed then?</p>

<p><strong>Update</strong></p>

<p>I believe I've found the <a href="http://multires.caltech.edu/pubs/ImplicitFairing.pdf" rel="noreferrer">original paper</a>, and I quote the bit I find interesting for my question:</p>

<p>From section 5.4.:</p>

<blockquote>
  <p>We can however proceed in exactly the same way, as the changes induced in a time step will be small. We simply compute the non-zero coefficients of the matrix:
  <span class="math-container">$I - \lambda dt K$</span>
  where <span class="math-container">$K$</span> represents the matrix of the curvature normals...</p>
</blockquote>

<p>So I guess the idea is actually what has been suggested already in one of the comments, we compute the matrix <span class="math-container">$I - \lambda dt K$</span> for the current iteration (say <span class="math-container">$i$</span>) and we keep the very same matrix for a small number of iterations <span class="math-container">$n$</span> after such number of iterations we update again the matrix.</p>

<p>However in the paper the autors confirm that the Laplace-Beltrami operator isn't linear and therefore the system isn't in general linear.</p>

<p>Could anyone confirm if my understanding is correct?</p>

CC BY-SA 4.0
2018-10-18T17:23:19.000Z
You're more likely to get good answers if you can make your question make sense to people who aren't reading the same book chapter as you.
Dan Hulme
CC BY-SA 4.0
2018-10-18T17:36:32.147Z
If you point out what I should expound that would help me to formulate better the question.
user8469759
CC BY-SA 4.0
2018-10-18T17:41:27.650Z
What is the problem you're trying to solve and what is $\Delta^k$ in this context?
Dan Hulme
CC BY-SA 4.0
2018-10-19T09:26:38.600Z
Hopefully is a bit better now, please point out anything that can make it clearer.
user8469759
CC BY-SA 4.0
2018-10-19T20:34:55.503Z
If you had a reference that does that it would be great.
user8469759
<mesh>
Surface fairing, implementation detail
user8469759
# Surface fairing, implementation detail

(Asked by: username_0 on 2018-10-18)

I'm reading through [Polygon Mesh Processing](https://books.google.co.uk/books?id=8zX-2VRqBAkC&printsec=frontcover&dq=polygon%20mesh%20processing&hl=en&sa=X&ved=0ahUKEwip4tz8t5DeAhUrL8AKHZHQCcsQ6AEIKTAA#v=onepage&q=polygon%20mesh%20processing&f=false), section 4.3 (Surface fairing).

The reference I mentioned defined defines "Surface Fairing" as

> The goal of surface fairing is to compute shapes that are as smooth as possible... in general fair surface should follow the principle of simplest shapes, the surfaces should be free by any unnecessary details or oscillation".

It continues...

> This can be modeled by a suitable energy that penalizes unaesthetic behavior of the surface. A minimization of this fairness energy - subject to user defined constraints - eventually yields to the desired shape.

The energy defined is function of a surface $x : \Omega \to \mathcal{S}$ and the form is

$$
E\_M(x) = \int\_{\Omega} \sqrt{det(\bf{I})}dudv
$$

where $\bf{I}$ is the first fundamental form, however instead of this the dirichlet energy

$$
\tilde{E}\_M(x) = \int\_{\Omega} \lVert x\_u \rVert^2 + \lVert x\_v \rVert^2 dudv
$$

Is the one usually used. Higher order energies are also used, depending on how smooth is the shape we're looking for. Using the Euler-Lagrange equation (in a continuous setting) we end up with the simple PDE

$$
\Delta x = 0,
$$

subject to suitable boundary conditions. Discretazion to meshes is done by means of the Laplace-beltrami operator which I'll define below.

In general it is my understanding that everything boils down to solve a problem of the form

$$
\Delta^k x = 0
$$

where $x = x(u,v)$ is our surface. If we descretize to a mesh the equation above might be rewritten as

$$
L^k x = 0
$$

Where $L^k$ is the $k^{th}$ power of the Laplacian matrix, and $L$ itself is defined using the Discrete Laplace Beltrami operator, which is something like

$$
\Delta x = \frac{1}{2A\_T} \sum\_{j} (cotg (\alpha\_{i,j}) + cotg (\beta\_{i,j}))(x\_j - x\_i)
$$

However despite the fact it might be very simple to understand in principle at least I'm trying to figure how this would be implemented.

This is because to me the angles depends from the mesh vertices and therefore the system isn't actually linear. The other approach I was thinking is to use a gradient flow (diffusion) approach of the form

$$
\frac{\partial x}{\partial t} = -\Delta^k x
$$

the above, in principle can be discretized as follows

$$
x(t + h) = x(t) - h L^k x(t)
$$

And $L$ might be computed on the current known state of the surface $x(t)$, however the above is an explicit integrator which might lead to instability, and solving the issue would be done by solving

$$
(id + h L^k) x(t + h) = x(t)
$$

But in such case there's the problem of $L$ which would depend from the unknown state $x(t+h)$, so the problem I'm having here is actually similar to the first one I exposed...

So the actual question is:

If the mesh $x$ is unknown so is the discretization of $\Delta^k$, how is this operator constructed then?

**Update**

I believe I've found the [original paper](http://multires.caltech.edu/pubs/ImplicitFairing.pdf), and I quote the bit I find interesting for my question:

From section 5.4.:

> We can however proceed in exactly the same way, as the changes induced in a time step will be small. We simply compute the non-zero coefficients of the matrix:
> $I - \lambda dt K$
> where $K$ represents the matrix of the curvature normals...

So I guess the idea is actually what has been suggested already in one of the comments, we compute the matrix $I - \lambda dt K$ for the current iteration (say $i$) and we keep the very same matrix for a small number of iterations $n$ after such number of iterations we update again the matrix.

However in the paper the autors confirm that the Laplace-Beltrami operator isn't linear and therefore the system isn't in general linear.

Could anyone confirm if my understanding is correct?

**Tags:** mesh

> **Comments:**
> **username_1**: You're more likely to get good answers if you can make your question make sense to people who aren't reading the same book chapter as you.  
> **username_0**: If you point out what I should expound that would help me to formulate better the question.  
> **username_1**: What is the problem you're trying to solve and what is $\Delta^k$ in this context?  
> **username_0**: Hopefully is a bit better now, please point out anything that can make it clearer.  
> **username_0**: If you had a reference that does that it would be great.  

---

## 1 answers

### Answer by username_2 on 2018-11-04. Score: 3 | Accepted answer

Yes, your understanding it correct. The Laplace-Beltrami depends on the current state $x$ and you have to recompute $L$ if $x$ changes. Therefore you cannot write matix of $\Delta$ without knowing the state $x$.

---

To expand on solving $\frac{\partial x}{\partial t} = -\Delta x$:

Discretization in space turns $\Delta$ to $L$, discretization in time turns $\frac{\partial x}{\partial t}$ to $\frac{x\_{n+1} - x\_n}{h}$ and decide if we evaluate the right hand side at time $t\_n$ or $t\_{n+1}$

**Option 1 = Forward Euler**:
$$
x\_{n+1} - x\_n = h L\_n x\_n
$$
where $L\_n = L(x\_n)$, $L$ is considered as a function of the state $x$. As you have mentioned this method can be unstable if we take too big step $h$. Therefore, we often resort to the second option.

**Option 2 = Backward Euler**:
$$
x\_{n+1} - x\_n = h L\_{n+1} x\_{n+1}
$$
Because $L\_{n+1}$ depends on $x\_{n+1}$ the above equation is non-linear system, we cannot solve it with a simple linear solve.

The common approach is to Taylor expand $L$ around the point $x\_n$, general Taylor expansion around point $x$ is
$$
L(x + \Delta x) = L(x) + \nabla L(x) \cdot \Delta x + \dots
$$
therefore for $x= x\_n$ and $\Delta x = x\_{n+1} - x\_n$ we get
$$
L\_{n+1} = L(x\_{n+1}) = L(x\_n + \Delta x) = L(x\_n) + \nabla L(x\_n) \cdot \Delta x + \dots
$$

The simplest thing to do is to take just the first term of the Taylor expansion, i.e. we replace $L\_{n+1}$ with $L\_n$. This way, we end up with a linear problem
$$
x\_{n+1} - x\_n = h L\_{n} x\_{n+1}
$$

You can take higher order terms or solve the Backward Euler equation with non-linear solver. However, in your case that is probably unnecessary and the headache of computing the gradient of $L$ is not worth it(I might be wrong here, I do not have hands on experience with solving this equation).

> **Comments:**
> **username_0**: Ok, your answer gave me a better insight. One small observation, in your option 2 you said the equation isn't linear, actually the equation isn't linear even in the first case.  
> **username_2**: In the first case, it is linear in $x_{n+1}$ but not in $x_n$.  
> **username_0**: I think the linearity should be tested against the $\left\{x_n\right\}_{n \in \mathbb{N}}$ as a sequence, not against single samples.  

---

