diff --git a/GetDDM/Helmholtz.pro b/GetDDM/Helmholtz.pro index 5f3357109a6d2acb3008177fdff1d1af6d0a47ab..cda0145360fc5d550941840abd264246785bcc1d 100644 --- a/GetDDM/Helmholtz.pro +++ b/GetDDM/Helmholtz.pro @@ -124,7 +124,7 @@ FunctionSpace { } } } - If (TC_TYPE == 2) + If (TC_TYPE == 2 || TC_TYPE == 4) For h In {1:NP_OSRC} { Name Hgrad_phi~{h}~{i}~{j}; Type Form0; BasisFunction { @@ -148,7 +148,7 @@ Formulation { { Name u~{i}; Type Local; NameOfSpace Hgrad_u~{i}; } For jj In {0: #myD~{i}()-1} j = myD~{i}(jj); - If(TC_TYPE == 2) + If(TC_TYPE == 2 || TC_TYPE == 4) For h In{1:NP_OSRC} { Name phi~{h}~{i}~{j}; Type Local; NameOfSpace Hgrad_phi~{h}~{i}~{j}; } @@ -226,6 +226,28 @@ Formulation { EndFor EndIf + If(TC_TYPE == 4) // RCot + For jj In {0: #myD~{i}()-1} + j = myD~{i}(jj); + Galerkin { [ 1/Ld~{j} * Dof{u~{i}}, {u~{i}} ]; + In Sigma~{i}~{j}; Jacobian JSur; Integration I1; } + + For h In{1:NP_OSRC} + Galerkin { [ +2/Ld~{j} * Dof{phi~{h}~{i}~{j}}, {u~{i}} ]; + In Sigma~{i}~{j}; Jacobian JSur; Integration I1; } + Galerkin { [ -2/Ld~{j} * (1/k[])^2 * Dof{d phi~{h}~{i}~{j}}, {d u~{i}} ]; + In Sigma~{i}~{j}; Jacobian JSur; Integration I1; } + + Galerkin { [ -1 * Dof{u~{i}}, {phi~{h}~{i}~{j}} ]; + In Sigma~{i}~{j}; Jacobian JSur; Integration I1; } + Galerkin { [ (1-((h*Pi)/(k[]*Ld~{j}))^2) * Dof{phi~{h}~{i}~{j}}, {phi~{h}~{i}~{j}} ]; + In Sigma~{i}~{j}; Jacobian JSur; Integration I1; } + Galerkin { [ -(1/k[])^2 * Dof{d phi~{h}~{i}~{j}}, {d phi~{h}~{i}~{j}} ]; + In Sigma~{i}~{j}; Jacobian JSur; Integration I1; } + EndFor + EndFor + EndIf + // Bayliss-Turkel absorbing boundary condition Galerkin { [ - I[] * kInf[] * Dof{u~{i}} , {u~{i}} ]; In GammaInf~{i}; Jacobian JSur; Integration I1; } @@ -246,7 +268,7 @@ Formulation { { Name u~{i}; Type Local; NameOfSpace Hgrad_u~{i}; } { Name g_out~{i}~{j}; Type Local; NameOfSpace Hgrad_g_out~{i}~{j}; } - If(TC_TYPE == 2) + If(TC_TYPE == 2 || TC_TYPE == 4) For h In{1:NP_OSRC} { Name phi~{h}~{i}~{j}; Type Local; NameOfSpace Hgrad_phi~{h}~{i}~{j}; } @@ -290,6 +312,18 @@ Formulation { Galerkin { [ 2 * I[] * kInf[] * {u~{i}}, {g_out~{i}~{j}}]; In TrBndPmlSigma~{i}~{j}; Jacobian JSur; Integration I1;} EndIf + + If(TC_TYPE == 4) // RCot + Galerkin { [ -2 * 1/Ld~{j} * {u~{i}}, {g_out~{i}~{j}} ]; + In Sigma~{i}~{j}; Jacobian JSur; Integration I1; } + For h In{1:NP_OSRC} + Galerkin { [ -2 * 2/Ld~{j} * {u~{i}}, {g_out~{i}~{j}} ]; + In Sigma~{i}~{j}; Jacobian JSur; Integration I1; } + Galerkin { [ -2 * 2/Ld~{j} * ((h*Pi)/(k[]*Ld~{j}))^2 * {phi~{h}~{i}~{j}}, {g_out~{i}~{j}} ]; + In Sigma~{i}~{j}; Jacobian JSur; Integration I1; } + EndFor + EndIf + } } diff --git a/GetDDM/cavity2d.pro b/GetDDM/cavity2d.pro new file mode 100644 index 0000000000000000000000000000000000000000..3509f3e5859121249bfcfb3602909ed17cae5814 --- /dev/null +++ b/GetDDM/cavity2d.pro @@ -0,0 +1,133 @@ +Include "cavity2d_data.geo"; + +DefineConstant[ + TC_TYPE = {2, Name "Input/01Transmission condition", + Choices {0="Order 0", 1="Order 2", 2="Pade (OSRC)", 4="RCot"}}, + NP_OSRC = {4, Name "Input/02Number of auxiliary variables (OSRC or RCot)"}, + PRECONDITIONER = 0, + IS_OPEN = {0, Name "Input/03Is cavity open?", Choices{0, 1}}, + MODE = {-1, Name "Input/04Mode"}, + ListOfCuts() = { {0, N_DOM-1} } +]; + +Include "cavity2d_src.pro"; + +Function{ + // Imaginary value + I[] = Complex[0, 1]; + + // Domain length + Ld~{0} = LX * t2Dom; + Ld~{1} = LX * (1-t2Dom); + + // Wavenumber + k = WAVENUMBER; + k[] = k; + + // ABC + kInf[] = k; + alphaBT[] = 0; + betaBT[] = 0; + + // 0th order TC + kIBC[] = k[]; + + // 2nd order TC (Gander 2002, pp. 46-47) + xsimin = 0; + xsimax = Pi / LC; + deltak[] = Pi; + alphastar[] = I[] * ((k[]^2 - xsimin^2) * (k[]^2 - (k[]-deltak[])^2))^(1/4); + betastar[] = ((xsimax^2 - k[]^2) * ((k[]+deltak[])^2 - k[]^2))^(1/4); + a[] = - (alphastar[] * betastar[] - k[]^2) / (alphastar[] + betastar[]); + b[] = - 1 / (alphastar[] + betastar[]); + + // Pade TC + kepsI = 0; + keps[] = k[]*(1+kepsI*I[]); + theta_branch = Pi/4; + + // No PMLs + D[] = 1; + E[] = 1; +} + +Group{ + D() = {}; + For idom In {0:N_DOM-1} + left = (idom-1)%N_DOM; // left neighbour + right = (idom+1)%N_DOM; // right neighbour + + D() += idom; + + Omega~{idom} = Region[(idom + 1)]; + GammaD~{idom} = Region[{}]; + GammaD0~{idom} = Region[(idom + 2*N_DOM + 2)]; + GammaInf~{idom} = Region[{}]; + GammaN~{idom} = Region[{}]; + + If(idom == 0) + D~{idom} = {right}; + + Sigma~{idom}~{right} = Region[(idom + N_DOM + 2)]; + Sigma~{idom} = Region[{Sigma~{idom}~{right}}]; + GammaD~{idom} += Region[(idom + N_DOM + 1)]; + + BndSigma~{idom}~{right} = Region[{}]; + BndSigma~{idom} = Region[{BndSigma~{idom}~{right}}]; + BndGammaInf~{idom}~{right} = Region[{}]; + BndGammaInf~{idom} = Region[{BndGammaInf~{idom}~{right}}]; + + Pml~{idom}~{right} = Region[{}]; + PmlD0~{idom}~{right} = Region[{}]; + PmlInf~{idom}~{right} = Region[{}]; + EndIf + If(idom == N_DOM-1) + D~{idom} = {left}; + + Sigma~{idom}~{left} = Region[(idom + N_DOM + 1)]; + Sigma~{idom} = Region[{Sigma~{idom}~{left}}]; + + If(IS_OPEN) + GammaInf~{idom} += Region[(idom + N_DOM + 2)]; + Else + GammaD0~{idom} += Region[(idom + N_DOM + 2)]; + EndIf + BndSigma~{idom}~{left} = Region[{}]; + BndSigma~{idom} = Region[{BndSigma~{idom}~{left}}]; + BndGammaInf~{idom}~{left} = Region[{}]; + BndGammaInf~{idom} = Region[{BndGammaInf~{idom}~{left}}]; + + Pml~{idom}~{left} = Region[{}]; + PmlD0~{idom}~{left} = Region[{}]; + PmlInf~{idom}~{left} = Region[{}]; + EndIf + If(idom >= 1 && idom < N_DOM-1) + D~{idom} = {left, right}; + + Sigma~{idom}~{left} = Region[(idom + N_DOM + 1)]; + Sigma~{idom}~{right} = Region[(idom + N_DOM + 2)]; + Sigma~{idom} = Region[{Sigma~{idom}~{left}, + Sigma~{idom}~{right}}]; + + BndSigma~{idom}~{left} = Region[{}]; + BndSigma~{idom}~{right} = Region[{}]; + BndSigma~{idom} = Region[{BndSigma~{idom}~{left}, + BndSigma~{idom}~{right}}]; + BndGammaInf~{idom}~{left} = Region[{}]; + BndGammaInf~{idom}~{right} = Region[{}]; + BndGammaInf~{idom} = Region[{BndGammaInf~{idom}~{left}, + BndGammaInf~{idom}~{right}}]; + + Pml~{idom}~{left} = Region[{}]; + Pml~{idom}~{right} = Region[{}]; + PmlD0~{idom}~{left} = Region[{}]; + PmlD0~{idom}~{right} = Region[{}]; + PmlInf~{idom}~{left} = Region[{}]; + PmlInf~{idom}~{right} = Region[{}]; + EndIf + EndFor +} + +Include "Decomposition.pro"; +Include "Helmholtz.pro"; +Include "Schwarz.pro"; diff --git a/GetDDM/cavity2d_src.pro b/GetDDM/cavity2d_src.pro new file mode 100644 index 0000000000000000000000000000000000000000..8cb4a599c75a732a6e2e79e5860453eb10ea04b7 --- /dev/null +++ b/GetDDM/cavity2d_src.pro @@ -0,0 +1,115 @@ +Function{ + If(MODE > 0) + ky = MODE * Pi / LY; + uinc[] = Complex[Sin[ky * Y[]], 0]; + Else + uinc[] = Complex[Sin[ 1 * Pi/LY * Y[]], 0] + + Complex[Sin[ 2 * Pi/LY * Y[]], 0] + + Complex[Sin[ 3 * Pi/LY * Y[]], 0] + + Complex[Sin[ 4 * Pi/LY * Y[]], 0] + + Complex[Sin[ 5 * Pi/LY * Y[]], 0] + + Complex[Sin[ 6 * Pi/LY * Y[]], 0] + + Complex[Sin[ 7 * Pi/LY * Y[]], 0] + + Complex[Sin[ 8 * Pi/LY * Y[]], 0] + + Complex[Sin[ 9 * Pi/LY * Y[]], 0] + + + Complex[Sin[10 * Pi/LY * Y[]], 0] + + Complex[Sin[11 * Pi/LY * Y[]], 0] + + Complex[Sin[12 * Pi/LY * Y[]], 0] + + Complex[Sin[13 * Pi/LY * Y[]], 0] + + Complex[Sin[14 * Pi/LY * Y[]], 0] + + Complex[Sin[15 * Pi/LY * Y[]], 0] + + Complex[Sin[16 * Pi/LY * Y[]], 0] + + Complex[Sin[17 * Pi/LY * Y[]], 0] + + Complex[Sin[18 * Pi/LY * Y[]], 0] + + Complex[Sin[19 * Pi/LY * Y[]], 0] + + + Complex[Sin[20 * Pi/LY * Y[]], 0] + + Complex[Sin[21 * Pi/LY * Y[]], 0] + + Complex[Sin[22 * Pi/LY * Y[]], 0] + + Complex[Sin[23 * Pi/LY * Y[]], 0] + + Complex[Sin[24 * Pi/LY * Y[]], 0] + + Complex[Sin[25 * Pi/LY * Y[]], 0] + + Complex[Sin[26 * Pi/LY * Y[]], 0] + + Complex[Sin[27 * Pi/LY * Y[]], 0] + + Complex[Sin[28 * Pi/LY * Y[]], 0] + + Complex[Sin[29 * Pi/LY * Y[]], 0] + + + Complex[Sin[30 * Pi/LY * Y[]], 0] + + Complex[Sin[31 * Pi/LY * Y[]], 0] + + Complex[Sin[32 * Pi/LY * Y[]], 0] + + Complex[Sin[33 * Pi/LY * Y[]], 0] + + Complex[Sin[34 * Pi/LY * Y[]], 0] + + Complex[Sin[35 * Pi/LY * Y[]], 0] + + Complex[Sin[36 * Pi/LY * Y[]], 0] + + Complex[Sin[37 * Pi/LY * Y[]], 0] + + Complex[Sin[38 * Pi/LY * Y[]], 0] + + Complex[Sin[39 * Pi/LY * Y[]], 0] + + + Complex[Sin[40 * Pi/LY * Y[]], 0] + + Complex[Sin[41 * Pi/LY * Y[]], 0] + + Complex[Sin[42 * Pi/LY * Y[]], 0] + + Complex[Sin[43 * Pi/LY * Y[]], 0] + + Complex[Sin[44 * Pi/LY * Y[]], 0] + + Complex[Sin[45 * Pi/LY * Y[]], 0] + + Complex[Sin[46 * Pi/LY * Y[]], 0] + + Complex[Sin[47 * Pi/LY * Y[]], 0] + + Complex[Sin[48 * Pi/LY * Y[]], 0] + + Complex[Sin[49 * Pi/LY * Y[]], 0] + + + Complex[Sin[50 * Pi/LY * Y[]], 0] + + Complex[Sin[51 * Pi/LY * Y[]], 0] + + Complex[Sin[52 * Pi/LY * Y[]], 0] + + Complex[Sin[53 * Pi/LY * Y[]], 0] + + Complex[Sin[54 * Pi/LY * Y[]], 0] + + Complex[Sin[55 * Pi/LY * Y[]], 0] + + Complex[Sin[56 * Pi/LY * Y[]], 0] + + Complex[Sin[57 * Pi/LY * Y[]], 0] + + Complex[Sin[58 * Pi/LY * Y[]], 0] + + Complex[Sin[59 * Pi/LY * Y[]], 0] + + + Complex[Sin[60 * Pi/LY * Y[]], 0] + + Complex[Sin[61 * Pi/LY * Y[]], 0] + + Complex[Sin[62 * Pi/LY * Y[]], 0] + + Complex[Sin[63 * Pi/LY * Y[]], 0] + + Complex[Sin[64 * Pi/LY * Y[]], 0] + + Complex[Sin[65 * Pi/LY * Y[]], 0] + + Complex[Sin[66 * Pi/LY * Y[]], 0] + + Complex[Sin[67 * Pi/LY * Y[]], 0] + + Complex[Sin[68 * Pi/LY * Y[]], 0] + + Complex[Sin[69 * Pi/LY * Y[]], 0] + + + Complex[Sin[70 * Pi/LY * Y[]], 0] + + Complex[Sin[71 * Pi/LY * Y[]], 0] + + Complex[Sin[72 * Pi/LY * Y[]], 0] + + Complex[Sin[73 * Pi/LY * Y[]], 0] + + Complex[Sin[74 * Pi/LY * Y[]], 0] + + Complex[Sin[75 * Pi/LY * Y[]], 0] + + Complex[Sin[76 * Pi/LY * Y[]], 0] + + Complex[Sin[77 * Pi/LY * Y[]], 0] + + Complex[Sin[78 * Pi/LY * Y[]], 0] + + Complex[Sin[79 * Pi/LY * Y[]], 0] + + + Complex[Sin[80 * Pi/LY * Y[]], 0] + + Complex[Sin[81 * Pi/LY * Y[]], 0] + + Complex[Sin[82 * Pi/LY * Y[]], 0] + + Complex[Sin[83 * Pi/LY * Y[]], 0] + + Complex[Sin[84 * Pi/LY * Y[]], 0] + + Complex[Sin[85 * Pi/LY * Y[]], 0] + + Complex[Sin[86 * Pi/LY * Y[]], 0] + + Complex[Sin[87 * Pi/LY * Y[]], 0] + + Complex[Sin[88 * Pi/LY * Y[]], 0] + + Complex[Sin[89 * Pi/LY * Y[]], 0] + + + Complex[Sin[90 * Pi/LY * Y[]], 0] + + Complex[Sin[91 * Pi/LY * Y[]], 0] + + Complex[Sin[92 * Pi/LY * Y[]], 0] + + Complex[Sin[93 * Pi/LY * Y[]], 0] + + Complex[Sin[94 * Pi/LY * Y[]], 0] + + Complex[Sin[95 * Pi/LY * Y[]], 0] + + Complex[Sin[96 * Pi/LY * Y[]], 0] + + Complex[Sin[97 * Pi/LY * Y[]], 0] + + Complex[Sin[98 * Pi/LY * Y[]], 0] + + Complex[Sin[99 * Pi/LY * Y[]], 0]; + EndIf +}