Eigenfunctions of the Laplacian on an arbitrary mesh

Question

So, I've constructed a mesh over which I'd like to find eigenfunctions of Laplace's equation with a free boundary (a zero Neumann boundary condition along the edge).

Mostly because I figured an electric guitar–shaped Chladni plate would look really cool (but also because I'd like to try to find the resonant frequencies of a lamina this shape).

Normally, to do this sort of eigenvalue problem, I'd do something like this:

ParametricNDSolve[
  {Laplacian[phi[x, y], {x, y}] == -w^2*phi[x, y]}, 
  {x, y} ∈ mesh, 
  w
]

with some boundary conditions in there.

However, while I feel like I understand (well enough) the theory behind what I want to accomplish (at the very least, I can derive the eigenfunctions in coordinate systems where Laplace's equation is solvable by separation of variables), I don't quite know how to go about implementing it in this case.

How do I impose a free boundary for an arbitrary mesh?

This is relatively easy on a rectangle or circle, but on this mesh, I think would have to somehow determine the vector normal to the mesh boundary at each point (so that I can specify $\mathbf{n}\cdot\nabla \phi = 0$). Furthermore, once I've specified the Neumann boundary condition, are there any other boundary conditions I have to keep in mind (I imagine that if I only specify a zero derivative on the edge, NDSolve will give me a constant solution)? Thank you so much in advance for your help, and I promise to post cool plots of irregularly-shaped guitar Chladni plates when I'm done.

Edit: Alright, I guess the question I was asking wasn't the question I should have been asking. To boil it down, we don't really have to worry about the Neumann boundary condition in this case, as it's zero by default. However, Mathematica doesn't actually have built-in eigenvalue analysis for these sorts of problems (thus the necessity of using some low-level code, as utilized to great effect by user21 here).

Answer 1

As promised:

And here are the nodes:

Now to outline the process. [This first part has no Mathematica.] First, I found an image of a guitar using Google's Image Search. I then went into GIMP (although you can use any image editing software, or even draw the image yourself) and used it as a template to create a silhouette (it's okay if the edges are a little rough). Then, I created a parametric plot using this technique.

img = Binarize[Import["/home/michael/Downloads/test.jpg"]~ColorConvert~"Grayscale"~ImageResize~500~Blur~3]~Blur~3;
lines = Cases[Normal@ListContourPlot[Reverse@ImageData[img], Contours -> {0.5}], _Line, -1];

param[x_, m_, t_] := Module[{f, n = Length[x], nf}, 
  f = Chop[Fourier[x]][[;; Ceiling[Length[x]/2]]]; nf = Length[f];
  Total[Rationalize[2 Abs[f]/Sqrt[n] Sin[Pi/2 - Arg[f] + 2. Pi Range[0, nf - 1] t],.01][[;;Min[m, nf]]]]]

tocurve[Line[data_], m_, t_] := param[#, m, t] & /@ Transpose[data]

parplot = ParametricPlot[Evaluate[tocurve[#, 40, t] & /@ lines], {t, 0, 0.998},
 Frame -> True, Axes -> False, PlotPoints -> 3] /. 
 Line[l_List] :> {{Blue, Polygon[l]}, {White, Line[l]}}

I recommend setting PlotPoints to something very low in order to decrease the need for mesh refinement near the edges of your lamina (although this will decrease the resolution of your boundary curve; I did not do this in my original question, which gave me a mesh that was far too complex near the edge). Also, closing the loop (letting t go to 1 in our plot bounds) was for some reason giving me a weird irregularity (jaggedness of the boundary) that was almost invisible but was causing the mesh to become extremely (and stubbornly) fine near that point.

Then, use DiscretizeGraphics[] and ToElementMesh[] to convert this polygon to an element mesh.

g = DiscretizeGraphics[parplot]

<<NDSolve`FEM`
mesh = ToElementMesh[g, MeshQualityGoal -> 0.8, MaxCellMeasure -> 30, "MeshOrder" -> 1]
mesh["Wireframe"]

Using the answer provided here by user21 (check out his answer for a better explanation of what each piece of this code does), we can then find the eigenfunctions of the Helmholtz differential equation (which is the eigenvalue equation for the Laplacian) over our lamina.

pde = D[u[t, x, y], t] - Laplacian[u[t, x, y], {x, y}] + u[t, x, y] == 0;
Γ = DirichletCondition[u[t, x, y] == 0, True];

nr = ToNumericalRegion[mesh];
{state} = NDSolve`ProcessEquations[{pde, Γ, u[0, x, y] == 0}, u, {t, 0, 1}, {x, y} ∈ nr];

femdata = state["FiniteElementData"]
initBCs = femdata["BoundaryConditionData"];
methodData = femdata["FEMMethodData"];
initCoeffs = femdata["PDECoefficientData"];

vd = methodData["VariableData"];
sd = NDSolve`SolutionData[{"Space" -> nr, "Time" -> 0.}];

discretePDE = DiscretizePDE[initCoeffs, methodData, sd];
discreteBCs = DiscretizeBoundaryConditions[initBCs, methodData, sd];

load = discretePDE["LoadVector"];
stiffness = discretePDE["StiffnessMatrix"];
damping = discretePDE["DampingMatrix"];

DeployBoundaryConditions[{load, stiffness, damping}, discreteBCs]

nDiri = First[Dimensions[discreteBCs["DirichletMatrix"]]];
numEigenToCompute = 10;
numEigen = numEigenToCompute + nDiri;

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;

The above plots can then be created with this:

densityplot[i_] := DensityPlot[Evaluate[evIF[[i]][x, y]], {x, y} \[Element] mesh,
  PlotPoints -> 256, PlotRange -> All, ColorFunction -> "TemperatureMap", PlotLabel -> i]
nodeplot[i_] := ContourPlot[Evaluate[evIF[[i]][x, y]] == 0, {x, y} \[Element] mesh, PlotPoints -> 256, PlotRange -> All, PlotLabel -> i]

Show[GraphicsGrid[Table[{densityplot[2i-1], densityplot[2i]}, {i, 1, 5}]], ImageSize -> Full]
Show[GraphicsGrid[Table[{nodeplot[2i-1], densityplot[2i]}, {i, 1, 5}]], ImageSize -> Full]

Thanks so much to user21 for his brilliant code.