Select Git revision
ABC_main.pro
ABC_main.pro 12.58 KiB
DefineConstant[
//Input
alpha = {0, Min 0, Max 2*Pi, Step 0.1, Name "Input/Input/5Incident angle"},
// ABC
ABC_ORDER = {2, Choices {0="Ordre 0", 2="Ordre 2"}, Name "Input/ABC/0ABC order condition"},
CORNER_CONDITION = {2, Choices {1 ="Homog. Neumann", 2="Corner Correction"}, Name "Input/ABC/3Corner Condition", ReadOnly (ABC_ORDER == 0)},
// DDM
beta = {-kap*Cos[theta], ReadOnly 1, Name "Input/DDM/8beta parameter"},
// Solver
SOLVER = {"jacobi", Choices {"jacobi", "gmres", "print"}, Name "Input/IterativeSolver/0Solver"},
TOLlog10 = {-6, Max -1, Min -16, Step -1, Name "Input/IterativeSolver/1Tolerance (log10)"},
TOL = 10^(TOLlog10),
MAXIT = {10, Min 10, Max 10000, Step 1, Name "Input/IterativeSolver/2Max it"},
RESTART = {0, Min 0, Max 1000, Step 1, Name "Input/IterativeSolver/3Restart"},
//Output
ABC_DIRREL = {"out_ABC/", Name "Input/Output/01Output Directory"},
ABC_DIR = {StrCat[WorkingDir, ABC_DIRREL], ReadOnly 1, Visible 0, Name "Input/Output/Output Directory (Absolute)"}
ABC_PREFIX = {"", Name "Input/Output/05Prefix for filename"}
];
Function {
I[] = Complex[0.,1.]; // Pure imaginary "i"
alpha_vect[] = Vector[Cos[alpha],Sin[alpha],0]; // incidence vector
uInc[] = Complex[Cos[kap*(XYZ[]*alpha_vect[])],Sin[kap*(XYZ[]*alpha_vect[])]]; // Incident plane wave
//Exact solution for a cylindrical obstacle located at (0,0,0) of radius = rs and incident angle = alpha
uex[] = AcousticFieldSoftCylinder[XYZ[]]{kap, rs, alpha};
}
// Mono domain case
Group{
Omega = Region[{(0:(NDom-1))}];
GammaExt = Region[{30:35}];
Sigma = Region[{20:25}];
dOmegaInt = Region[{10:15}];
dOmegaExt = Region[{GammaExt}];
ListGammaExt() = {}; // indices k of the segments
// Fourier or Neumann condition
For k In {0:NGammaExt-1}
ListGammaExt() +=k;
GammaExt~{k} = Region[(30+k)];
// End-points
l_left = (k-1+NGammaExt)%NGammaExt;
l_right = (k+1)%NGammaExt;
ListdGammaExt~{k}() = {l_left,l_right};
Apt~{k}~{l_left} = Region[(40 + k)];
Apt~{k}~{l_right} = Region[(40 + (k+1)%6)];
// Every end-points of Gamma_k
dGammaExt~{k} = Region[{Apt~{k}~{l_left},Apt~{k}~{l_right}}];
EndFor
dGammaExt = Region[{40:45}]; //every end-points
TrdOmegaInt = ElementsOf[Omega, OnOneSideOf dOmegaInt];
}
//DDM Case
If (NDom == 2)
Group {
D() = {};
For i In {0:(NDom-1)}
Omega~{i} = Region[{i}];
D() += i;
// left = right = {1, 0}
j = (i+1)%NDom; // The other subdomain
D~{i} = {j};
GammaInt~{i} = Region[{(10+i)}];
NGammaExt~{i} = (NGammaExt/NDom); // Number of segments Gamma^i_k = dOmega_i cap Gamma
ListGammaExt~{i}() = {}; //Set of indices k such that Gamma^i_k is not empty
// Every subdomain has 3 Gamma^i_k :
// if i == 0 then k = {0,1,2}
// else i == 1 then k = {3,4,5} For kk In {0:(NGammaExt~{i}-1)} // kk = 0,1,2
k = kk + 3*i; // k = 0,1,2 or 3,4,5
ListGammaExt~{i}() += k;
ListdGammaExt~{i}~{k}() = {}; // set of indices (j,l) of every end-points of Gamma^i_k
GammaExt~{i}~{k} = Region[{(30+k)}]; //Gamma^i_k
// Subdomain number of neightbor segment
If(kk == 0)
j_left = j;
Else
j_left = i;
EndIf
If(kk == 2)
j_right = j;
Else
j_right = i;
EndIf
// Indice of the previous/next segment (independant of subdomain number)
l_left = (k-1+NGammaExt)%NGammaExt;
l_right = (k+1)%NGammaExt;
// first end-point A^ij_kl of Gamma^i_k
Apt~{i}~{j_left}~{k}~{l_left} = Region[{(40+k)}];
dGammaExt~{i}~{j_left}~{k}~{l_left} = Region[{(40+k)}];
ListdGammaExt~{i}~{k}() += j_left; // j
ListdGammaExt~{i}~{k}() += l_left; // l
//second end-point A^ij_kl of Gamma^i_k
Apt~{i}~{j_right}~{k}~{l_right} = Region[{(40+(k+1)%NGammaExt)}];
dGammaExt~{i}~{j_right}~{k}~{l_right} = Region[{(40+(k+1)%NGammaExt)}];
ListdGammaExt~{i}~{k}() += j_right; // j
ListdGammaExt~{i}~{k}() += l_right; // l
//Set C^i_k of every end-points of Gamma^i_k
Cset~{i}~{k} = Region[{(40+k),(40+(k+1)%NGammaExt)}];
dGammaExt~{i}~{k} = Region[{(40+k),(40+(k+1)%NGammaExt)}];
EndFor
dOmegaInt~{i}= Region[{GammaInt~{i}}];
dOmegaExt~{i}= Region[{GammaExt~{i}~{3*i}, GammaExt~{i}~{1 + 3*i}, GammaExt~{i}~{2 + 3*i}}];
//Sigma_{i,j} = dOmega_i cap dOmega_j
Sigma~{i}~{j} = Region[{20}];
//2 end points on Sigma cap GammaExt
ListdSigmaExt~{i}~{j}() = {};
//First
k = 0 + 3*i;
l_left = (k - 1 + NGammaExt)%NGammaExt;
dSigmaExt~{i}~{j}~{0} = Region[{Apt~{i}~{j}~{k}~{l_left}}];
ListdSigmaExt~{i}~{j}() += j;
ListdSigmaExt~{i}~{j}() += l_left;
//Second
k = 2 + 3*i;
l_right = (k +1 )%NGammaExt;
dSigmaExt~{i}~{j}~{1} = Region[{Apt~{i}~{j}~{k}~{l_right}}];
ListdSigmaExt~{i}~{j}() += j;
ListdSigmaExt~{i}~{j}() += l_right;
// union
dSigmaExt~{i}~{j} = Region[{dSigmaExt~{i}~{j}~{0}, dSigmaExt~{i}~{j}~{1}}];
//2 end points on GammaInt
l = j;
ListdSigmaInt~{i}~{j}() = {};
dSigmaInt~{i}~{j}~{0} = Region[{(50)}];
ListdSigmaInt~{i}~{j}() += i; // k == i (6 sub dom)
ListdSigmaInt~{i}~{j}() += l; // l == j (6 sub dom)
dSigmaInt~{i}~{j}~{1} = Region[{(53)}];
ListdSigmaInt~{i}~{j}() += i; // k == i (6 sub dom)
ListdSigmaInt~{i}~{j}() += l; // l == j (6 sub dom)
dSigmaInt~{i}~{j} = Region[{dSigmaInt~{i}~{j}~{0}, dSigmaInt~{i}~{j}~{1}}];
dSigma~{i}~{j} = Region[{dSigmaExt~{i}~{j},dSigmaInt~{i}~{j}}];
Sigma~{i} = Region[{Sigma~{i}~{j}}];
EndFor
}
EndIf
If (NDom==3)
Group {
D() = {};
For i In {0:(NDom-1)}
Omega~{i} = Region[{i}];
D() += i;
left = (i-1+NDom)%NDom; //previous subdomain (counter-clowise orientation)
right = (i+1)%NDom; //next subdomain (counter-clowise orientation)
D~{i} = {left, right}; // neighbors subdomains
dOmegaInt~{i} = Region[{(10+i)}];
NGammaExt~{i} = (NGammaExt/NDom); // Number of Gamma^i_k (= 2)
ListGammaExt~{i}() = {}; //Indices k of these Gamma^i_k
For kk In {0:(NGammaExt~{i}-1)} // kk = 0,1
k = kk + 2*i; // k = [0,1] or [2,3] or [4,5]
ListGammaExt~{i}() += k;
ListdGammaExt~{i}~{k}() = {}; // indices (j,l) of the neightbor Gamma^i_k
GammaExt~{i}~{k} = Region[{(30+k)}]; //Gamma^i_k
// Subdomain number of neightbor segment
If(kk == 0)
j_left = left;
j_right = i;
EndIf
If(kk == 1)
j_left = i;
j_right = right;
EndIf
l_left = (k-1+NGammaExt)%NGammaExt;
l_right = (k+1)%NGammaExt;
//first end-point A^ij_kl of Gamma^i_k
Apt~{i}~{j_left}~{k}~{l_left} = Region[{(40+k)}];
ListdGammaExt~{i}~{k}() += j_left; // j
ListdGammaExt~{i}~{k}() += l_left; // l
//second noeud A^ij_kl de Gamma^i_k
Apt~{i}~{j_right}~{k}~{l_right} = Region[{(40+(k+1)%NGammaExt)}];
ListdGammaExt~{i}~{k}() += j_right; // j
ListdGammaExt~{i}~{k}() += l_right; // l
//ensemble C^i_k = noeuds bord Gamma^i_k
Cset~{i}~{k} = Region[{(40+k),(40+(k+1)%NGammaExt)}];
dGammaExt~{i}~{k} = Region[{(40+k),(40+(k+1)%NGammaExt)}];
EndFor
// Union
dOmegaExt~{i} = Region[{GammaExt~{i}~{2*i}, GammaExt~{i}~{1 + 2*i}}];
// Transmission boundaries
Sigma~{i}~{left} = Region[{(20 + i)}];
Sigma~{i}~{right} = Region[{(20 + right)}];
Sigma~{i} = Region[{Sigma~{i}~{left}, Sigma~{i}~{right}}];
// Unique End-point of Sigma_{ij} that also belons to Gamma
ListdSigmaExt~{i}~{left}() = {};
ListdSigmaExt~{i}~{right}() = {};
k_left = 2*i;
l_left = (k_left - 1 +NGammaExt )%NGammaExt;
k_right = 2*i +1 ;
l_right = (k_right + 1 )%NGammaExt;
dSigmaExt~{i}~{left}~{0} = Region[{Apt~{i}~{left}~{k_left}~{l_left}}];
ListdSigmaExt~{i}~{left}() += k_left;
ListdSigmaExt~{i}~{left}() += l_left;
dSigmaExt~{i}~{left} = Region[{dSigmaExt~{i}~{left}~{0}}];
dSigmaExt~{i}~{right}~{0} = Region[{Apt~{i}~{right}~{k_right}~{l_right}}];
ListdSigmaExt~{i}~{right}() += k_right;
ListdSigmaExt~{i}~{right}() += l_right;
dSigmaExt~{i}~{right} = Region[{dSigmaExt~{i}~{right}~{0}}];
ListdSigmaInt~{i}~{left}() = {};
dSigmaInt~{i}~{left}~{0} = Region[{(50 + 2*i)}];
ListdSigmaInt~{i}~{left}() += i; // Useless ?
ListdSigmaInt~{i}~{left}() += left; // Useless ?
dSigmaInt~{i}~{left} = Region[{dSigmaInt~{i}~{left}~{0}}];
ListdSigmaInt~{i}~{right}() = {};
dSigmaInt~{i}~{right}~{0} = Region[{(50 + 2*i + 1)}];
ListdSigmaInt~{i}~{right}() += i; // Useless ?
ListdSigmaInt~{i}~{right}() += right; // Useless ?
dSigmaInt~{i}~{right} = Region[{dSigmaInt~{i}~{right}~{0}}];
dSigma~{i}~{left} = Region[{dSigmaExt~{i}~{left},dSigmaInt~{i}~{left}}];
dSigma~{i}~{right} = Region[{dSigmaExt~{i}~{right},dSigmaInt~{i}~{right}}];
//Sigma i = Sigma i i-1 U Sigma i i+1
Sigma~{i} = Region[{Sigma~{i}~{left},Sigma~{i}~{right}}];
dSigmaExt~{i} = Region[{dSigmaExt~{i}~{left},dSigmaExt~{i}~{right}}];
dSigmaInt~{i} = Region[{dSigmaInt~{i}~{left},dSigmaInt~{i}~{right}}];
dSigma~{i} = Region[{dSigmaExt~{i},dSigmaInt~{i}}];
EndFor
}
EndIf
If (NDom==6)
Group {
D() = {};
For i In {0:(NDom-1)}
Omega~{i} = Region[{i}]; //Omega_i
D() += i;
left = (i-1+NDom)%NDom; // previous subdomain (counter-clockwise orientation)
right = (i+1)%NDom; //next subdomain (counter-clockwise orientation)
D~{i} = {left,right};
NGammaExt~{i} = (NGammaExt/NDom); // Number of segments Gamma^i_k (=1)
dOmegaInt~{i} = Region[{(10+i)}]; //Gamma i int
ListGammaExt~{i}() = {}; // index k of this Gamma^i_k and in this case k=i
k = i;
ListGammaExt~{i}() += k;
GammaExt~{i}~{k} = Region[{(30+k)}]; //Gamma i k ext
ListdGammaExt~{i}~{k}() = {}; // list of indices (j,l) for end-points A^{ij}_{kl} of Gamma^i_k
j_left = left;
l_left = left; // l = j
j_right = right;
l_right = right; //l =j
// first end-point A^ij_kl of Gamma^i_k
Apt~{i}~{j_left}~{k}~{l_left} = Region[{(40+k)}];
dGammaExt~{i}~{j_left}~{k}~{l_left} = Region[{(40+k)}];
ListdGammaExt~{i}~{k}() += j_left; // j
ListdGammaExt~{i}~{k}() += l_left; // l
//second end-point A^ij_kl of Gamma^i_k
Apt~{i}~{j_right}~{k}~{l_right} = Region[{(40+(k+1)%NGammaExt)}];
dGammaExt~{i}~{j_right}~{k}~{l_right} = Region[{(40+(k+1)%NGammaExt)}];
ListdGammaExt~{i}~{k}() += j_right; // j
ListdGammaExt~{i}~{k}() += l_right; // l
//Set C^i_k = end-points of Gamma^i_k
Cset~{i}~{k} = Region[{(40+k),(40+(k+1)%NGammaExt)}];
dGammaExt~{i}~{k} = Region[{(40+k),(40+(k+1)%NGammaExt)}];
// Union (here =) of gamma^i_k
dOmegaExt~{i}= Region[{GammaExt~{i}~{k}}];
//Sigma_{ij} = dOmega_i cap dOmega_j
// For the left subdomain
Sigma~{i}~{left} = Region[{(20+i)}]; ListdSigmaExt~{i}~{left}() = {}; // get k and l for provided i and j
dSigmaExt~{i}~{left}~{0} = Region[{Apt~{i}~{left}~{i}~{left}}];
ListdSigmaExt~{i}~{left}() += k; // k == i (6 sub dom)
ListdSigmaExt~{i}~{left}() += l_left; // l == left (6 sub dom)
dSigmaExt~{i}~{left} = Region[{dSigmaExt~{i}~{left}~{0}}];
ListdSigmaInt~{i}~{left}() = {}; // get k and l for provided i and j
dSigmaInt~{i}~{left}~{0} = Region[{(50+i)}];
ListdSigmaInt~{i}~{left}() += i; // k == i (6 sub dom)
ListdSigmaInt~{i}~{left}() += left; // l == left (6 sub dom)
dSigmaInt~{i}~{left} = Region[{dSigmaInt~{i}~{left}~{0}}];
dSigma~{i}~{left} = Region[{dSigmaExt~{i}~{left},dSigmaInt~{i}~{left}}];
// For the right subdomain
Sigma~{i}~{right} = Region[{(20+right)}];
ListdSigmaExt~{i}~{right}() = {}; // get k and l for provided i and j
dSigmaExt~{i}~{right}~{0} = Region[{Apt~{i}~{right}~{i}~{right}}];
ListdSigmaExt~{i}~{right}() += i; // k == i (6 sub dom)
ListdSigmaExt~{i}~{right}() += right; // l == right (6 sub dom)
dSigmaExt~{i}~{right} = Region[{dSigmaExt~{i}~{right}~{0}}];
ListdSigmaInt~{i}~{right}() = {}; // get k and l for provided i and jerer k et l pour i et j donnes
dSigmaInt~{i}~{right}~{0} = Region[{(50+(i+1)%NDom)}];
ListdSigmaInt~{i}~{right}() += i; // k == i (6 sub dom)
ListdSigmaInt~{i}~{right}() += right; // l == right (6 sub dom)
dSigmaInt~{i}~{right} = Region[{dSigmaInt~{i}~{right}~{0}}];
dSigma~{i}~{right} = Region[{dSigmaExt~{i}~{right},dSigmaInt~{i}~{right}}];
//Sigma i = Sigma_{i, left} U Sigma_{i,right}
Sigma~{i} = Region[{Sigma~{i}~{left},Sigma~{i}~{right}}];
dSigmaExt~{i} = Region[{dSigmaExt~{i}~{left},dSigmaExt~{i}~{right}}];
dSigmaInt~{i} = Region[{dSigmaInt~{i}~{left},dSigmaInt~{i}~{right}}];
dSigma~{i} = Region[{dSigmaExt~{i},dSigmaInt~{i}}];
EndFor
}
EndIf
//===========================
Include "ABC_MonoDomain.pro";
Include "ABC_DDM.pro";
//===========================