• API
• FAQ
• Tools
• Archive
SHARE
TWEET

# Untitled

a guest Jun 15th, 2019 63 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. << NDSolve`FEM`
2.
3. r = 0.8;
4. outerCirclePoints =
5.   With[{r = 2.},
6.    Table[{r Cos[[Theta]], r Sin[[Theta]]}, {[Theta],
7.      Range[0, 2 [Pi], 0.05 [Pi]] // Most}]];
8. innerCirclePoints =
9.   With[{r = r},
10.    Table[{r Cos[[Theta]], r Sin[[Theta]]}, {[Theta],
11.      Range[0, 2 [Pi], 0.08 [Pi]] // Most}]];
12.
13. bmesh = ToBoundaryMesh[
14.    "Coordinates" -> Join[outerCirclePoints, innerCirclePoints],
15.    "BoundaryElements" -> {LineElement[
16.       Riffle[Range[Length@outerCirclePoints],
17.         RotateLeft[Range[Length@outerCirclePoints], 1]] //
18.        Partition[#, 2] &],
19.      LineElement[
20.       Riffle[Range[Length@outerCirclePoints + 1,
21.          Length@Join[outerCirclePoints, innerCirclePoints]],
22.         RotateLeft[
23.          Range[Length@outerCirclePoints + 1,
24.           Length@Join[outerCirclePoints, innerCirclePoints]], 1]] //
25.        Partition[#, 2] &]}];
26.
27. mesh = ToElementMesh[bmesh];
28. {bmesh["Wireframe"], mesh["Wireframe"]}
29.
30. glass = 1.45; air = 1.; k0 = (2 [Pi])/1.55;
31. n[x_, y_] := If[x^2 + y^2 <= r^2, glass, air]
32.
33. helm = !(
34. *SubsuperscriptBox[([Del]), ({x, y}), (2)](u[x, y])) +
35.    n[x, y]^2*k0^2*u[x, y];
36. boundary = DirichletCondition[u[x, y] == 0., True];
37.
38. (*region=ImplicitRegion[x^2+y^2[LessEqual]2.^2,{x,y}];*)
39.
40. {vals, funs} =
41.   NDEigensystem[{helm, boundary}, u[x, y], {x, y} [Element] mesh, 1,
42.    Method -> {"Eigensystem" -> {"FEAST",
43.        "Interval" -> {k0^2, glass^2 k0^2}}}];
44. vals
45.
46.  Table[Plot3D[funs[[i]], {x, y} [Element] mesh, PlotRange -> All,
47.   PlotLabel -> vals[[i]]], {i, Length[vals]}]
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.

Top