DiscreteExteriorCalculus.jl is a package implementing Discrete Exterior Calculus. Data structures for cell complexes, primal, and dual meshes are provided along with implementations of the exterior derivative, hodge star, codifferential, and Laplace-de Rham operators.
Clone the repository from GitHub and install Julia 1.1. No build is required beyond the default Julia compilation.
Example usage: modes of the Laplace-de Rham operator on a rectangle with Dirichlet boundary conditions
See test/test_laplacian_rectangle.jl for a more complete version of this example. Also see DiscretePDEs.jl for more examples.
Import packages.
using DiscreteExteriorCalculus
const DEC = DiscreteExteriorCalculus
using LinearAlgebra: eigen
using AdmittanceModels: sparse_nullbasis
using Base.Iterators: product
using PlotlyJS: plot, heatmapCreate a rectangular grid that is r1×r2, then subdivide each rectangle into two triangles. Collect these into a cell complex.
r1, r2 = .5, .4
num = 40
points, tcomp = DEC.triangulated_lattice([r1, 0], [0, r2], num, num)Orient the cell complex positively and compute the dual mesh.
orient!(tcomp.complex)
m = Metric(2)
mesh = Mesh(tcomp, circumcenter(m))Compute the Laplace-de Rham operator for 0-forms and a sparse nullbasis for the constraint that the 0-form goes to 0 on the boundary.
laplacian = differential_operator(m, mesh, "Δ", 1, true)
comp = tcomp.complex
_, exterior = boundary_components_connected(comp)
constraint = zero_constraint(comp, exterior.cells[1], 1)
nullbasis = sparse_nullbasis(constraint)Compute the eigenvalues and eigenvectors of the restricted Laplace-de Rham operator.
vals, vects = eigen(collect(transpose(nullbasis) * laplacian * nullbasis))
inds = sortperm(vals)
vals, vects = vals[inds], vects[:, inds]Lift the eigenvectors back to the original space and reshape for plotting.
comp_points = [c.points[1] for c in comp.cells[1]]
ordering = [findfirst(isequal(p), comp_points) for p in points]
vs = [collect(transpose(reshape((nullbasis * vects[:,i])[ordering],
num+1, num+1))) for i in 1:size(vects, 2)]
for i in 1:length(vs)
j = argmax(abs.(vs[i]))
vs[i] /= vs[i][j]
endPlot the first eigenvector.
plot(heatmap(z=vs[1]))
Plot the second eigenvector.
plot(heatmap(z=vs[2]))
Plot the third eigenvector.
plot(heatmap(z=vs[3]))
Plot the fourth eigenvector.
plot(heatmap(z=vs[4]))
And so on.