Tutorial 9 – The Dirichlet problem

In the lecture we saw that the Dirichlet energy has a unique minimizer among all functions with prescribed boundary values. In this tutorial we want to visualize these minimizers in the discrete setting.

Let MR2 be a triangulated surface with boundary M and let V, E  and F denote the set of vertices, edges and triangles of the underlying simplicial complex Σ. Further, let piR2 denote the position of the vertex iV.

To each discrete function fRV we assigned a corresponding piecewise affine function f^=iVfiϕi. Here ϕi denotes the hat function corresponding to the vertex i, i.e. ϕi is the unique piecewise affine function such that ϕi(pj)=δij.

Further we have computed the Dirichlet energy of an affine function f^ on a single triangle: If TσR2 denotes the triangle corresponding to σ={i,j,k}F and αjki denotes the interior angle in Tσ at the vertex i, thenTσ|gradf^|2=12(cotαjki(fjfk)2+cotαkij(fkfi)2+cotαijk(fifj)2).The Dirichlet energy of a discrete function fRV is then defined as the Dirichlet energy of the corresponding piecewise affine function f^:2ED(f)=σΣTσ|gradf^|2={i,j,k}F12(cotαjki(fjfk)2+cotαkij(fkfi)2+cotαijk(fifj)2)={i,j}E{i,j,k}F12cotαkij(fifj)2.As a quadratic form the Dirichlet energy has a corresponding symmetric bilinear form: Letωij=12{i,j,k}Fcotαijk and define L by Lij={i,kEωik,if i=j,ωij,if {i,j}E0,else,thenED(f)=12fTLf.Note, that on a surface each edge is contained in exactly one or two triangles. Thus the weights ωij are sums of only one or two cotangents as illustrated in the picture below for the case of an inner edge.

cotanweights

Now let V˚ denote the set of interior vertices and set Vbd:=VV˚. To each g:VbdR we define an affine spaceWg:={fRVf|Vbd=g}.Given g we are looking for a critical points of the restriction ED:WgR.  Therefore we compute the derivative of ED: Let fWg and let ft be a curve in Wg through f with X=ddt|t=0ft. Then we getdfED(X)=ddt|t=0ED(ft)=12ddt|t=0(ftTLft)=12fTLX12XTLf=XTLf.Since ftWg  for all t we obtain that XW0, i.e. X|Vbd=0. Conversely each XW0 yields a curve in Wg through f. Thus f is a critical point of ED if and only if the entries of Lf corresponding to interior vertices vanish: For each vertex iV˚ we have {i,j}Eωijfi{i,j}Eωijfj={i,j}Eωij(fifj)=0.Since ED is convex there is always a global minimum and thus there is always a solution. As in the smooth case the minimizer is unique.

To summarize: Given any function g:VbdR, the unique minimizer fWg of the Dirichlet energy is the solution of the following linear system{{i,j}Eωij(fifj)=0,for each iV˚,fi=gi,for each iVbd.

Now we want to build up a network that computes the minimizer for given a function given on the boundary of a surface in R2.

First we need to specify a surface in R2 and a function g on its boundary. One possible way to do this is to use a space curve: First we draw with a curve node a space curve and probably resample it. Then we use the say y-component as the function g on the boundary and project it to a polygon in the plane xz-plane and remesh it. Here the network

boundary_datanetwork

The pointwrangle node contains the following 3 lines of code:

Then we have a triangulated surface in the plane and a graph over its boundary (the curve we started with).

specifyingboundarydata

Now we want to build up the matrix A of the linear system. To solve the system we will use scipy. Though, for performance reasons, we will do as much as we can in VEX, which is much faster than python. In particular we will avoid python loops when filling the matrix.

We will use a compressed sparse row matrix (scipy.sparse.csr_matrix). The constructor needs 3 arrays of equal length – two integer arrays (rows,colums) specifying the position of non-zero entries and a double array that specifies the non-zero entries – and two integers that specify the format of the matrix.

The non-zero entries in A  are exactly the diagonal entries Aii and the off-diagonal entries Aij that belong to oriented edges ij of the surface.  Accessing Houdini point attributes and Houdini primitive attributes in python is fast, but unfortunately that’s not the case for the Houdini vertex attributes. Fortunately, we can store the vertex attribute in vectors sitting on primitives: For a surface the Houdini vertices of the primitives are one-to-one with its halfedges – each vertex is the start vertex of a halfedge in the boundary of the Houdini primitive.

By this observation above we can get a complete array of non-zero value indices as follows: we can store an attribute i@id = @ptnum (@ptnum itself is not accessible from python) on points and two attributes v@start and v@end (encoding row and column resp.) on Houdini primitive that contain for each edge in the primitive (starting from primhedge(geo,@primnum)) the point number of the source point (@start) and destination point (@end) of the current edge in the Houndini primitive.

extractcombinatorics

Here the VEX code of the primitivewrangle node:

The pointwrangle node just contains one line (i@id = @ptnum;).

Before we can finally start to compute the cotangent weights and write the matrix entries to attributes we still need one more information – whether a point lies on the boundary or not. This we can do e.g. as follows:

boundarypoints

The group node has in its edge tab a toggle to detect ‘unshared edges’, i.e. edge which appear in only one triangle. With this we build a group boundary over whose points we then run with a pointwrangle node setting an attribute i@bd = 1 (the other points then have @bd = 0 per default).

Now let us start to compute the matrix entries. We will do this here in several steps.

matrixentries

Here VEX code of the wrangle nodes in their order:

Now that we have set up the values we can build up the matrix in a python node: The attributes can read out of the Houdini geometry (geo) using methods as pointIntAttribValues(“attributename”) or primFloatAttribValues(“attributename”).  The arrays are concatenated by numpy.append and then handed over to the csr_matrix-constructor. To access the matrix in the next node we can store it as cachedUserData in the node (compare the code below).

We still need to prepare the right hand side of the system. Therefore we extend the function g defined on the boundary vertices to the interior vertices by zero (pointwrangle node) and store it as cachedUserData as well (python node).

rhs

Here the corresponding code

Now we stick both python node – the one in which we stored the matrix and the one in which we stored the right hand side – into another one where we finally solve the linear system.

solvesystem

The cached user data can be accessed by its name. The easiest way to solve the system is to use scipy.sparse.linalg.spsolve. After that the solution is stored as a point attribute. Below the corresponding python code:

Though this can be speeded up when rewriting the system as an inhomogeneous symmetric system, which can then be solved by scipy.sparse.linalg.cg. Here the alternative code:

Now we are back in VEX land and we can finally draw the graph of the solution:

harmonicgraph

Homework (due 5/7 July). Build a network which allows to prescribe boundary values and then computes the corresponding minimizer of the Dirichlet energy.

This entry was posted in Tutorial. Bookmark the permalink.