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
test/test_laplacian_rectangle.jl for a more complete version of this example. Also see DiscretePDEs.jl for more examples.
using DiscreteExteriorCalculus const DEC = DiscreteExteriorCalculus using LinearAlgebra: eigen using AdmittanceModels: sparse_nullbasis using Base.Iterators: product using PlotlyJS: plot, heatmap
Create 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) 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 for c in comp.cells] 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] end
Plot the first eigenvector.
Plot the second eigenvector.
Plot the third eigenvector.
Plot the fourth eigenvector.
And so on.