Skip to main content
Mathematica

Return to Answer

added 168 characters in body
Source Link
user21
user21
  • 42.2k
  • 8
  • 116
  • 176

Version 11 has both symbolic and numeric eigensolvers, see here for an overview

Version 11 has both symbolic and numeric eigensolvers, see here for an overview

Fixed typo.
Source Link
user21
user21
  • 42.2k
  • 8
  • 116
  • 176
Helmholz2D[bdry_Helmholtz2D[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} ∈ 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} ]

Get the eigen valueseigenvalues and vectors:

Example 3: As a self contained example, let us encapsulate the HelmholzHelmholtz solver

 Helmholz3D[bdry_Helmholtz3D[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} ∈ 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} ] 
RR = ImplicitRegion[ x^4 + y^4 + z^4 < 1, {{x, -1, 1}, {y, -1, 1}, {z, -1, 1}}] {es, ifs, mesh} = Helmholz3D[RRHelmholtz3D[RR, nm=4]; mm = MinMax[ifs[[nm]]["ValuesOnGrid"]]; Map[{Opacity[0.4], PointSize[0.01], ColorData["Heat"][0.3 + 1/mm[[2]] ifs[[nm]][Sequence @@ #]], Point[#]} &, mesh["Coordinates"]] // Graphics3D[#, Boxed -> False] & 
{es, ifs, mesh} = Helmholz3D[RRHelmholtz3D[RR, nm = 4]; 
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} ∈ 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} ]

Get the eigen values and vectors:

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} ∈ 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} ] 
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"]]; Map[{Opacity[0.4], PointSize[0.01], ColorData["Heat"][0.3 + 1/mm[[2]] ifs[[nm]][Sequence @@ #]], Point[#]} &, mesh["Coordinates"]] // Graphics3D[#, Boxed -> False] & 
{es, ifs, mesh} = Helmholz3D[RR, nm = 4]; 
Helmholtz2D[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} ∈ 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} ]

Get the eigenvalues and vectors:

Example 3: As a self contained example, let us encapsulate the Helmholtz solver

 Helmholtz3D[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} ∈ 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} ] 
RR = ImplicitRegion[ x^4 + y^4 + z^4 < 1, {{x, -1, 1}, {y, -1, 1}, {z, -1, 1}}] {es, ifs, mesh} = Helmholtz3D[RR, nm=4]; mm = MinMax[ifs[[nm]]["ValuesOnGrid"]]; Map[{Opacity[0.4], PointSize[0.01], ColorData["Heat"][0.3 + 1/mm[[2]] ifs[[nm]][Sequence @@ #]], Point[#]} &, mesh["Coordinates"]] // Graphics3D[#, Boxed -> False] & 
{es, ifs, mesh} = Helmholtz3D[RR, nm = 4];
deleted 129 characters in body
Source Link
chris
chris
  • 23.4k
  • 5
  • 63
  • 154
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} ] 
 {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 
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} ] 
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} 
 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} ] 
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]]Cos[θ] v1[t] + Sin[\[Theta]]Sin[θ] 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θscale = 0.5;(*scale roughly proportional to \ speeds*)dom = ToElementMesh[ FullRegion[2], {{0, tscale}, {0, \[Theta]scaleθscale}},(*domain*) MaxCellMeasure -> {"Area" -> 0.001}]; coords = g[4 Pi #1/tscale, 2 Pi #2/\[Theta]scale]θscale] & @@@ dom["Coordinates"];(*apply g*)bmesh2 = ToBoundaryMesh["Coordinates" -> coords, "BoundaryElements" -> dom["MeshElements"]]; emesh2 = ToElementMesh@bmesh2; RR = MeshRegion@emesh2 
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} ] 
 {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 
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} ] 
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} 
 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} ] 
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 
Needs["NDSolve`FEM`"] PDEtoMatrix[{pde_, Γ___}, u_, r__, o : OptionsPattern[NDSolve`ProcessEquations]] := Module[{ndstate, feData, sd, bcData, methodData, pdeData}, {ndstate} = NDSolve`ProcessEquations[Flatten[{pde, Γ}], 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} ] 
 {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} ∈ 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} ∈ mesh, Frame -> False, ColorFunction -> ColorData["RedBlueTones"]] & /@ ifs 
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} ∈ 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} ] 
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} ∈ mesh]; ] {6.24463, Null} 
 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} ∈ 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} ] 
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_, θ_] := f[t] + (Cos[θ] v1[t] + Sin[θ] v2[t])/2 gr = ParametricPlot3D[ Evaluate@g[t, θ], {t, 0, 4 Pi}, {θ, 0, 2 Pi}, Mesh -> None, MaxRecursion -> 4, Boxed -> False, Axes -> False]; tscale = 4; θscale = 0.5;(*scale roughly proportional to \ speeds*)dom = ToElementMesh[ FullRegion[2], {{0, tscale}, {0, θscale}},(*domain*) MaxCellMeasure -> {"Area" -> 0.001}]; coords = g[4 Pi #1/tscale, 2 Pi #2/θscale] & @@@ dom["Coordinates"];(*apply g*)bmesh2 = ToBoundaryMesh["Coordinates" -> coords, "BoundaryElements" -> dom["MeshElements"]]; emesh2 = ToElementMesh@bmesh2; RR = MeshRegion@emesh2
deleted 4 characters in body
Source Link
chris
chris
  • 23.4k
  • 5
  • 63
  • 154
deleted 4 characters in body
Source Link
chris
chris
  • 23.4k
  • 5
  • 63
  • 154
added 4 characters in body
Source Link
chris
chris
  • 23.4k
  • 5
  • 63
  • 154
deleted 228 characters in body
Source Link
chris
chris
  • 23.4k
  • 5
  • 63
  • 154
deleted 228 characters in body
Source Link
chris
chris
  • 23.4k
  • 5
  • 63
  • 154
deleted 21 characters in body
Source Link
chris
chris
  • 23.4k
  • 5
  • 63
  • 154
added 98 characters in body
Source Link
chris
chris
  • 23.4k
  • 5
  • 63
  • 154
addes DirchletCondition code
Source Link
user21
user21
  • 42.2k
  • 8
  • 116
  • 176
added 90 characters in body
Source Link
chris
chris
  • 23.4k
  • 5
  • 63
  • 154
added 1359 characters in body
Source Link
chris
chris
  • 23.4k
  • 5
  • 63
  • 154
deleted 870 characters in body
Source Link
chris
chris
  • 23.4k
  • 5
  • 63
  • 154
added 876 characters in body
Source Link
chris
chris
  • 23.4k
  • 5
  • 63
  • 154
edited body
Source Link
chris
chris
  • 23.4k
  • 5
  • 63
  • 154
added 44 characters in body
Source Link
chris
chris
  • 23.4k
  • 5
  • 63
  • 154
added 1218 characters in body
Source Link
chris
chris
  • 23.4k
  • 5
  • 63
  • 154
added pic
Source Link
user21
user21
  • 42.2k
  • 8
  • 116
  • 176
Bounty Awarded with 100 reputation awarded by chris
Source Link
user21
user21
  • 42.2k
  • 8
  • 116
  • 176