Here is a slightly different way to do it. We write a function that converts any PDE (1D/2D/3D) into discretized system matices:
Needs["NDSolve`FEM`"] PDEtoMatrix[{pde_, \[CapitalGamma]___}, u_, r__, o : OptionsPattern[NDSolve`ProcessEquations]] := Module[{ndstate, feData, sd, bcData, methodData, pdeData}, {ndstate} = NDSolve`ProcessEquations[Flatten[{pde, \[CapitalGamma]}], u, Sequence @@ {r}, o]; sd = ndstate["SolutionData"][[1]]; feData = ndstate["FiniteElementData"]; pdeData = feData["PDECoefficientData"]; bcData = feData["BoundaryConditionData"]; methodData = feData["FEMMethodData"]; {DiscretizePDE[pdeData, methodData, sd], DiscretizeBoundaryConditions[bcData, methodData, sd], sd, methodData} ] Example 1: An eigensolver is then something like this:
{dPDE, dBC, sd, md} = PDEtoMatrix[{D[u[t, x, y], t] == Laplacian[u[t, x, y], {x, y}], u[0, x, y] == 0, DirichletCondition[u[t, x, y] == 0, True]}, u, {t, 0, 1}, {x, y} \[Element] Rectangle[]]; l = dPDE["LoadVector"]; s = dPDE["StiffnessMatrix"]; d = dPDE["DampingMatrix"]; constraintMethod = "Remove"; DeployBoundaryConditions[{l, s, d}, dBC, "ConstraintMethod" -> "Remove"]; First[es = Reverse /@ Eigensystem[{s, d}, -4, Method -> "Arnoldi"]] If[constraintMethod === "Remove", es[[2]] = NDSolve`FEM`DirichletValueReinsertion[#, dBC] & /@ es[[2]];]; ifs = ElementMeshInterpolation[sd, #] & /@ es[[2]]; mesh = ifs[[2]]["ElementMesh"]; ContourPlot[#[x, y], {x, y} \[Element] mesh, Frame -> False, ColorFunction -> ColorData["RedBlueTones"]] & /@ ifs 
This can be encapsulated as follows:
Helmholz2D[bdry_, order_] := Module[{dPDE, dBC, sd, md, l, s, d, ifs, es, mesh, constraintMethod}, {dPDE, dBC, sd, md} = PDEtoMatrix[{D[u[t, x, y], t] == Laplacian[u[t, x, y], {x, y}], u[0, x, y] == 0, DirichletCondition[u[t, x, y] == 0, True]}, u, {t, 0, 1}, {x, y} \[Element] bdry]; l = dPDE["LoadVector"]; s = dPDE["StiffnessMatrix"]; d = dPDE["DampingMatrix"]; constraintMethod = "Remove"; DeployBoundaryConditions[{l, s, d}, dBC, "ConstraintMethod" -> "Remove"]; First[es = Reverse /@ Eigensystem[{s, d}, -4, Method -> "Arnoldi"]] If[constraintMethod === "Remove", es[[2]] = NDSolve`FEM`DirichletValueReinsertion[#, dBC] & /@ es[[2]];]; ifs = ElementMeshInterpolation[sd, #] & /@ es[[2]]; mesh = ifs[[2]]["ElementMesh"]; {es, ifs, mesh} ] Example 2: The the remaining problem in the question can then be solved like this:
RR = ImplicitRegion[ x^6 - 5 x^4 y z + 3 x^4 y^2 + 10 x^2 y^3 z + 3 x^2 y^4 - y^5 z + y^6 + z^6 <= 1, {{x, -1.25, 1.25}, {y, -1.25, 1.25}, {z, -1.25, 1.25}}]; mesh = ToElementMesh[RR, "BoundaryMeshGenerator" -> {"RegionPlot", "SamplePoints" -> 31}] mesh["Wireframe"] 
This creates a second order mesh with about 80T tets and 140T nodes. To discretize the PDE we use:
AbsoluteTiming[{dPDE, dBC, sd, md} = PDEtoMatrix[{D[u[t, x, y, z], t] == Laplacian[u[t, x, y, z], {x, y, z}], u[0, x, y, z] == 0, DirichletCondition[u[t, x, y, z] == 0, True]}, u, {t, 0, 1}, {x, y, z} \[Element] mesh]; ] {6.24463, Null} Get the eigen values and vectors:
l = dPDE["LoadVector"]; s = dPDE["StiffnessMatrix"]; d = dPDE["DampingMatrix"]; DeployBoundaryConditions[{l, s, d}, dBC, "ConstraintMethod" -> "Remove"]; AbsoluteTiming[ First[es = Reverse /@ Eigensystem[{s, d}, -4, Method -> "Arnoldi"]] ] {13.484131`, {8.396796994677874`, 16.044484716974942`, 17.453692912770126`, 17.45703443132916`}} Post process / visualize:
ifs = ElementMeshInterpolation[sd, #, "ExtrapolationHandler" -> {(Indeterminate &), "WarningMessage" -> False}] & /@ es[[2]]; Generate slices of the eigenfunctions in the region:
ctrs = Range @@ Join[mm = MinMax[ifs[[2]]["ValuesOnGrid"]], {Abs[Subtract @@ mm]/50}]; levels = Range[-1.25, 1.25, 0.25]; res = ContourPlot[ ifs[[2]][x, y, #], {x, -1.25, 1.25}, {y, -1.25, 1.25}, Frame -> False, ColorFunction -> ColorData["RedBlueTones"], Contours -> ctrs] & /@ levels; Show[{ RegionPlot3D[RR, PlotPoints -> 31, PlotStyle -> Directive[Opacity[0.25]]], Graphics3D[{Opacity[0.25], Flatten[MapThread[ Function[{gr, l}, Cases[gr, _GraphicsComplex] /. GraphicsComplex[coords_, rest__] :> GraphicsComplex[ Join[coords, ConstantArray[{l}, {Length[coords]}], 2], rest]], {res, levels}]]}] }, Boxed -> False, Background -> Gray] 
Imgur http://i.imgur.com/2Pls8K2.gif
Example 3: As a self contained example, let us encapsulate the Helmholz solver
Helmholz3D[bdry_, order_] := Module[{dPDE, dBC, sd, md, l, s, d, ifs, es, mesh, constraintMethod}, {dPDE, dBC, sd, md} = PDEtoMatrix[{D[u[t, x, y, z], t] == Laplacian[u[t, x, y, z], {x, y, z}], u[0, x, y, z] == 0, DirichletCondition[u[t, x, y, z] == 0, True]}, u, {t, 0, 1}, {x, y, z} \[Element] bdry]; l = dPDE["LoadVector"]; s = dPDE["StiffnessMatrix"]; d = dPDE["DampingMatrix"]; constraintMethod = "Remove"; DeployBoundaryConditions[{l, s, d}, dBC, "ConstraintMethod" -> "Remove"]; First[es = Reverse /@ Eigensystem[{s, d}, -4, Method -> "Arnoldi"]] If[constraintMethod === "Remove", es[[2]] = NDSolve`FEM`DirichletValueReinsertion[#, dBC] & /@ es[[2]];]; ifs = ElementMeshInterpolation[sd, #] & /@ es[[2]]; mesh = ifs[[2]]["ElementMesh"]; {es, ifs, mesh} ] and consider
RR = ImplicitRegion[ x^4 + y^4 + z^4 < 1, {{x, -1, 1}, {y, -1, 1}, {z, -1, 1}}] {es, ifs, mesh} = Helmholz3D[RR, nm=4]; res = Table[ If[x^4 + y^4 + z^4 < 1, ifs[[nm]][x, y, z], mm[[1]]], {x, -1.25, 1.25, 0.0125}, {y, -1.25, 1.25, 0.0125}, {z, -1.25, 1.25, 0.0125}]; Image3D[res, ColorFunction -> "RainbowOpacity"] // ImageAdjust 
Example 4 Eigen modes on 3D Knot
Needs["NDSolve`FEM`"] f[t_] := With[{s = 3 t/2}, {(2 + Cos[s]) Cos[t], (2 + Cos[s]) Sin[t], Sin[s]} - {2, 0, 0}] v1[t_] := Cross[f'[t], {0, 0, 1}] // Normalize v2[t_] := Cross[f'[t], v1[t]] // Normalize g[t_, \[Theta]_] := f[t] + (Cos[\[Theta]] v1[t] + Sin[\[Theta]] v2[t])/2 gr = ParametricPlot3D[ Evaluate@g[t, \[Theta]], {t, 0, 4 Pi}, {\[Theta], 0, 2 Pi}, Mesh -> None, MaxRecursion -> 4, Boxed -> False, Axes -> False]; tscale = 4; \[Theta]scale = 0.5;(*scale roughly proportional to \ speeds*)dom = ToElementMesh[ FullRegion[2], {{0, tscale}, {0, \[Theta]scale}},(*domain*) MaxCellMeasure -> {"Area" -> 0.001}]; coords = g[4 Pi #1/tscale, 2 Pi #2/\[Theta]scale] & @@@ dom["Coordinates"];(*apply g*)bmesh2 = ToBoundaryMesh["Coordinates" -> coords, "BoundaryElements" -> dom["MeshElements"]]; emesh2 = ToElementMesh@bmesh2; RR = MeshRegion@emesh2 
{es, ifs, mesh} = Helmholz3D[RR, nm = 4]; then
mm = MinMax[ifs[[nm]]["ValuesOnGrid"]]; Map[{Opacity[0.4], PointSize[0.01], ColorData["Heat"][0.3 + 1/mm[[2]] ifs[[nm + 1]][Sequence @@ #]], Point[#]} &, emesh2["Coordinates"]] // Graphics3D[#, Boxed -> False] & 
