DEC using gptoolbox

Alec Jacobson

April 14, 2015

weblog/

Speaking with the Caltech Discrete Exterior Calculus (DEC) promulgators, it seems that despite using DEC notation to explain their math they advocate for building the final algebraic entities (discrete cotangent Laplacian) directly. That is, rather than by multiplying the core DEC stars and d's.

Nonetheless, when re-implementing these papers I find it useful to follow the texts step-by-step. To do that, I build DEC operators (differentials) using functions in my gptoolbox for MATLAB. I'll use D0 for d1 and S0 for ★0 and so on:

% List of all "half"-edges: 3*#F by 2
allE = [F(:,[2 3]); F(:,[3 1]); F(:,[1 2])];
% Sort each edge and remember if order changed
[sortallE,I] = sort(allE,2);
% 1 if edge matches half-edge orientation, -1 otherwise
flip = 3-2*I(:,1);
% find unique edges. EMAP(i) tells us where to
% find sortallE(i,:) in uE: so that sortallE(i,:) = uE(EMAP(i),:)
[E,~,EMAP] = unique(sortallE,'rows');
% number of edges
ne = size(E,1);
% number of vertices
nv = size(V,1);
% number of faces
nf = size(F,1);
% S0  nv by nv lumped diagonal mass matrix
S0 = massmatrix(V,F);
% S1  ne by ne diagonal matrix of edge dual-primal ratios (cotangents)
C = cotangent(V,F);
S1 = sparse(EMAP,EMAP,-C(:),ne,ne);
% S2  nf by nf diagonal matrix of inverted face areas
S2 = diag(sparse(doublearea(V,F).^-1))/2;
% D0  ne by nv edge to vertex signed incidence matrix
D0 = sparse(repmat(1:ne,2,1)',E,repmat([1 -1],ne,1),ne,nv);
% D1  nf by ne edge to face signed incidence matrix
D1 = sparse(repmat(1:nf,1,3)',EMAP,flip,nf,ne);

Then you can build familiar matrices like the discrete Laplacian:

L = D0' * S1 * D0;

The analogous code for building these in C++ libigl should hopefully be obvious.