//======================================================== // Benchmark "3D unit cell" // File: Equations // Contributors: // Miguel Rojas S. // Pablo Stuardo O. //======================================================== Include "unitcell.dat"; Group { Port1 = Region[MasterA]; Port2 = Region[SlaveA]; Port3 = Region[MasterB]; Port4 = Region[SlaveB]; Contour = Region[PEC]; Vol = Region[Vol]; Wall = Region[{Contour, Port1, Port2, Port3, Port4}]; Vol_S = Region[{Wall}] ; Vol_C = Region[{Vol}] ; Tot = Region[{Vol_S,Vol_C}] ; } DefineConstant[ Excitation = {1,Name "2Signal/1Mode Type", Choices{1="TM01", 2="TM02", 3="TM03", 4="TE01", 5="TE02", 6="TE03", 7="TM11", 8="TM12", 9="TM13", 10="TE11", 11="TE12", 12="TE13"}}, NMODES = {1, Min 1, Max 10, Step 0.1, Name "3Parameters/1Numbers of Modes to find"}, FREQ = { 1.2e+9, Min 1e+6, Max 10e+9, Step 1, Name "3Parameters/2Frequency Shift[Hz]"} ]; Function { nu [] = 1./mu0 ; epsilon [] = ep0 ; sigma [] = 0; nu_cmplx[] = Complex[0,nu[]*(1/($EigenvalueReal))]; ALPHA_MIN = -(ALPHA/2); ALPHA_MAX = (ALPHA/2); ALPHA_DIV = Pi/180; R_MIN = 0.95; R_MAX = 1.05; R_DIV = 0.01; H_MIN = 0; H_MAX = H; H_DIV = H/50; X_01 = 2.40482555769577276862; // TM01 Mode X_02 = 5.52007811028631064959; // TM02 Mode X_03 = 8.65372791291101221695; // TM03 Mode X_11 = 3.83170597020751231561; // TM11 Mode X_12 = 7.01558666981561875353; // TM12 Mode X_13 = 10.1734681350627220771; // TM13 Mode X_21 = 5.13562230184068255630; // TM21 Mode Xp_01 = 3.8317059702075123156; // TE01 Mode Xp_02 = 7.0155866698156187535; // TE02 Mode Xp_03 = 10.173468135062722077; // TE03 Mode Xp_11 = 1.8412; // TE11 Mode Xp_12 = 5.3314; // TE12 Mode Xp_13 = 8.5363; // TE13 Mode Xp_21 = 3.0542; // TE21 Mode Xp_31 = 4.2012; // TE31 Mode Xp_41 = 5.3176; // TE41 Mode k = (2*Pi*FREQ)/c0; If (Excitation == 1 ) kc = X_01/R; ElseIf (Excitation == 2) kc = X_02/R; ElseIf (Excitation == 3) kc = X_03/R; ElseIf (Excitation == 4) kc = Xp_01/R; ElseIf (Excitation == 5) kc = Xp_02/R; ElseIf (Excitation == 6) kc = Xp_03/R; ElseIf (Excitation == 7) kc = X_11/R; ElseIf (Excitation == 8) kc = X_12/R; ElseIf (Excitation == 9) kc = X_13/R; ElseIf (Excitation == 10) kc = Xp_11/R; ElseIf (Excitation == 11) kc = Xp_12/R; ElseIf (Excitation == 12) kc = Xp_13/R; EndIf beta[] = Sqrt[(k)^2 - (kc)^2]; gam = Sqrt[(k)^2-(kc)^2]; MPI_Printf("---------frequency: %g", FREQ); MPI_Printf("---------Kc: %g", kc); MPI_Printf("---------Ko: %g", k); MPI_Printf("---------Gam: %g", gam); } DefineConstant[ FREQ_CUT = { kc*c0/(2*Pi), ReadOnly 1, Highlight "LightGrey", Name "3Parameters/3Cutoff frequency [Hz]"}, LAMB = { c0/FREQ*100, ReadOnly 1, Highlight "LightGrey", Name "3Parameters/4Wavelength [cm]"} ]; Constraint { { Name PEC ; Case { { Region Contour; Value 0. ; } } } { Name periodic ; Case { { Region Contour ; Value 0.; } { Region Port2; Type LinkCplx; RegionRef Port1; Coefficient Complex[Cos[gam*H],-Sin[gam*H]]; Function Vector[$X,$Y,$Z-H]; } // Here we have the problem, we need to create constrains for every axis in the same plane { Region Port3; Type Link; RegionRef Port4; Coefficient Vector[1,-1,1]; Function Vector[$X,$Y-2*$Y,$Z]; } } } } /* Standard volume Jacobian */ Jacobian { { Name JVol; Case { { Region Vol_C; Jacobian Vol; } { Region Vol_S; Jacobian Sur; } } } } /* Standard integration */ Integration { { Name I1 ; Case { { Type Gauss ; Case { { GeoElement Triangle ; NumberOfPoints 3; } { GeoElement Tetrahedron ; NumberOfPoints 5; } { GeoElement Pyramid ; NumberOfPoints 8; } { GeoElement Prism ; NumberOfPoints 9; } } } } } } FunctionSpace { { Name Hcurl_e_3D; Type Form1; BasisFunction { { Name se; NameOfCoef e; Function BF_Edge; Support Tot; Entity EdgesOf[All]; } } Constraint { { NameOfCoef e; EntityType EdgesOf ; NameOfConstraint periodic; } } } } Formulation { { Name MW_e_3D; Type FemEquation; Quantity { { Name e; Type Local; NameOfSpace Hcurl_e_3D; } } Equation { Galerkin { DtDt [ epsilon[] * Dof{e} , {e} ]; In Vol_C; Integration I1; Jacobian JVol; } Galerkin { [ nu[] * Dof{d e} , {d e} ]; In Vol_C; Integration I1; Jacobian JVol; } } } } Resolution { { Name MW_e_3D; System { If (FREQ_CUT <= FREQ) { Name A; NameOfFormulation MW_e_3D; Type ComplexValue; } ElseIf (FREQ_CUT > FREQ) MPI_Printf(" *** ERROR: Fc > FREQ ***"); EndIf } Operation { CreateDir["res"] ; GenerateSeparate [A]; EigenSolve[A,NMODES,(2*Pi*FREQ)^2,0]; SaveSolutions[A]; } } } PostProcessing { { Name MW_e_3D; NameOfFormulation MW_e_3D; NameOfSystem A; Quantity { // ********************* E Field ************************ { Name e; Value{ Local{ [ {e} ] ; In Vol; } } } { Name e2; Value{ Local{ [ SquNorm[{e}] ] ; In Vol; } } } // ********************* H Field ************************ { Name h; Value{ Local{ [ nu_cmplx[]*{Curl e} ] ; In Vol; } } } { Name h2; Value{ Local{ [ SquNorm[nu_cmplx[]*{Curl e}] ] ; In Vol; } } } // ****************************************************** } } } PostOperation { { Name MW_e_3D; NameOfPostProcessing MW_e_3D; Operation { // ********************* E Field Result ************************ Print[ e, OnElementsOf Wall, EigenvalueLegend, Format TimeTable, File "res/EigenValues.txt"]; Print[ e, OnElementsOf Vol , File "res\e.pos" ] ; Print[ e2, OnElementsOf Vol , File "res\e2.pos" ] ; Print[ e, OnGrid { (R-Q)*$A*Cos[$B], (R-Q)*$A*Sin[$B], $C} { 1, ALPHA_MIN:ALPHA_MAX:ALPHA_DIV, H_MIN:H_MAX:H_DIV}, File "res\e_S.pos" ] ; Print[ e, OnGrid { (R-Q)*$A*Cos[$B], (R-Q)*$A*Sin[$B], $C} { 1, ALPHA_MIN:ALPHA_MAX:ALPHA_DIV, H_MIN:H_MAX:H_DIV }, Format SimpleTable, File "res\e_S.txt" ] ; // ********************* H Field Result ************************ Print[ h, OnElementsOf Vol , File "res\h.pos" ] ; Print[ h2, OnElementsOf Vol , File "res\h2.pos" ] ; Print[ h, OnGrid { (R-Q)*$A*Cos[$B], (R-Q)*$A*Sin[$B], $C} { 1, ALPHA_MIN:ALPHA_MAX:ALPHA_DIV, H_MIN:H_MAX:H_DIV}, File "res\h_S.pos" ] ; Print[ h, OnGrid { (R-Q)*$A*Cos[$B], (R-Q)*$A*Sin[$B], $C} { 1, ALPHA_MIN:ALPHA_MAX:ALPHA_DIV, H_MIN:H_MAX:H_DIV}, Format SimpleTable, File "res\h_S.txt" ] ; // ************************************************************* } } }