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}, -order, 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"]
![enter image description here][2]
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]
![enter image description here][3]
**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];
mm = MinMax[ifs[[nm]]["ValuesOnGrid"]];
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]][Sequence @@ #]],
Point[#]} &, emesh2["Coordinates"]] //
Graphics3D[#, Boxed -> False] &
[1]: https://i.sstatic.net/Vznyb.png
[2]: https://i.sstatic.net/Vt2dv.png
[3]: https://i.sstatic.net/X54De.png