PDEs on Meshes

$\newcommand{\dotp}[2]{\langle #1, #2 \rangle}$ $\newcommand{\enscond}[2]{\lbrace #1, #2 \rbrace}$ $\newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} }$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\umax}[1]{\underset{#1}{\max}\;}$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\uargmin}[1]{\underset{#1}{argmin}\;}$ $\newcommand{\norm}[1]{\|#1\|}$ $\newcommand{\abs}[1]{\left|#1\right|}$ $\newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. }$ $\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\diag}[1]{{diag}\left( #1 \right)}$ $\newcommand{\qandq}{\quad\text{and}\quad}$ $\newcommand{\qwhereq}{\quad\text{where}\quad}$ $\newcommand{\qifq}{ \quad \text{if} \quad }$ $\newcommand{\qarrq}{ \quad \Longrightarrow \quad }$ $\newcommand{\ZZ}{\mathbb{Z}}$ $\newcommand{\CC}{\mathbb{C}}$ $\newcommand{\RR}{\mathbb{R}}$ $\newcommand{\EE}{\mathbb{E}}$ $\newcommand{\Zz}{\mathcal{Z}}$ $\newcommand{\Ww}{\mathcal{W}}$ $\newcommand{\Vv}{\mathcal{V}}$ $\newcommand{\Nn}{\mathcal{N}}$ $\newcommand{\NN}{\mathcal{N}}$ $\newcommand{\Hh}{\mathcal{H}}$ $\newcommand{\Bb}{\mathcal{B}}$ $\newcommand{\Ee}{\mathcal{E}}$ $\newcommand{\Cc}{\mathcal{C}}$ $\newcommand{\Gg}{\mathcal{G}}$ $\newcommand{\Ss}{\mathcal{S}}$ $\newcommand{\Pp}{\mathcal{P}}$ $\newcommand{\Ff}{\mathcal{F}}$ $\newcommand{\Xx}{\mathcal{X}}$ $\newcommand{\Mm}{\mathcal{M}}$ $\newcommand{\Ii}{\mathcal{I}}$ $\newcommand{\Dd}{\mathcal{D}}$ $\newcommand{\Ll}{\mathcal{L}}$ $\newcommand{\Tt}{\mathcal{T}}$ $\newcommand{\si}{\sigma}$ $\newcommand{\al}{\alpha}$ $\newcommand{\la}{\lambda}$ $\newcommand{\ga}{\gamma}$ $\newcommand{\Ga}{\Gamma}$ $\newcommand{\La}{\Lambda}$ $\newcommand{\si}{\sigma}$ $\newcommand{\Si}{\Sigma}$ $\newcommand{\be}{\beta}$ $\newcommand{\de}{\delta}$ $\newcommand{\De}{\Delta}$ $\newcommand{\phi}{\varphi}$ $\newcommand{\th}{\theta}$ $\newcommand{\om}{\omega}$ $\newcommand{\Om}{\Omega}$

This tour explores the numerical computation of PDE on meshes. The spacial derivatives are replaced by discrete Laplacian operators on meshes. The time derivatives are computed with explicit time stepping.

In [2]:

Spacial Differential Operator

We first define differential operator using sparse matrices.

Load a mesh of $n$ vertices.

In [3]:
name = 'bunny';
clear options;
options.name = name;
[X0,F] = read_mesh(name);
n = size(X0,2);
options.name = name;

Display it.

In [4]:

Exercise 1

Compute the cotangent Laplacian $W$ of the mesh $$ W_{i,j} = \text{cot}(\al_{i,j}) + \text{cot}(\be_{i,j}) $$ where $\al_{i,j}$ and $\be_{i,j}$ are the two angles oposite the the edge $(i,j)$.

In [5]:
In [6]:
%% Insert your code here.

Compute the symmetric Laplacian matrix $L=D-W$ where $D = \text{diag}_i( \sum_j W_{i,j} )$

In [7]:
d = full( sum(W,1) );
D = spdiags(d(:), 0, n,n);
L = D - W;

Compute normalized operators $\tilde W = D^{-1}W$ and $\tilde L = D^{-1}L$.

In [8]:
iD = spdiags(d(:).^(-1), 0, n,n);
tW = iD * W;
tL = iD * L;

The Heat Equation on a Function

The heat equation is the prototype of parabolic PDE. It produces a smoothing of the initial data that increases with time.

This PDE computes a function $f(x,t)$ of space and time that solves $$ \pd{f}{t} = -\tilde L f \qandq f(\cdot,t=0) = f_0. $$

Load a texture image.

In [9]:
M = load_image('lena',256);

Compute spherical coordinates $(\th,\phi)$.

In [10]:
v = X0 - repmat(mean(X0,2), [1 n]);
theta = acos(v(1,:)./sqrt(sum(v.^2)))/pi;
phi = (atan2(v(2,:),v(3,:))/pi+1)/2;

Interpolate the texture on the mesh to obtain $f_0$.

In [11]:
x = linspace(0,1,size(M,1));
f0 = interp2(x,x,M',theta,phi)';

Display the initial condition $f_0$ as a function on the mesh.

In [12]:
options.face_vertex_color = f0(:);
plot_mesh(X0,F, options);
lighting none;

The PDE is discretized as $$ f^{(\ell+1)} = f^{(\ell)} - \tau \tilde L f^{(\ell)}. $$ Here $f^{(\ell)}$ is intended to approximate $f(\cdot,t)$ at time $t=\ell\tau$.

Final time.

In [13]:
Tmax = 50;

The time step $\tau$ should be smaller than 1 to ensure stability of the scheme.

In [14]:
tau = .4;

Number of iterations of the method.

In [15]:
niter = ceil(Tmax/tau);

Initial solution at time t=0.

In [16]:
f = f0;

One step of discretization.

In [17]:
f = f - tau*tL*f;

Exercise 2

Compute the linear heat diffusion.

In [18]:
In [19]:
%% Insert your code here.

The Heat Equation on a Surface

It is possible to apply the heat flow to each coordinate of the vertices $X$.

Initial solution at time t=0.

In [20]:
X = X0;

We use an explicit discretization in time of the PDE. Here is one iteration. Note that since the matrix |X| stores the coordinates as row vectors, the application of $\tilde L$ is obtained by multiplication on the right by |tL'|.

In [21]:
X = X - tau*X*tL';

Exercise 3

Compute the linear heat diffusion by iterating this gradient descent.

In [22]: