How would I implement linear FEM 2D Dirichlet energy stiffness matrix?

[alecjacobson]

Here's what I might like to write:

given
V ∈ ℝ^(3×n)
F ∈ ℤ^(3×m)

Matrix computing gradient in a 2D triangle with corners (v₁,v₂,v₃)
G(v₁,v₂,v₃) = [v₃-v₁ v₂-v₃]⁻ᵀ [-1 1 0
                               -1 0 1]
Area of a 2D triangle with corners (v₁,v₂,v₃)
a(v₁,v₂,v₃) = ½ |v₃-v₁ v₂-v₃|
Quadratic form of Dirichlet energy in a single 2D traignle with corners (v₁,v₂,v₃)
K(v₁,v₂,v₃) = ½ a(v₁,v₂,v₃)G(v₁,v₂,v₃)ᵀG(v₁,v₂,v₃)
Sparse quadratic form of Dirichlet energy for a triangle mesh (V,F)
L_(F_*,i),(F_*,i) += K( v_*,(F_1,i) , v_*,(F_2,i) ,  v_*,(F_3,i) )

where
L∈ ℝ^(n×n)

assuming that we implement inline comments, we treat subscripts as marks to avoid excessive backticking, and we figure out using variables as indices, we get a vanity ½, inverse transpose ⁻ᵀ and local functions.

This candidate requires some new notation for += as a scatter into multiple entries of sliced out of the sparse matrix L. I'm using notation similar to our diagonal matrix += definition, but instead of _i,i I have written:

L_I,J +=  ...

where I ∈ ℤ^n_I and J ∈ ℤ^n_J are vectors of indices. This kind of gather-scatter construction is fairly common in FEM. In my candidate, L_(F_*,i),(F_*,i) += recognizes the i and creates a loop over from i in 1..m.

I could not see a clean way to do this using our current (scalar) sparse constructors. As a concrete feature request, we should have block sparse matrix constructors .