-- FreeFem++ v4.12 (Wed Feb 15 11:36:54 CET 2023 - git v4.12-27-g58d32866) file : C:\Users\Ming\Downloads\script (20).edp Load: lg_fem lg_mesh lg_mesh3 eigenvalue 1 : load "msh3"(load: loadLibary C:\Program Files (x86)\FreeFem++\\.\msh3 = 0) 2 : load "medit"(load: loadLibary C:\Program Files (x86)\FreeFem++\\.\medit = 0) 3 : load "Element_Mixte3d"(load: loadLibary C:\Program Files (x86)\FreeFem++\\.\Element_Mixte3d = 0) //for Edge13d 4 : //load "P012_3d_Modif" //for Edge13d (my file) 5 : // The boundary value problem: 6 : // (sigma = 0, here k is the wavenumber) 7 : // -k^2*E + curl(curl E) = 0 in Omega 8 : // E x n = 0 on x = 0, x = a, y = 0, y = b 9 : // Curl(E) x n + i*beta n x (E x n) = Ginc on z = 0 10 : // Curl(E) x n + i*beta n x (E x n) = 0 on z = c 11 : 12 : // Mesh data 13 : int nloc = 4; // number of segments on the smallest dimension 14 : real a = 0.00254, b = 0.00127, c = 0.01; // dimensions of the waveguide 15 : //real a = 0.00254, b = 0.00127, c = 0.005; 16 : 17 : // Build the mesh 18 : include "cube.idp"load "msh3" (already loaded: msh3) 2 : load "medit" (already loaded: medit) 3 : // ! basic functions to build regular mesh of a cube 4 : /* 5 : mesh3 Cube(NN,BB,L); 6 : -- build the surface mesh of a 3d box 7 : where: for exqmple: 8 : int[int] NN=[nx,ny,nz]; // the number of seg in the 3 direction 9 : real [int,int] BB=[[xmin,xmax],[ymin,ymax],[zmin,zmax]]; // bounding bax 10 : int [int,int] L=[[1,2],[3,4],[5,6]]; // the label of the 6 face left,right, front, back, down, right 11 : */ 12 : func mesh3 Cube(int[int] & NN,real[int,int] &BB ,int[int,int] & L) 13 : { 14 : // first build the 6 faces of the hex. 15 : real x0=BB(0,0),x1=BB(0,1); 16 : real y0=BB(1,0), *** Warning The identifier y0 hide a Global identifier y1=BB(1,1); *** Warning The identifier y1 hide a Global identifier 17 : real z0=BB(2,0),z1=BB(2,1); 18 : 19 : int nx=NN[0],ny=NN[1],nz=NN[2]; 20 : mesh Thx = square(nx,ny,[x0+(x1-x0)*x,y0+(y1-y0)*y]); 21 : 22 : int[int] rup=[0,L(2,1)], rdown=[0,L(2,0)], 23 : rmid=[1,L(1,0), 2,L(0,1), 3, L(1,1), 4, L(0,0) ]; 24 : mesh3 Th=buildlayers(Thx,nz, zbound=[z0,z1], 25 : labelmid=rmid, labelup = rup, 26 : labeldown = rdown); 27 : 28 : return Th; 29 : } 30 : func mesh3 Cube(int Nx,int Ny,int Nz) 31 : { 32 : int[int] NN=[Nx,Ny,Nz]; 33 : real [int,int] BB=[[0,1],[0,1],[0,1]]; 34 : int[int,int] LL=[[1,2],[3,4],[5,6]]; 35 : return Cube(NN,BB,LL); 36 : } 37 : 38 : 39 : 40 : 19 : int mx, my, mz; // to decide the number of seg in the 3 directions 20 : mx = a/min(a,b); 21 : my = b/min(a,b); 22 : mz = c/min(a,b); 23 : int[int] NN = [mx*nloc, my*nloc, mz*nloc]; // the number of seg in the 3 directions 24 : int guide = 1, in = 2, out = 3; // labels for the waveguide 25 : real [int,int] BB = [[0,a],[0,b],[0,c]]; // bounding box 26 : int [int,int] L = [[guide,guide],[guide,guide],[in,out]]; // labels of the 6 parallelipiped faces 27 : mesh3 Th = Cube(NN,BB,L); // build the mesh 28 : //medit("mesh", Th); // plot the mesh 29 : 30 : // Sol data 31 : real f = 94*10^9; // frequence (I think omega=2*pi*f) 32 : real er = 1; // dielectric constant 33 : real c0 = 299792458; // speed of light in vacuum 34 : real k = 2*pi*f*sqrt(er)/c0; // it's the wavenumber if mu_r=1 35 : int m = 1, n = 0; 36 : real beta = sqrt(k^2-(m*pi/a)^2-(n*pi/b)^2); // it comes from the dispersive relation 37 : // (we assume E(x,y,z) = Etilde(x,y)*exp(-i*beta*z)) 38 : real Z = sqrt(er)*120*pi; // ? 39 : 40 : real ukb = 1/(k^2-beta^2); 41 : func expbz = exp(-1i*beta*z); 42 : 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) 43 : func EyTE = -(1i*k*Z)*ukb*(m*pi)/a*sin(m*pi*x/a)*cos(n*pi*y/b)*expbz; 44 : // For the impedance condition at the waveguide entrance: 45 : func Gix = +2*1i*beta*ExTE; 46 : func Giy = +2*1i*beta*EyTE; 47 : func Giz = 0; 48 : // the sign here is + and in the variational formulation it is - int2d(Th,in)(Ginc*v) (all is written on the lhs) 49 : 50 : // Finite element space 51 : fespace Nh(Th, Edge23d); 52 : // Edge13d: edge finite elements of degree 2 53 : // (I called the space I introduced like this because the Nedelec elements of degree 1 are called Edge03d) 54 : Nh [Ex,Ey,Ez], [vx,vy,vz]; // define the vector field and the test function 55 : // (edge elements are vector elements and they give a subspace of Hcurl) 56 : 57 : // Macros 58 : macro Curl(ux,uy,uz) [dy(uz)-dz(uy),dz(ux)-dx(uz),dx(uy)-dy(ux)] ) // EOM 59 : 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 60 : macro Curlabs(ux,uy,uz) [abs(dy(uz)-dz(uy)),abs(dz(ux)-dx(uz)),abs(dx(uy)-dy(ux))] ) //EOM 61 : 62 : // Variational formulation of the problem to solve 63 : // (sigma = 0, here k is the wavenumber) 64 : // -k^2*E + curl(curl E) = 0 in Omega 65 : // E x n = 0 on x = 0, x = a, y = 0, y = b 66 : // Curl(E) x n + i*beta n x (E x n) = Ginc on z = 0 67 : // Curl(E) x n + i*beta n x (E x n) = 0 on z = c 68 : 69 : problem waveguide([Ex,Ey,Ez], [vx,vy,vz], solver=sparsesolver) = 70 : int3d(Th)(Curl(Ex,Ey,Ez) [dy(Ez)-dz(Ey),dz(Ex)-dx(Ez),dx(Ey)-dy(Ex)] '*Curl(vx,vy,vz) [dy(vz)-dz(vy),dz(vx)-dx(vz),dx(vy)-dy(vx)] ) 71 : - int3d(Th)(k^2*[Ex,Ey,Ez]'*[vx,vy,vz]) 72 : + int2d(Th,in,out)(1i*beta*Nvec(Ex,Ey,Ez) [Ey*N.z-Ez*N.y,Ez*N.x-Ex*N.z,Ex*N.y-Ey*N.x] '*Nvec(vx,vy,vz) [vy*N.z-vz*N.y,vz*N.x-vx*N.z,vx*N.y-vy*N.x] ) 73 : - int2d(Th,in)([vx,vy,vz]'*[Gix,Giy,Giz]) 74 : + on(guide,Ex=0,Ey=0,Ez=0); 75 : waveguide; 76 : 77 : Nh [Eex,Eey,Eez] = [ExTE,EyTE,0]; // the exact solution 78 : Nh [Errx,Erry,Errz]; // the error wrt the exact solution 79 : [Errx,Erry,Errz] = [Eex,Eey,Eez]-[Ex,Ey,Ez]; 80 : 81 : // Norm of the exact solution 82 : real Hcurlerrsqex, Hcurlerrex, L2errsqex, L2errex; 83 : L2errsqex = int3d(Th)(abs(Eex)^2+abs(Eey)^2+abs(Eez)^2); 84 : Hcurlerrsqex = int3d(Th)(Curlabs(Eex,Eey,Eez) [abs(dy(Eez)-dz(Eey)),abs(dz(Eex)-dx(Eez)),abs(dx(Eey)-dy(Eex))] '*Curlabs(Eex,Eey,Eez) [abs(dy(Eez)-dz(Eey)),abs(dz(Eex)-dx(Eez)),abs(dx(Eey)-dy(Eex))] )+L2errsqex; 85 : Hcurlerrex = sqrt(Hcurlerrsqex); 86 : L2errex = sqrt(L2errsqex); 87 : cout << "Hcurl norm of the exact solution = " << Hcurlerrex << endl; 88 : cout << "L2 of the exact solution = " << L2errex << endl << endl; 89 : 90 : // Norm of the error 91 : real Hcurlerrsq, Hcurlerr, L2errsq, L2err; 92 : L2errsq = int3d(Th)(abs(Errx)^2+abs(Erry)^2+abs(Errz)^2); 93 : Hcurlerrsq = int3d(Th)(Curlabs(Errx,Erry,Errz) [abs(dy(Errz)-dz(Erry)),abs(dz(Errx)-dx(Errz)),abs(dx(Erry)-dy(Errx))] '*Curlabs(Errx,Erry,Errz) [abs(dy(Errz)-dz(Erry)),abs(dz(Errx)-dx(Errz)),abs(dx(Erry)-dy(Errx))] )+L2errsq; 94 : Hcurlerr = sqrt(Hcurlerrsq); 95 : L2err = sqrt(L2errsq); 96 : cout << "Hcurl norm of the error = " << Hcurlerr << endl; 97 : cout << "L2 norm of the error = " << L2err << endl << endl; 98 : 99 : // Relative errors 100 : cout << "relative Hcurl norm of the error = " << Hcurlerr/Hcurlerrex << endl; 101 : cout << "relative L2 norm of the error = " << L2err/L2errex << endl << endl; 102 : 103 : 104 : // Plot the real part of the solution 105 : medit("real",Th,[real(Ex),real(Ey),real(Ez)]); // in the medit window press h=help, m=data!! 106 : sizestack + 1024 =15008 ( 13984 ) -- Square mesh : nb vertices =45 , nb triangles = 64 , nb boundary edges 24 -- Build Nodes/DF on mesh : n.v. 1305, n. elmt. 5376, n b. elmt. 1472 nb of Nodes 24280 nb of DoF 107304 DFon=0363 -- FESpace: Nb of Nodes 24280 Nb of DoF 107304 Error Umfpack -1 : out_of_memory current line = 75 Exec error : Error Umfpack -1 : out_of_memory -- number :1 catch an erreur in solve => set sol = 0 !!!!!!! Exec error : Error Umfpack -1 : out_of_memory -- number :1 err code 8 , mpirank 0 try getConsole C:\Users\Ming\Downloads\script (20).edp