// (M. Bonazzoli, Nov 2015)

load "medit"
load "Element_Mixte3d" //for Edge13d
//load "P012_3d_Modif" //for Edge13d (my file)

// The boundary value problem:
// (sigma = 0, here k is the wavenumber) 
//  -k^2*E + curl(curl E) = 0	in Omega
//  E x n = 0	on x = 0, x = a, y = 0, y = b
//  Curl(E) x n + i*beta n x (E x n) = Ginc   on z = 0
//  Curl(E) x n + i*beta n x (E x n) = 0      on z = c

// Mesh data
int nloc = 4;        // number of segments on the smallest dimension
real a = 0.00254, b = 0.00127, c = 0.01;    // dimensions of the waveguide
//real a = 0.00254, b = 0.00127, c = 0.005;   

// Build the mesh
include "cube.idp"
int mx, my, mz;      // to decide the number of seg in the 3 directions 
mx = a/min(a,b);
my = b/min(a,b);
mz = c/min(a,b);
int[int] NN = [mx*nloc, my*nloc, mz*nloc];  // the number of seg in the 3 directions 
int guide = 1, in = 2, out = 3;             // labels for the waveguide
real [int,int] BB = [[0,a],[0,b],[0,c]];    // bounding box
int [int,int] L = [[guide,guide],[guide,guide],[in,out]];  // labels of the 6 parallelipiped faces
mesh3 Th = Cube(NN,BB,L);                   // build the mesh
//medit("mesh", Th);                        // plot the mesh

// Sol data
real f = 94*10^9;                  // frequence (I think omega=2*pi*f)
real er = 1;                       // dielectric constant
real c0 = 299792458;               // speed of light in vacuum
real k = 2*pi*f*sqrt(er)/c0;       // it's the wavenumber if mu_r=1
int m = 1, n = 0;
real beta = sqrt(k^2-(m*pi/a)^2-(n*pi/b)^2); // it comes from the dispersive relation
// (we assume E(x,y,z) = Etilde(x,y)*exp(-i*beta*z))
real Z = sqrt(er)*120*pi; // ?                        

real ukb = 1/(k^2-beta^2);
func expbz = exp(-1i*beta*z);
func ExTE = (1i*k*Z)*ukb*(n*pi)/b*cos(m*pi*x/a)*sin(n*pi*y/b)*expbz;  // (ricorda: k*Z = mu*omega)
func EyTE = -(1i*k*Z)*ukb*(m*pi)/a*sin(m*pi*x/a)*cos(n*pi*y/b)*expbz;
// For the impedance condition at the waveguide entrance:
func Gix = +2*1i*beta*ExTE;                       
func Giy = +2*1i*beta*EyTE;    
func Giz = 0;
// the sign here is + and in the variational formulation it is - int2d(Th,in)(Ginc*v) (all is written on the lhs)

// Finite element space
fespace Nh(Th, Edge23d);
// Edge13d: edge finite elements of degree 2
// (I called the space I introduced like this because the Nedelec elements of degree 1 are called Edge03d)
Nh<complex> [Ex,Ey,Ez], [vx,vy,vz];   // define the vector field and the test function
// (edge elements are vector elements and they give a subspace of Hcurl)

// Macros
macro Curl(ux,uy,uz) [dy(uz)-dz(uy),dz(ux)-dx(uz),dx(uy)-dy(ux)] // EOM
macro Nvec(ux,uy,uz) [uy*N.z-uz*N.y,uz*N.x-ux*N.z,ux*N.y-uy*N.x] // EOM //uxN 
macro Curlabs(ux,uy,uz) [abs(dy(uz)-dz(uy)),abs(dz(ux)-dx(uz)),abs(dx(uy)-dy(ux))] //EOM

// Variational formulation of the problem to solve
// (sigma = 0, here k is the wavenumber) 
//  -k^2*E + curl(curl E) = 0	in Omega
//  E x n = 0	on x = 0, x = a, y = 0, y = b
//  Curl(E) x n + i*beta n x (E x n) = Ginc   on z = 0
//  Curl(E) x n + i*beta n x (E x n) = 0      on z = c

problem waveguide([Ex,Ey,Ez], [vx,vy,vz], solver=sparsesolver) =
                  int3d(Th)(Curl(Ex,Ey,Ez)'*Curl(vx,vy,vz))
                - int3d(Th)(k^2*[Ex,Ey,Ez]'*[vx,vy,vz])
                + int2d(Th,in,out)(1i*beta*Nvec(Ex,Ey,Ez)'*Nvec(vx,vy,vz))
                - int2d(Th,in)([vx,vy,vz]'*[Gix,Giy,Giz])   
                + on(guide,Ex=0,Ey=0,Ez=0);      
waveguide;  // solve the problem

Nh<complex> [Eex,Eey,Eez] = [ExTE,EyTE,0]; // the exact solution
Nh<complex> [Errx,Erry,Errz]; // the error wrt the exact solution
[Errx,Erry,Errz] = [Eex,Eey,Eez]-[Ex,Ey,Ez];

// Norm of the exact solution
real Hcurlerrsqex, Hcurlerrex, L2errsqex, L2errex;
L2errsqex = int3d(Th)(abs(Eex)^2+abs(Eey)^2+abs(Eez)^2);
Hcurlerrsqex = int3d(Th)(Curlabs(Eex,Eey,Eez)'*Curlabs(Eex,Eey,Eez))+L2errsqex;
Hcurlerrex = sqrt(Hcurlerrsqex);
L2errex = sqrt(L2errsqex);
cout << "Hcurl norm of the exact solution = " << Hcurlerrex << endl;
cout << "L2 of the exact solution = " << L2errex << endl << endl;

// Norm of the error
real Hcurlerrsq, Hcurlerr, L2errsq, L2err;
L2errsq = int3d(Th)(abs(Errx)^2+abs(Erry)^2+abs(Errz)^2);
Hcurlerrsq = int3d(Th)(Curlabs(Errx,Erry,Errz)'*Curlabs(Errx,Erry,Errz))+L2errsq;
Hcurlerr = sqrt(Hcurlerrsq);
L2err = sqrt(L2errsq);
cout << "Hcurl norm of the error = " << Hcurlerr << endl;
cout << "L2 norm of the error = " << L2err << endl << endl;

// Relative errors
cout << "relative Hcurl norm of the error = " << Hcurlerr/Hcurlerrex << endl;
cout << "relative L2 norm of the error = " << L2err/L2errex << endl << endl;

// Plot the real part of the solution
medit("real",Th,[real(Ex),real(Ey),real(Ez)]); // in the medit window press h=help, m=data!!

