This notebook contains additional material to the Mathematica GuideBooks  (Springer-Verlag, New York, 2004) http://www.MathematicaGuideBooks.org . It uses the stylesheet of the Mathematica GuideBooks, which is embedded.

Symbolics Volume, Chapter 1

Calculate and visualize some of the eigenvalues and eigenfunctions of the Helmholtz equation Δψ(x,y)=-(ω^2)ψ(x,y) in the interior of the two spectral twins bilby and hawk [✶], [✶], [✶], [✶], [✶], [✶]. Here are the shapes of the bilby and hawk. All horizontal and vertical edges have length 2 and all non-90° angles of the boundary are 45° .

Use higher order triangular finite elements (up to order 12) to discretize the problem. Compare the quality of higher order refinements with more triangles in a regular homogeneous refinement.

For a given order o, the function ShapeFunctionList returns a list of the O[o] shape functions, such that the value of the kth shape function is one at the point basePoint[o][[k]] and zero at the other base points.

𝒪[o_] := (o + 1)(o + 2)/2;

basePoints[o_] := Flatten[Table[{i, j}/o, {j, 0, o}, {i, 0, o - j}], 1]

(* calculate all shape functions of order o at once;
   store calculated list of shape functions *)
ShapeFunctionList[o_, {ξ_, η_}] := ShapeFunctionList[o, {ξ, η}] =
Module[{monomials, ℒ, sols},
  monomials[{x_, y_}] = Flatten[Table[x^i y^(j - i), 
                                      {j, 0, o}, {i, 0, j}]];
  (* form and solve equations efficiently *)
  ℒ = LinearSolve[monomials /@ basePoints[o]];
  (* use all 𝒪[o] right-hand side *)
  sols = Table[ℒ @ Table[KroneckerDelta[i, j], {i, 𝒪[o]}], {j, 𝒪[o]}];               
  (* form polynomials *)
  Evaluate[monomials[{ξ, η}].#]& /@ sols]

Calculating the 66 shape functions of order 10 can be done in a fraction of a second.

ShapeFunctionList[10, {ξ, η}]; // Timing

The function 𝒯ℐ integrates the bivariate polynomial 𝓅 over the unit triangle.

(* integrating a bivariate polynomial over the unit triangle *)
𝒯ℐ[𝓅_, {x_, y_}] := Plus @@ 
   (Function[{p, q, c}, c (p! q!) /(p + q + 2)!] @@@ (Flatten /@ 
             Internal`DistributedTermsList[𝓅, {x, y}][[1]]))

And here are the definitions for the mass and stiffness matrices. For uniformity of notations and the arguments, both have the geometry parameters J, A, B , and C as parameters.

MassMatrix[o_, {J_, A_, B_, C_}] := MassMatrix[o, {J, A, B, C}] =
Module[{M = Table[0, {𝒪[o]}, {𝒪[o]}], SL = ShapeFunctionList[o, {ξ, η}]},
       (* calculate matrix elements and fill in symmetrically *)
       Do[m = 𝒯ℐ[Expand[SL[[i]] SL[[j]]], {ξ, η}];
          M[[i, j]] = m; M[[j, i]] = m, {i, 𝒪[o]}, {j, i}]; J M]

StiffnessMatrix[o_, {J_, A_, B_, C_}] := StiffnessMatrix[o, {J, A, B, C}] =
Module[{M = Table[0, {𝒪[o]}, {𝒪[o]}], SL = ShapeFunctionList[o, {ξ, η}]},
       (* calculate matrix elements and fill in symmetrically *)
       Do[m = (A 𝒯ℐ[Expand[D[SL[[i]], ξ] D[SL[[j]], ξ]], {ξ, η}] +
               C 𝒯ℐ[Expand[D[SL[[i]], η] D[SL[[j]], η]], {ξ, η}] +
               (* for the following, the next term vanishes *)
               B 𝒯ℐ[Expand[D[SL[[i]], ξ] D[SL[[j]], η] +
                           D[SL[[i]], η] D[SL[[j]], ξ]], {ξ, η}]);
          M[[i, j]] = m; M[[j, i]] = m, {i, 𝒪[o]}, {j, i}]; M]

Here are some timings for calculating some higher-order mass and stiffness matrices that we will use later. Each requires the calculation of  (o+1) (o+2) ((o^2)+3 o+4)/8 integrations over the unit triangle. Calculating the 3×4186 integrals for the order 12 stiffness matrices takes a few minutes.

Table[{o, Timing[MassMatrix[o, {J, A, B, C}];][[1]],
          Timing[StiffnessMatrix[o, {J, A, B, C}];][[1]]}, {o, 12}]

Next let us deal with the geometrical shapes. Here are the outlines of bilby and hawk.

{triangles[bilby], triangles[hawk]} =
(* triangles expressed through vertices *)
{Evaluate[Map[Slot, {{1, 2, 3}, {3, 4, 1}, {5, 1, 4}, {1, 5, 9}, 
                     {6, 9, 5}, {6, 8, 9}, {6, 7, 8}}, {2}]]&[
     { 1, -1}, { 1, -3}, { 3, -1}, { 3,  1}, { 1,  1}, 
     {-1,  1}, {-1,  3}, {-3,  1}, {-1, -1}],
 Evaluate[Map[Slot, {{2, 9, 1}, {2, 3, 4}, {2, 4, 9}, {5, 9, 4}, 
                     {5, 8, 9}, {5, 6, 8}, {7, 8, 6}}, {2}]]&[
     { 1, -3}, { 1, -1}, { 3, -1}, { 1,  1}, {-1,  1}, 
     {-1,  3}, {-3,  3}, {-3,  1}, {-1, -1}]};

Because later we will need smaller triangles, we define a function triangleMesh that contains the o-times recursively subdivided triangles of bilby and hawk.

subdivideTriangle[{p1_, p2_, p3_}] :=
Module[{p4 = (p2 + p3)/2}, {{p4, p1, p2}, {p4, p3, p1}}]

triangleMesh[bilbyHawk_][o_] := triangleMesh[bilbyHawk][o] =
 Flatten[subdivideTriangle /@ triangleMesh[bilbyHawk][o - 1], 1];
   
triangleMesh[bilbyHawk_][0] = triangles[bilbyHawk];

Also, for later use, we define the bounding line segments of bilby and hawk.

borderLine[bilby] = {{1, -3}, {3, -1}, {3, 1}, {1, 1}, {-1, 1}, 
                     {-1, 3}, {-3, 1}, {-1, -1}, {1, -1}, {1, -3}};
                     
borderLine[hawk] = {{1, -3}, {1, -1}, {3, -1}, {1, 1}, {-1, 1}, 
                    {-1, 3}, {-3, 3}, {-3, 1}, {-1, -1}, {1, -3}};

The functions edgePointQ[bilbyHawk] calculate if a point {x, y} lies on the boundary of bilby or hawk.

edgePointQ[bilby][{x_, y_}] := (* check all edges *)
-1 <= x <= 3 && y == 1 || x == 3 && -1 <= y <= 1 ||
y == x - 4 && -3 <= y <= -1 || x == 1 && -3 <= y <= -1 ||
-1 <= x <= 1 && y == -1 || y == -x - 2 && -1 <= y <= 1 ||
y == x + 4 && 1 <= y <= 3 || x == -1 && 1 <= y <= 3

edgePointQ[hawk][{x_, y_}] := (* check all edges *)
-1 <= x <= 1 && y == 1 || x == -1 && 1 <= y <= 3 ||
-3 <= x <= -1 && y == 3 || x == - 3 && 1 <= y <= 3 ||
y == -x - 2 && -3 <= y <= 1 || x == 1 && -3 <= y <= -1 ||
1 <= x <= 3 && y == -1 || y == -x + 2 && -1 <= y <= 1

Here are triangulated versions of bilby and hawk and their outlines.

(* f0r better visibility, contract a triangle *)
contract[l_, f_] := Function[m, (m + f (# - m))& /@ l][Plus @@ l/3]

Show[GraphicsArray[
Graphics[{(* the border lines segments *)
          {GrayLevel[0], Thickness[0.02], Line[borderLine[#]]}, 
           (* triangles forming bilby and hawk *)
          {GrayLevel[0.6], Polygon[contract[#, 0.96]]}& /@ 
                                            triangleMesh[#][0],
          (* recursively subdivided triangles *)
          {Hue[Random[]], Polygon[contract[#, 0.5]]& /@ #}& /@ 
                             Partition[triangleMesh[#][4], 8]},
         AspectRatio -> Automatic, Frame -> True,
         FrameTicks -> None]& /@ {bilby, hawk}]]

For getting a couple of correct digits for the eigenvalues, we might have to subdivide the triangles. The function meshData returns a list of the finite element vertices obtained after order subdivisions of all triangles of bilby or hawk for a given finite element order femOrder. The second argument DN indicates the use of Dirichlet boundary conditions (D) or Neumann boundary conditions (N). For Dirichlet boundary conditions, all points at the boundary are eliminated because they have the value zero by definition and their values do not have to be calculated, and for Neumann boundary conditions all points are returned (Neumann boundary conditions are the most natural boundary boundary conditions for finite element discretizations). The numbered finite element vertices are returned in lists, each list containing the vertices of the triangle. Each vertex is specified in the form P[coordinates, {triangleNumber, vertexNumberInsideTriangle}, globalvertexNumber].

(* map the point {x, y} into the triangle {p1, p2, p3} *)
toTriangle[{p1_, p2_, p3_}, {x_, y_}] := p1 + (p2 - p1) x + (p3 - p1) y

meshData[{bilbyHawk_, order_}, DN_, femOrder_] :=
Module[{(* the triangulated bilby and hawk *)
        tm = triangleMesh[bilbyHawk][order], λ = 𝒪[femOrder],
        trianglePoints, relevantPointsInTriangles, numberedPoints, 
        vertexNumber}, 
 (* points in the form {coordinate, {triangleNumber, pointNumber}};
    grouped by triangles *) 
 trianglePoints = Table[Table[P[toTriangle[tm[[j]], 
                                basePoints[femOrder][[k]]], {j, k}], 
                              {k, λ}], {j, Length[tm]}]; 
 (* all shape functions base points in all triangles *)
 relevantPointsInTriangles =
 (* keep only inner points for Dirichlet boundary conditions *)
 If[DN === D, DeleteCases[trianglePoints, 
                          _?(edgePointQ[bilbyHawk][#[[1]]]&), {2}],
    (* use all points for Neumann case *) trianglePoints];                       
 (* the numbered inner points, no multiplicity *)
 numberedPoints = MapIndexed[{#2[[1]], #1}&, Union[First /@ 
                             Flatten[relevantPointsInTriangles]]];
 (* constant-time look-up for vertex numbers *)
 (vertexNumber[#2] = #1)& @@@ numberedPoints;        
 Map[Append[#, vertexNumber[#[[1]]]]&, relevantPointsInTriangles, {2}]]

Here are the numbered vertices for the simplest nontrivial case—the seventh order zero bilby triangles for Dirichlet boundary conditions with degree one basis functions. The triangles at the cusp vertices have no nontrivial vertices and appear as empty lists.

meshData[{bilby, 1}, D, 1]

Next, we have to assemble the global mass and stiffness matrices. For a given list of numbered finite element vertices numberedTrianglePoints, and a given element mass or stiffness matrix massStiffnessMatrix, the function globalMatrix does most of this. Because the resulting global matrices are rather sparse in average, we will return the matrices as sparse arrays. If needed, the sparse array can easily converted to a dense one.

globalMatrix[numberedTrianglePoints_, massStiffnessMatrix_] :=
Module[{dim, thePoints, 𝒯, j, m, a, n, b},
(* make definition for the inner points *)
 dim = Max[Last /@ Flatten[numberedTrianglePoints]];
 (* indices m, n for the triangles;
    indices a, b for the global matrix *)
 (* all points of the mesh *)
 thePoints = Flatten[numberedTrianglePoints];
 (* now loop over all points of the triangle mesh *) 
 𝒯 = Table[{j, m} = thePoints[[i, 2]]; a = thePoints[[i, 3]];
           (* loop over all vertices from the same triangle *)
           Table[n = numberedTrianglePoints[[j, l, 2, 2]]; 
                 b = numberedTrianglePoints[[j, l, 3]]; 
                 {a, b} -> massStiffnessMatrix[[m, n]],
                 {l, Length[numberedTrianglePoints[[j]]]}],
          {i, 1, Length[thePoints]}];                   
 (* return global matrix as a sparse array *)
 SparseArray[(* sum entries belonging to same vertices *)
             ((#[[1, 1]] -> (Plus @@ (Last /@ #)))& /@ 
                 Split[Sort[Flatten[𝒯]], #1[[1]] === #2[[1]]&])]]

The next graphics show the nonzero elements of the resulting elements of the global mass matrix for the original triangles for finite element orders 2 to 6. We clearly see the overlapping square block-type matrices from the individual elements,

Show[GraphicsArray[
Block[{$DisplayFunction = Identity}, 
Table[ListDensityPlot[1 - Sign[Abs[Normal[
      globalMatrix[meshData[{bilby, 0}, D, femOrder], 
                            MassMatrix[femOrder, {4, 1, 0, 1}]]]]],
             FrameTicks -> False, Mesh -> False], {femOrder, 2, 6}]]]]

Here are two matrices from higher order subdivision (three and four). The overall banded nature of the matrices is clearly visible.

Show[GraphicsArray[
Block[{$DisplayFunction = Identity}, 
Table[ListDensityPlot[1 - Sign[Abs[Normal[
      globalMatrix[meshData[{bilby, ℴ}, D, ℴ], 
                            MassMatrix[ℴ, {2^(2 - ℴ), 1, 0, 1}]]]]],
             FrameTicks -> False, Mesh -> False], {ℴ, 3, 4}]]]]

For further use, we calculate the mass and stiffness matrices ℳ and 𝒮 for various orders. For all the triangles and subdivision orders of our subdivisions of bilby and hawk, we have A=1, B=0, and C=1.

(* use values of J, A, B, and C corresponding to 
   uniform triangle subdivision implemented above *)
Do[JABCRule = Dispatch[{J -> 2^(2 - order), A -> 1, B -> 0, C -> 1}];
   (* reuse already calculated mass and stiffness matrices *)
   ℳ[femOrder][order_] = MassMatrix[femOrder, 
                                     {J, A, B, C}] //. JABCRule;
    𝒮[femOrder][order_] = StiffnessMatrix[femOrder, 
                                          {J, A, B, C}] //. JABCRule, 
   {femOrder, 1, 12}]

Here are some assembly timings, dimensions, and number of nonvanishing matrix elements for the resulting global matrices for various subdivision and finite element orders. The subdivision order 6, finite element order 6 global matrices are of dimensions 7849×7849. And a subdivision order 4, finite element order 12 global matrix is also of dimensions 7849×784, but has 818353 nonvanishing elements.

Table[{order, femOrder,
       Timing[{Dimensions[#][[1]], Length[ArrayRules[#]] - 1}& @   
       globalMatrix[meshData[{bilby, order}, D, femOrder], 
                    ℳ[femOrder][order]]]},
      {order, 6}, {femOrder, 6}]

Now we are ready to test the isospectrality numerically. For a given subdivision order order and a given finite element order femOrder, we now construct the pair of global mass and stiffness matrices bilbyHawkFEMMatrices.

bilbyHawkFEMMatrices[{bilbyHawk_, order_}, DN_, femOrder_] :=
(globalMatrix[meshData[{bilbyHawk, order}, DN, femOrder], 
                               #[femOrder][order]]& /@ {𝒮, ℳ})

We will use the Arnoldi method to calculate the eigenvalues (and the eigenvectors) to make efficient use of the sparse structure of the matrices and to deal with larger matrices. At this stage, we numericalize the global mass and stiffness matrices. The function bilbyHawkEigenvalues calculates the eigenvalues. We allow for options opts to extract higher eigenvalues conveniently.

bilbyHawkEigenvalues[{{bilbyHawk_, order_}, DN_, femOrder_}, 
                     n_, opts___] :=
Module[{M, S},
(* global mass and stiffness matrix *)
{S, M} = bilbyHawkFEMMatrices[{bilbyHawk, order}, DN, femOrder];
(* return eigenvalues in increasing order *)
Reverse @ Eigenvalues[N[{S, M}], n,
                Method -> {"Arnoldi", opts, Tolerance -> 10^-8,
                BasisSize -> Min[Abs[20 n], Dimensions[M]],
                                  MaxIterations -> 10^5}]]

We start with the original seven triangles of bilby and hawk and use a sixth order finite element approximation.

bilbyHawkEigenvalues[{{bilby, 0}, D, 6}, -6]

bilbyHawkEigenvalues[{{hawk, 0}, D, 6}, -6]

The first six eigenvalues agree quite well. Using the sixth subdivision and linear finite elements gives considerably worse results (although the matrix dimension is nearly twice as large). But the actual values of bilby and hawk eigenvalues agree perfectly with each other. The weak performance of low-order finite elements for domains with reentrant corners is a well-known weakness of the finite element method.

bilbyHawkEigenvalues[{{bilby, 6}, D, 1}, -6]

bilbyHawkEigenvalues[{{hawk, 6}, D, 1}, -6]

Because the finite element approximation is basically a variational solution, the resulting values bound the Dirichlet eigenvalues from above. We obtain a few correct digits of the eigenvalues by using higher order triangle subdivisions. In addition, we can extrapolate the sequence of eigenvalues [✶]. The correct value of the first eigenvalue to ten digits is 2.537949998… [✶]. With just seven of the order twelve elements, we obtain the lowest eigenvalue accurate to more 0.0001%.

Table[bilbyHawkEigenvalues[{{bilby, 0}, D, fo}, -1][[1]], {fo, 2, 12}]

Table[bilbyHawkEigenvalues[{{bilby, 1}, D, fo}, -1][[1]], {fo, 2, 12}]

Table[bilbyHawkEigenvalues[{{bilby, 2}, D, fo}, -1][[1]], {fo, 2, 12}]

Assuming that the eigenvalues λ obey a femOrder dependence λ=(Subscript[λ, 0])+(∑_k)(Subscript[c, k])  (femOrder^(-k)), we can refine our predicted value of the lowest eigenvalue. The resulting error is now less than 0.0001%.

Module[{λs = %, l = Length[%], eqs, sol, c, λ0},
 eqs = Table[λ0 + Sum[c[j] fo^-(j + 2), {j, l - 1}] == 
             Rationalize[λs[[k]], 0] /. fo -> k, {k, l}];
 sol = Solve[eqs, Join[{λ0}, Table[c[j], {j, l - 1}]]];
 N[λ0 /. sol[[1]]]]

Because bilby and hawk can be constructed from right angles isosceles triangles, there exist eigenfunctions that are finite linear combinations of products of sine functions and have nodes along coinciding triangle edges. The smallest such eigenvalue for bilby and hawk is ((1^2)+(((1+1))^2))(((π/2))^2)=5/4(π^2). Here we find the numerical eigenvalue nearest to this value and its relative error.

{#, # - 5/4 Pi^2}&[
  bilbyHawkEigenvalues[{{bilby, 1}, D, 12}, 1, Shift -> 12.337][[1]]]

And here are the first few eigenvalues for Neumann boundary conditions. For the solution ψ(x,y)=constant, we have a zero eigenvalue. Again, we see the isospectrality of the bilby and hawk.

bilbyHawkEigenvalues[{{bilby, 0}, N, 6}, -6]

bilbyHawkEigenvalues[{{hawk, 0}, N, 6}, -6]

For mixed boundary conditions (some of the edges have Dirichlet, and some of the edges have Neumann boundary conditions), see [✶].

For visualizing the eigenmodes we need the eigenvectors. The next input defines the function bilbyHawkEigenvectors in an analogous way to the function bilbyHawkEigenvalues. Again, we use the Arnoldi method option in Eigenvectors.

DownValues[bilbyHawkEigenvectors] = 
(* reuse code from above *) DownValues[bilbyHawkEigenvalues] /. 
    {bilbyHawkEigenvalues -> bilbyHawkEigenvectors,
     Eigenvalues -> Eigenvectors};

We start with some 3D plots of the eigenmodes. To obtain good graphics for low subdivision orders of bilby and hawk, we define a triangle subdivision mask that will be applied to each triangle.

subdividedTriangleMask[o_] := subdividedTriangleMask[o] =
Module[{sd, ts, allPoints, dp},
       (* subdivide a triangle into four *)
       sd[{p1_, p2_, p3_}] := 
       With[{p12 = (p1 + p2)/2, p23 = (p2 + p3)/2, p31 = (p3 + p1)/2},
            {{p1, p12, p31}, {p12, p2, p23}, {p23, p3, p31}, 
             {p12, p23, p31}}];
       (* iterate subdivision *)      
       ts = Nest[Flatten[sd /@ #, 1]&, {{{0, 0}, {1, 0}, {0, 1}}}, o];
       (* all points of the subdivision *) 
       allPoints = Union[Flatten[ts, 1]];
       (* identify identical points in the subdivision *)
       dp = Dispatch[MapIndexed[(#1 -> #2[[1]])&, allPoints]];
       (* return list of relative point coordinates and triangles *)
       {allPoints,
        Function[Evaluate[Map[Slot, Polygon /@ ts /. dp, {-1}]]]}]

Now, we define the function eigenfunctionsGraphics3D that generates a 3D plot of the eigenfunctions over bilby and hawk. The argument opp denotes the order of the triangle subdivisions for visualizations and colorFunction is a pure function that will be applied to assign a color to each triangle.

eigenfunctionsGraphics3D[{{bilbyHawk_, order_}, DN_, femOrder_}, 
                         evec_, opp_, colorFunction_, opts___] :=
Module[{(* get triangles *)
        tm = triangleMesh[bilbyHawk][order], λ = numKnots[femOrder],
        (* scale eigenvector *)
        evecN = evec/Plus @@ evec, numberedTrianglePoints = 
                       meshData[{bilbyHawk, order}, DN, femOrder],
        basePointsT, polygonFunction, polys, ST, STF}, 
 (* form the approximate eigenfunctions within a triangle *)
 Do[ST[j][{ξ_, η_}] = Join[toTriangle[tm[[j]], {ξ, η}], {Plus @@ 
    ((evecN[[#[[3]]]]*
      ShapeFunctionList[femOrder, {ξ, η}][[#[[2, 2]]]])& /@ 
                                   numberedTrianglePoints[[j]])}],
    {j, Length[tm]}];
 STF = ST;                 
 (* refinement template *)
 {basePointsT, polygonFunction} = subdividedTriangleMask[opp];    
 (* construct polygons *)
 polys = Table[polygonFunction @@ (STF[j] /@ basePointsT), 
               {j, Length[tm]}];      
 (* return Graphics3D object *)
 {min, max} = {Min[evecN], Max[evecN]};
 (* scaled z-coordinate of polygons *)
 scaledZ[Polygon[l_]] := ((Plus @@ (Last /@ l)/4) - min)/(max - min);
 (* make graphics *)
 Graphics3D[{EdgeForm[], Map[{colorFunction[scaledZ[#]], #}&, polys, {2}],
            {Thickness[0.008], GrayLevel[0], 
             Line[Append[#, 0]& /@ borderLine[bilbyHawk]]}}, 
            opts, PlotRange -> All, BoxRatios -> {1, 1, 0.33},
            BoxStyle -> {Thickness[0.002], GrayLevel[0.5]}]]

Here are some examples of the eigenfunctions.

evecsDBilby = bilbyHawkEigenvectors[{{bilby, 3}, D, 4}, -9];
evecsDHawk = bilbyHawkEigenvectors[{{hawk, 3}, D, 4}, -9];

Do[Show[GraphicsArray[
  eigenfunctionsGraphics3D[{{#1, 3}, D, 4}, #2[[j]], 3, 
                SurfaceColor[Hue[0.78 (1 - #)], 
                             Hue[0.78 # + 0.12], 2.2]&]& @@@
       {{bilby, evecsDBilby}, {hawk, evecsDHawk}}]], {j, 1, 9}];

As already mentioned, the relatively poor performance of the low-order finite element eigenvalues above can be traced back to the reentrant corners of bilby and hawk. There the gradient diverges. With a slight change of the function eigenfunctionsGraphics3D, we can easily display the magnitudes of the gradient vectors over bilby and hawk.

DownValues[eigenfunctionsGradGraphics3D] = 
DownValues[eigenfunctionsGraphics3D] /. 
HoldPattern[eigenfunctionsGraphics3D] :> eigenfunctionsGradGraphics3D /.
(* modify code *)
HoldPattern[Module[{args__}, body_]] :> Module[{args, ξηdx, ξηdy}, body] /.
HoldPattern[STF = ST] :>
 Do[(* chain rule parts in derivative change from {x, y} -> {ξ, η} *)
   {ξηdx, ξηdy} = {D[#, x], D[#, y]}&[{ξ, η} /. 
    Solve[{x, y} == toTriangle[tm[[j]], {ξ, η}], {ξ, η}][[1]]];
    (* norm of gradient *)
    STF[j][{ξ_, η_}] = Join[Most[ST[j][{ξ, η}]], 
         {Sqrt[#.#]&[{ξηdx.#, ξηdy.#}&[({D[#, ξ], D[#, η]}&[
                     ST[j][{ξ, η}][[-1]]])]]}], {j, Length[tm]}];

While the Lagrange element-based calculations give a continuous solution, it is not continuously differentiable. As a result, the following gradient plots have discontinuities along the finite element boundaries. But the following graphics nevertheless clearly show the divergence of the gradient at the reentrant boundary vertices for the first three of the just-visualized eigenfunctions.

Do[Show[GraphicsArray[
  eigenfunctionsGradGraphics3D[{{#1, 3}, D, 4}, #2[[j]], 3, 
                SurfaceColor[Hue[0.78 (1 - #)], 
                             Hue[0.78 # + 0.12], 2.2]&]& @@@
       {{bilby, evecsDBilby}, {hawk, evecsDHawk}}]], {j, 1, 3}];

And here are the corresponding eight lowest nontrivial eigenfunctions eigenfunctions for the Neumann case (the lowest eigenfunction is just a constant).

evecsNBilby = bilbyHawkEigenvectors[{{bilby, 3}, N, 4}, -9];
evecsNHawk = bilbyHawkEigenvectors[{{hawk, 3}, N, 4}, -9];

Do[Show[GraphicsArray[
  eigenfunctionsGraphics3D[{{#1, 3}, N, 4}, #2[[j]], 3, 
                SurfaceColor[Hue[0.78 (1 - #)], 
                             Hue[0.78 # + 0.12], 2.2]&]& @@@
       {{bilby, evecsNBilby}, {hawk, evecsNHawk}}]], {j, 2, 9}];

Contour plots are a better way to visualize higher states. The complicated pattern of minima and maxima is visible at once, without some of the minima and maxima covered by others. To make a contour plot inside the bilby and hawk polygons, we cannot use ListContourPlot directly because it expects a rectangular array of data. While the points of the triangle subdivisions and the finite element base points lie on a Cartesian tensor product grid, they do not fill a rectangle. We could either supplement with function values zero or make individual contour plots in each (subdivided) triangle. First, we will follow the last possibility.

(* function to insert points along the edges of a polygon *)
insertAdditionalPoints[Polygon[l_], ε_] := 
  Module[{n, λ}, Polygon[Join @@ (Function[s, 
               (* segment small enough? *)
               If[(λ = Sqrt[#. #]&[Subtract @@ s]) <= ε, s, 
         (* insert more points *)
         n = Floor[λ/ε] + 1; (* form segments *)
         Table[# + i/n (#2 - #1), {i, 0, n - 1}]& @@ s]] /@ 
                   Partition[Append[l, First[l]], 2, 1])]]

eigenfunctionsContourGraphics[{{bilbyHawk_, order_}, DN_, femOrder_}, 
                               evec_, pp_, zeroContourQ_] :=
Module[{tm = triangleMesh[bilbyHawk][order], λ = numKnots[femOrder],
        evecN = evec/Plus @@ evec, numberedTrianglePoints = 
                         meshData[{bilbyHawk, order}, DN, femOrder],
       mp = {1/3, 1/3}, ε = 10^-12, 𝒸, ST, SF}, 
{min, max} = {Min[evecN], Max[evecN]};        
evecN1 = (evecN - Min[evecN])/(Max[evecN] - Min[evecN]);      
(* form the approximate eigenfunctions within a triangle;
   compile for speed *) 
Do[ST[j] = Compile[{ξ, η}, Evaluate[Plus @@ ((evecN[[#[[3]]]]*
         ShapeFunctionList[femOrder, {ξ, η}][[#[[2, 2]]]])& /@ 
                                   numberedTrianglePoints[[j]])]],
   {j, Length[tm]}];
(* make contour plots within each triangle *)
cps = Table[(* parametrize triangle with rectangular domain *)
            tab = Table[If[ξ === 1, Table[ST[j][1, 0], {pp + 1}],
                        Table[ST[j][ξ, η], {η, 0, 1 - ξ, (1 - ξ)/pp}]], 
                        {ξ, 0, 1, 1/pp}];
            (* options for the two ListContourPlot calls *)
            copts =  Sequence[MeshRange -> {{0, 1}, {0, 1}}, 
                PlotRange -> {min, max}, ColorFunctionScaling -> False,
                ContourSmoothing -> True, DisplayFunction -> Identity];
             (* make contour plot; use some absolute coloring function *)
            cp = ListContourPlot[Transpose[tab], copts, 
                 Contours -> 50, ContourLines -> False, 
                 ColorFunction -> (GrayLevel[(max - #)/(max - min)]&)];
             (* make contour plot of zero0value curve *)
            cpZeroLines = If[zeroContourQ, 
            ListContourPlot[Transpose[tab], copts, Contours -> {0}, 
                 ContourLines -> True, ContourShading -> False,
                 ContourStyle -> {{Thickness[0.004], RGBColor[0, 0, 1]}}], {}];    
            (* make Graphics object *)
            gr = Show[Graphics /@ {cp, cpZeroLines}]; 
            gr1 = gr /. p_Polygon  :> 
                        (* insert points along edges *)
                        insertAdditionalPoints[p, 2^-order/12];    
            (* map to triangle; expand upper polygons *)
            cM = Count[gr, _Polygon, Infinity];            
            ξη[{ξ_, η_}] := {ξ, η (1 - ξ)};
            xy[{ξ_, η_}] = toTriangle[tm[[j]], {ξ, η}];
            expandT[{ξ_, η_}] := mp + (1 + (𝒸++) ε) ({ξ, η} - mp);
            gr1 /. Polygon[l_] :> 
                   (𝒸 = 0; Polygon[xy /@ expandT /@ ξη /@ l]) /.
                   Line[l_] :> Line[xy /@ ξη /@ l],
          {j, Length[tm]}];   
 (* the border *)
 borderGraphics = Graphics[{Hue[0], Thickness[0.008], 
                            Line[borderLine[bilbyHawk]]}];               
 (* form one Graphics object from all contour plots *)
 Show[Append[cps, borderGraphics],
      PlotRange -> All, Frame -> False, AspectRatio -> Automatic] /.  
                            (DisplayFunction -> Identity) -> 
                            (DisplayFunction -> $DisplayFunction)]

Here are the contour plot equivalents of the above 3D plots for bilby.

Show[GraphicsArray[#]]& /@ Partition[
  Block[{$DisplayFunction = Identity}, 
  eigenfunctionsContourGraphics[{{bilby, 3}, D, 4}, #, 4, False]& /@
                                          Take[evecsDBilby, 6]], 3]

The zeroth Neumann eigenfunction is constant and the contour plot shows just the numerical uncertainties. Now, we also show the nodal lines [✶].

Show[GraphicsArray[#]]& /@ Partition[
  Block[{$DisplayFunction = Identity}, 
  eigenfunctionsContourGraphics[{{bilby, 3}, D, 4}, #, 4, True]& /@
                                    Take[evecsDBilby, {2, 7}]], 3]

For applying methods more adapted to the reentrant corners of bilby and hawk to calculate their eigensystems, see [✶]; for experimental data about their eigensystems, see [✶].

For higher states, it is faster and more convenient to use one call to ListContourPlot. For even orders of the triangle subdivision, all triangle vertices lie on a Cartesian x,y tensor product grid (for odd orders, rotating bilby and hawk by π/4 will make the vertices lie on such a grid). We start by creating a square matrix with the function values to be used in ListContourPlot.

eigenfunctionsMat[{{bilbyHawk_, order_}, DN_, femOrder_}, evec_, pp_] :=
Module[{tm = triangleMesh[bilbyHawk][order], λ = numKnots[femOrder],
        evecN = evec/Plus @@ evec, numberedTrianglePoints = 
                         meshData[{bilbyHawk, order}, DN, femOrder],
       ST, toInteger, toIntegerCoordinates, lcpM, δ = 2/2^(order/2)/pp}, 
(* form the approximate eigenfunctions within a triangle;
   compile for speed *) 
Do[ST[j] = Compile[{ξ, η}, Evaluate[Plus @@ ((evecN[[#[[3]]]]*
         ShapeFunctionList[femOrder, {ξ, η}][[#[[2, 2]]]])& /@ 
                                   numberedTrianglePoints[[j]])]],
   {j, Length[tm]}];
(* physical coordinates to grid numberings *)
toIntegerCoordinates[{x_, y_}] := ({x, y} + 3)/δ + 1;
(* empty matrix *)
lcpM = Table[0, {6/δ + 1}, {6/δ + 1}]; 
(* fill in matrix *)
Do[(* go through all triangles *)
   Do[(* pp (1 + pp) points inside each triangle *)
      {a, b} = toIntegerCoordinates[toTriangle[tm[[j]], {ξ, η}]];
      lcpM[[a, b]] = ST[j][ξ, η], {ξ, 0, 1, 1/pp}, {η, 0, 1 - ξ, 1/pp}],
      {j, Length[tm]}]; lcpM]

The function eigenfunctionsContourGraphicsHolistic takes the matrix returned by eigenfunctionsMat and generates a contour plot of the eigenfunctions of bilby and hawk. We cover the parts of the matrix that represent points outside of bilby or hawk with the white polygons outerCoverPolygon.

{outerCoverPolygon[bilby], outerCoverPolygon[hawk]} = Polygon /@ 
 {{{3, 1}, {-1, 1}, {-1, 3}, {-3, 1}, {-1, -1}, {1, -1}, {1, -3}, {3, -1}, 
   {3, -3}, {-3, -3}, {-3, 3}, {3, 3}}, {{-1, 3}, {-1, 1}, {1, 1}, {3, -1}, 
   {1, -1}, {1, -3}, {-3, 1}, {-3, -3}, {3, -3}, {3, 3}, {-1, 3}}};

eigenfunctionsContourGraphicsHolistic[{{bilbyHawk_, order_}, DN_, femOrder_}, 
                                      {evec_, F_}, pp_] :=
Module[{M0, lcp, coverGraphics, borderGraphics}, 
 (* matrix with elongation values *)
 M0 = eigenfunctionsMat[{{bilbyHawk, order}, DN, femOrder}, evec, pp];
 (* make contour plot *)
 lcp = ListContourPlot[Map[F, Transpose[M0], {2}], 
                  PlotRange -> All, DisplayFunction -> Identity,
                  MeshRange -> {{-3, 3}, {-3, 3}}, ContourLines -> False,
                  Contours -> 36];  
 (* cover the outside with white polygons *)
 coverGraphics = Graphics[{GrayLevel[1], outerCoverPolygon[bilbyHawk]}];
 borderGraphics = Graphics[{Hue[0], Thickness[0.008], 
                            Line[borderLine[bilbyHawk]]}];       
 (* make resulting graphic *)
 Show[{lcp, coverGraphics, borderGraphics},
      DisplayFunction -> $DisplayFunction, Frame -> False]]

Here are some of the higher states for the Dirichlet case shown. This time, we show the fourth root of the absolute value of the amplitudes to mimick the appearance of the sand in Chladny figures. The second pair of graphics shows again a state that can be written as a product of sine functions inside each of the seven order zero triangles.

evecsDBilby = bilbyHawkEigenvectors[{{bilby, 4}, D, 10}, 8, Shift -> 101];
evecsDHawk  = bilbyHawkEigenvectors[{{hawk,  4}, D, 10}, 8, Shift -> 101];

Do[Show[GraphicsArray[
     Block[{$DisplayFunction = Identity}, 
     eigenfunctionsContourGraphicsHolistic[{{#1, 4}, D, 10}, 
                                {#2[[j]], Abs[#]^(1/4)&}, 12]& @@@
                             {{bilby, evecsDBilby}, {hawk, evecsDHawk}}]]],
   {j, 8}];

We end with the corresponding graphics for the Neumann case.

evecsNBilby = bilbyHawkEigenvectors[{{bilby, 4}, N, 10}, 8, Shift -> 101];
evecsNHawk  = bilbyHawkEigenvectors[{{hawk,  4}, N, 10}, 8, Shift -> 101];

Do[Show[GraphicsArray[
     Block[{$DisplayFunction = Identity}, 
     eigenfunctionsContourGraphicsHolistic[{{#1, 4}, N, 10}, 
                                {#2[[j]], Abs[#]^(1/4)&}, 12]& @@@
                             {{bilby, evecsNBilby}, {hawk, evecsNHawk}}]]],
   {j, 8}];

After dealing with the inside of these two isospectral drums, one could continue and investigate the Helmholtz (scattering) problem for the outside. See [✶] for some initial results.