I've encapsulated the code of the mysterious user21 into a helmholzSolve command. The code is at the end of this post. It adds very little to user21's code but it does allow us to examine multiple examples quite easily, though it has certainly not been tested extensively and could be improved quite a lot I'm sure. It should be called as follows:
{ev,if,mesh} = helmholzSolve[g_Graphics, n_Integer, opts:OptionsPattern[]];
In this code, g can be a Graphics object, an ImplicitRegion, or a ParametricRegion defining the region in question, n is an integer determining the number of eigenvalues that will be computed, and opts is a list of options to be passed to the discretization functions. It returns ev a list of the computed eigenvalues, if a list of corresponding eigenfunctions represented as InterpolatingFunctions and the mesh for plotting purposes. Using this, we can compute the eigenfunctions of the unit disk is as easy as follows:
{ev, if, mesh} = helmholzSolve[Disk[], 6];
ev
(* Out: {6.80538, 15.7385, 15.7385, 27.477, 27.477, 31.5901} *)
We can visualize the eigenfunctions as follows:
GraphicsGrid[Partition[Table[ContourPlot[if[[k]][x, y], Element[{x, y}, mesh],
PlotRange -> All, PlotPoints -> 50], {k, 1, 6}], 3]]
Here's a semi-interesting region:
n = 20;
vertices = Table[(1 + (-1)^k/5) {Cos[2 Pi*k/n], Sin[2 Pi*k/n]}, {k, 1, n}];
g = Graphics[{EdgeForm[Black], Gray, Polygon[vertices]}]
And the plot of an eigenfunction:
{ev, if, mesh} = helmholzSolve[g, 6, "MaxCellMeasure" -> 0.005];
Plot3D[-if[[6]][x, y], Element[{x, y}, mesh],
PlotRange -> All, PlotPoints -> 20, Mesh -> All,
MeshStyle -> Opacity[0.3]]
Here's an implicitly defined region with a hole:
{ev, if, mesh} = helmholzSolve[
ImplicitRegion[1/4 < x^2 + y^2 && x^4 + y^6 <= 1, {x, y}],
4];
ContourPlot[if[[4]][x, y], Element[{x, y}, mesh],
PlotRange -> All, PlotPoints -> 40]
Finally, here's the definition of helmholzSolve.
Needs["NDSolve`FEM`"];
helmholzSolve[g_, numEigenToCompute_Integer,
opts : OptionsPattern[]] := Module[
{u, x, y, t, pde, dirichletCondition, mesh, boundaryMesh,
nr, state, femdata, initBCs, methodData, initCoeffs, vd, sd,
discretePDE, discreteBCs, load, stiffness, damping, pos, nDiri,
numEigen, res, eigenValues, eigenVectors, evIF},
(* Discretize the region *)
If[Head[g] === ImplicitRegion || Head[g] === ParametricRegion,
mesh = ToElementMesh[DiscretizeRegion[g], opts],
mesh = ToElementMesh[DiscretizeGraphics[g], opts]
];
boundaryMesh = ToBoundaryMesh[mesh];
(* Set up the PDE and boundary condition *)
pde = D[u[t,x,y], t] - Laplacian[u[t,x,y], {x, y}] + u[t,x,y] == 0;
dirichletCondition = DirichletCondition[u[t,x,y] == 0, True];
(* Pre-process the equations to obtain the FiniteElementData in StateData *)
nr = ToNumericalRegion[mesh];
{state} = NDSolve`ProcessEquations[{pde, dirichletCondition,
u[0, x, y] == 0}, u, {t, 0, 1}, Element[{x, y}, nr]];
femdata = state["FiniteElementData"];
initBCs = femdata["BoundaryConditionData"];
methodData = femdata["FEMMethodData"];
initCoeffs = femdata["PDECoefficientData"];
(* Set up the solution *)
vd = methodData["VariableData"];
sd = NDSolve`SolutionData[{"Space" -> nr, "Time" -> 0.}];
(* Discretize the PDE and boundary conditions *)
discretePDE = DiscretizePDE[initCoeffs, methodData, sd];
discreteBCs = DiscretizeBoundaryConditions[initBCs, methodData, sd];
(* Extract the relevant matrices and deploy the boundary conditions *)
load = discretePDE["LoadVector"];
stiffness = discretePDE["StiffnessMatrix"];
damping = discretePDE["DampingMatrix"];
DeployBoundaryConditions[{load, stiffness, damping}, discreteBCs];
(* Set the number of eigenvalues ignoring the Dirichlet positions *)
pos = discreteBCs["DirichletMatrix"]["NonzeroPositions"][[All, 2]];
nDiri = Length[pos];
numEigen = numEigenToCompute + nDiri;
(* Solve the eigensystem *)
res = Eigensystem[{stiffness, damping}, -numEigen];
res = Reverse /@ res;
eigenValues = res[[1, nDiri + 1 ;; Abs[numEigen]]];
eigenVectors = res[[2, nDiri + 1 ;; Abs[numEigen]]];
evIF = ElementMeshInterpolation[{mesh}, #] & /@ eigenVectors ;
(* Return the relevant information *)
{eigenValues, evIF, mesh}
]
