NDeigensystem returns error due to mesh discretization when calculating vibrations of a cantilever

Question

There was an transcription error in the code I provided in a previous post:

https://mathematica.stackexchange.com/questions/105713/ndeigensystem-returns-complex-non-hermitian-error-for-the-calculation-of-vibrati

This resulted in it being marked as a duplicate. My apologies. I am still facing the issue I outlined in the post and would appreciate the assistance. I have corrected the question, it is provided below.

We are trying to use NDEigensystem to solve for the vibration of a cantilever with a triangular cross-section. The relevant PDE is $\mu \nabla^2 \vec u + (\lambda + \mu) \nabla(\nabla \cdot \vec u) = \rho \omega^2 \vec u$ where $\vec u(x,y,z)$ is a 3D vector function representing the vibrational displacement field, $\rho$ is the material density, and $\lambda$ and $\mu$ are elastic constants.

The geometry is depicted below.

TriangularPrism[W_, L_, H_, MAT_: DiagonalMatrix[{1, 1, 1}]] := 
Prism[MAT.# & /@ {{0, W*Sqrt[3]/3, 0}, {-W*1/2, -W*Sqrt[3]/6, 0}, {W*1/2, -W*Sqrt[3]/6, 0}, 
    {0, W*Sqrt[3]/3, L + H}, {-W*1/2, -W*Sqrt[3]/6, L + H}, {W*1/2, -W*Sqrt[3]/6, L + H}}];
CircularSubstrate[R_, H_] := ImplicitRegion[x^2 + y^2 <= R^2 && 0 <= z <= H, {x, y, z}];
CompoundStructure[A_, B_] := RegionUnion[A, B]

In the above example, the long axis of the cantilever extends from 0.5 to 1. The cantilever has an equilateral triangle cross-section of side length 0.4 The cylindrical substrate has radius 4 and height 1/2.

The code used to calculate the vibration of a cantilever of an arbitrary side length and height is given below (L=Cavity Height, W=Cavity Width, R=Substrate Radius, H=Substrate Height)

{L, W, R, H, GridNo, EignNo} = {0.8, 0.7, 4, 1/2, 0.01, 17};

LIxx = FullSimplify[({{λ + 2 μ, 0, 0}, {0, μ, 0}, {0, 0, μ}})];
LIxy = FullSimplify[({{0, λ, 0}, {μ, 0, 0}, {0, 0, 0}})];
LIxz = FullSimplify[({{0, 0, λ}, {0, 0, 0}, {μ, 0, 0}})];
LIyx = FullSimplify[({{0, μ, 0}, {λ, 0, 0}, {0, 0, 0}})];
LIyy = FullSimplify[({{μ, 0, 0}, {0, λ + 2 μ, 0}, {0, 0, μ}})];
LIyz = FullSimplify[({{0, 0, 0}, {0, 0, λ}, {0, μ, 0}})];
LIzx = FullSimplify[({{0, 0, μ}, {0, 0, 0}, {λ, 0, 0}})];
LIzy = FullSimplify[({{0, 0, 0}, {0, 0, μ}, {0, λ, 0}})];
LIzz = FullSimplify[({{μ, 0, 0}, {0, μ, 0}, {0, 0, λ + 2 μ}})];

{{LIxx, LIxy, LIxz}, {LIyx, LIyy, LIyz}, {LIzx, LIzy, LIzz}} // MatrixForm
PhononEquation3D = {
Inactive[Div][(LIxx.Inactive[Grad][u[x, y, z], {x, y, z}]), {x, y, z}]+
Inactive[Div][(LIxy.Inactive[Grad][v[x, y, z], {x, y, z}]), {x, y, z}]+
Inactive[Div][(LIxz.Inactive[Grad][w[x, y, z], {x, y, z}]), {x, y, z}],

Inactive[Div][(LIyy.Inactive[Grad][v[x, y, z], {x, y, z}]), {x, y, z}]+
Inactive[Div][(LIyx.Inactive[Grad][u[x, y, z], {x, y, z}]), {x, y, z}]+
Inactive[Div][(LIyz.Inactive[Grad][w[x, y, z], {x, y, z}]), {x, y, z}],

Inactive[Div][(LIzz.Inactive[Grad][w[x, y, z], {x, y, z}]), {x, y, z}]+
Inactive[Div][(LIzx.Inactive[Grad][u[x, y, z], {x, y, z}]), {x, y, z}]+
Inactive[Div][(LIzy.Inactive[Grad][v[x, y, z], {x, y, z}]), {x, y, z}]} /. 
    {λ  -> -85*10^(3)/((2*\[Pi])^2*3512), μ -> -536*10^(3)/((2*\[Pi])^2*3512)};
Cavity = CompoundStructure[TriangularPrism[W, L, H], CircularSubstrate[R, H]];
Pillar = TriangularPrism[W, L, H];

ΓCircularSubstrateClampedEdge = DirichletCondition[{u[x, y, z] == 0, 
    v[x, y, z] == 0, w[x, y, z] == 0}, x^2 + y^2 == R^2 && 0 <= z <= H];
ΓSubstrateClampedBase = DirichletCondition[{u[x, y, z] == 0, 
    v[x, y, z] == 0, w[x, y, z] == 0}, z == 0];
ΓPrismClampedBase = DirichletCondition[{u[x, y, z] == 0, 
    v[x, y, z] == 0, w[x, y, z] == 0}, z == 1/2];

mparams = {"Eigensystem" -> {"FEAST","Interval" -> {(50 - 0.145)^2, (50 + 0.145)^2}, 
    "SubspaceSize" -> 50}, "SpatialDiscretization" -> {"FiniteElement", 
    {"MeshOptions" -> {"MaxCellMeasure" -> GridNO}}}};

ωf =  NDEigensystem[{PhononEquation3D, ΓSubstrateClampedBase}, 
    {u, v, w}, {x, y, z} ∈ Cavity, EignNo, Method -> mparams]

In general this code works, except for some cavity geometries. For example, the height and width specified in the code above returns the following error:

We have diagnosed that the error is the result of the spatial discretization of the cantilever region by the NDEigensystem function. This was indicated by the sensitivity of the solver to the dimensions of the cavity. For example, the solver fails when applied to cantilever with height 0.7 and width 0.8 whilst it works for cantilever dimensions 0.7002 and 0.8002. This is particularly problematic for our purposes, as we seek to automate the calculation of the properties of a cantilevers with many differing heights and widths. In particular, in a parallel calculation this error causes the entire calculation to crash.

Any ideas on how to approach this problem? Thanks, again.