1// (M. Bonazzoli, Nov 2015)
2
3load "msh3"
4load "medit"
5load "Element_Mixte3d" //for Edge13d
6//load "P012_3d_Modif" //for Edge13d (my file)
7
8// The boundary value problem:
9// (sigma = 0, here k is the wavenumber)
10//  -k^2*E + curl(curl E) = 0	in Omega
11//  E x n = 0	on x = 0, x = a, y = 0, y = b
12//  Curl(E) x n + i*beta n x (E x n) = Ginc   on z = 0
13//  Curl(E) x n + i*beta n x (E x n) = 0      on z = c
14
15// Mesh data
16int nloc = 4;        // number of segments on the smallest dimension
17real a = 0.00254, b = 0.00127, c = 0.01;    // dimensions of the waveguide
18//real a = 0.00254, b = 0.00127, c = 0.005;
19
20// Build the mesh
21include "cube.idp"
22int mx, my, mz;      // to decide the number of seg in the 3 directions
23mx = a/min(a,b);
24my = b/min(a,b);
25mz = c/min(a,b);
26int[int] NN = [mx*nloc, my*nloc, mz*nloc];  // the number of seg in the 3 directions
27int guide = 1, in = 2, out = 3;             // labels for the waveguide
28real [int,int] BB = [[0,a],[0,b],[0,c]];    // bounding box
29int [int,int] L = [[guide,guide],[guide,guide],[in,out]];  // labels of the 6 parallelipiped faces
30mesh3 Th = Cube(NN,BB,L);                   // build the mesh
31//medit("mesh", Th);                        // plot the mesh
32
33// Sol data
34real f = 94*10^9;                  // frequence (I think omega=2*pi*f)
35real er = 1;                       // dielectric constant
36real c0 = 299792458;               // speed of light in vacuum
37real k = 2*pi*f*sqrt(er)/c0;       // it's the wavenumber if mu_r=1
38int m = 1, n = 0;
39real beta = sqrt(k^2-(m*pi/a)^2-(n*pi/b)^2); // it comes from the dispersive relation
40// (we assume E(x,y,z) = Etilde(x,y)*exp(-i*beta*z))
41real Z = sqrt(er)*120*pi; // ?
42
43real ukb = 1/(k^2-beta^2);
44func expbz = exp(-1i*beta*z);
45func ExTE = (1i*k*Z)*ukb*(n*pi)/b*cos(m*pi*x/a)*sin(n*pi*y/b)*expbz;  // (ricorda: k*Z = mu*omega)
46func EyTE = -(1i*k*Z)*ukb*(m*pi)/a*sin(m*pi*x/a)*cos(n*pi*y/b)*expbz;
47// For the impedance condition at the waveguide entrance:
48func Gix = +2*1i*beta*ExTE;
49func Giy = +2*1i*beta*EyTE;
50func Giz = 0;
51// the sign here is + and in the variational formulation it is - int2d(Th,in)(Ginc*v) (all is written on the lhs)
52
53// Finite element space
54fespace Nh(Th, Edge13d);
55// Edge13d: edge finite elements of degree 2
56// (I called the space I introduced like this because the Nedelec elements of degree 1 are called Edge03d)
57Nh<complex> [Ex,Ey,Ez], [vx,vy,vz];   // define the vector field and the test function
58// (edge elements are vector elements and they give a subspace of Hcurl)
59
60// Macros
61macro Curl(ux,uy,uz) [dy(uz)-dz(uy),dz(ux)-dx(uz),dx(uy)-dy(ux)] // EOM
62macro Nvec(ux,uy,uz) [uy*N.z-uz*N.y,uz*N.x-ux*N.z,ux*N.y-uy*N.x] // EOM //uxN
63macro Curlabs(ux,uy,uz) [abs(dy(uz)-dz(uy)),abs(dz(ux)-dx(uz)),abs(dx(uy)-dy(ux))] //EOM
64
65// Variational formulation of the problem to solve
66// (sigma = 0, here k is the wavenumber)
67//  -k^2*E + curl(curl E) = 0	in Omega
68//  E x n = 0	on x = 0, x = a, y = 0, y = b
69//  Curl(E) x n + i*beta n x (E x n) = Ginc   on z = 0
70//  Curl(E) x n + i*beta n x (E x n) = 0      on z = c
71
72problem waveguide([Ex,Ey,Ez], [vx,vy,vz], solver=sparsesolver) =
73                  int3d(Th)(Curl(Ex,Ey,Ez)'*Curl(vx,vy,vz))
74                - int3d(Th)(k^2*[Ex,Ey,Ez]'*[vx,vy,vz])
75                + int2d(Th,in,out)(1i*beta*Nvec(Ex,Ey,Ez)'*Nvec(vx,vy,vz))
76                - int2d(Th,in)([vx,vy,vz]'*[Gix,Giy,Giz])
77                + on(guide,Ex=0,Ey=0,Ez=0);
78waveguide;  // solve the problem
79
80Nh<complex> [Eex,Eey,Eez] = [ExTE,EyTE,0]; // the exact solution
81Nh<complex> [Errx,Erry,Errz]; // the error wrt the exact solution
82[Errx,Erry,Errz] = [Eex,Eey,Eez]-[Ex,Ey,Ez];
83
84// Norm of the exact solution
85real Hcurlerrsqex, Hcurlerrex, L2errsqex, L2errex;
86L2errsqex = int3d(Th)(abs(Eex)^2+abs(Eey)^2+abs(Eez)^2);
87Hcurlerrsqex = int3d(Th)(Curlabs(Eex,Eey,Eez)'*Curlabs(Eex,Eey,Eez))+L2errsqex;
88Hcurlerrex = sqrt(Hcurlerrsqex);
89L2errex = sqrt(L2errsqex);
90cout << "Hcurl norm of the exact solution = " << Hcurlerrex << endl;
91cout << "L2 of the exact solution = " << L2errex << endl << endl;
92
93// Norm of the error
94real Hcurlerrsq, Hcurlerr, L2errsq, L2err;
95L2errsq = int3d(Th)(abs(Errx)^2+abs(Erry)^2+abs(Errz)^2);
96Hcurlerrsq = int3d(Th)(Curlabs(Errx,Erry,Errz)'*Curlabs(Errx,Erry,Errz))+L2errsq;
97Hcurlerr = sqrt(Hcurlerrsq);
98L2err = sqrt(L2errsq);
99cout << "Hcurl norm of the error = " << Hcurlerr << endl;
100cout << "L2 norm of the error = " << L2err << endl << endl;
101
102// Relative errors
103cout << "relative Hcurl norm of the error = " << Hcurlerr/Hcurlerrex << endl;
104cout << "relative L2 norm of the error = " << L2err/L2errex << endl << endl;
105
106// Plot the real part of the solution
107medit("real",Th,[real(Ex),real(Ey),real(Ez)]); // in the medit window press h=help, m=data!!
108
109