Skip to content
Snippets Groups Projects
Commit ac0efbe3 authored by Bertrand Thierry's avatar Bertrand Thierry
Browse files

Helmholtz Model : Corner Condition for ABC and DDM

parent 4e631e9e
No related branches found
No related tags found
No related merge requests found
DefineConstant[
L = {10, Min 10, Max 1000, Step 0.1, Name "Input/05Geometry/10X-width"},
l = {10, Min 5, Max 100, Step 0.1, Name "Input/05Geometry/11Y-width"},
TC_BROKEN_LINE = {1, Choices{1="Broken Line", 0="Straight Line"}, Name "Input/05Geometry/15Type of border line"},
dx = {1, Min 0.1, Max L/4, Step 0.1, Name "Input/05Geometry/X-coord. of the pick point", Visible TC_BROKEN_LINE},
dy = {2, Min 0.1, Max l/2-0.1, Step 0.1, Name "Input/05Geometry/Y-coord. of bottom point", Visible TC_BROKEN_LINE},
INCOMING_CONDITION = {"Fourier", Choices {"Dirichlet", "Fourier"}, Name "Input/10Boundary Conditions/1Incoming (left)"},
OUTGOING_CONDITION = {"Fourier", Choices {"Dirichlet", "Neumann", "Fourier"}, Name "Input/10Boundary Conditions/2Outgoing (right, =0)"},
TOP_CONDITION = {"Neumann", Choices {"Dirichlet", "Neumann"}, Name "Input/10Boundary Conditions/3Top (=0)"},
BOTTOM_CONDITION = {"Neumann", Choices {"Dirichlet", "Neumann"}, Name "Input/10Boundary Conditions/4Bottom (=0)"},
kap = {Pi, Min 1, Max 30, Step 0.1, Name "Input/15Input/0Wavenumber"},
INCIDENT_WAVE = {"Plane Wave", Choices {"Eigenmode", "Plane Wave"}, Name "Input/15Input/2Type of incident Wave"},
alpha = {0, Min 0, Max 2*Pi, Step 0.1, Name "Input/15Input/3Incident angle", Visible (!StrCmp(INCIDENT_WAVE,"Plane Wave"))},
m_mode = {1, Min 0, Max 1000, Step 1, Name "Input/15Input/3Mode number", Visible (!StrCmp(INCIDENT_WAVE, "Eigenmode"))},
NLambda = {15, Min 10, Max 100, Step 1, Name "Input/Mesh/NLambda"},
h = {2*Pi/NLambda/kap, ReadOnly 1, Name "Input/Mesh/h"},
TC_ORDER = {2, Choices {0="Ordre 0", 2="Ordre 2"}, Name "Input/01DDM/05TC order"},
CORNER_CONDITION = {(!TC_BROKEN_LINE)?0:2, Choices {0 = "Dirichlet", 1="Homog. Neumann", 2="Corner Correction"}, Name "Input/01DDM/10Corner Condition", ReadOnly (TC_ORDER == 0 || !TC_BROKEN_LINE)},
NDom = {2, Choices {2}, ReadOnly 1, Visible 0, Name "Input/01DDM/15N sub domain"},
TC_DIRREL = {"out_TC/", Name "Input/Output/0Directory"},
TC_DIR = {StrCat[WorkingDir,TC_DIRREL], Visible 0, ReadOnly 1, Name "Input/Output/Output Directory (Absolute)"},
TC_PREFIX = {"", Name "Input/Output/05Prefix for filename"},
PRINT_PHI = {0, Choices {0,1}, Name "Input/Output/10Print Every Phi" },
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"}
];
//Rectangle separe en deux avec un coin
//4 sommets
Point(1) = {L, l/2, 0,h};
Point(2) = {0, l/2, 0,h};
Point(3) = {0, -l/2,0,h};
Point(4) = {L, -l/2,0,h};
dx_geo = TC_BROKEN_LINE?dx:0;
dy_geo = TC_BROKEN_LINE?dy:l/2-1;
//5 points interface
Point(11) = {L/2, l/2, 0,h};
Point(12) = {L/2, dy_geo,0,h};
Point(13) = {L/2+dx_geo,0, 0,h}; //pointe
Point(14) = {L/2,-dy_geo, 0,h};
Point(15) = {L/2,-l/2, 0,h};
// Point{13} In Surface{1}; //force pointe a etre noeud du maillage
Line(2) = {2,3}; //segments bords rectangle
Line(4) = {4,1}; //segments bords rectangle
Line(21) = {11,2}; // Gamma Top inter. Omega_1
Line(24) = {3,15}; // Gamma Bottom inter. Omega_1
Line(22) = {1,11}; // Gamma Top inter. Omega_2
Line(25) = {15,4}; // Gamma Bottom inter. Omega_2
Line(30) = {11,12}; //interface Sigma
Line(31) = {12,13};
Line(32) = {13,14};
Line(33) = {14,15};
Curve Loop(1) = {-33,-32,-31,-30,21,2,24}; //bord Omega1
Plane Surface(1) = {1}; //Omega1
Curve Loop(2) = {4,22,30,31,32,33,25}; //bord Omega2
Plane Surface(2) = {2}; //Omega2
Physical Surface(0) = {1}; //Omega1
Physical Surface(1) = {2}; //Omega2
Physical Curve(31) = {2}; // Left Boundary (incoming)
Physical Curve(33) = {4}; // Right Boundary (outcoming)
Physical Curve(100) = {21}; // Gamma_Top on Omega_1
Physical Curve(101) = {22}; // Gamma_Top on Omega_2
Physical Curve(200) = {24}; // Gamma_Bottom on Omega_1
Physical Curve(201) = {25}; // Gamma_Bottom on Omega_2
For i In {0:3}
Physical Point(40+i) = {i+1}; //noeuds Gamma
Physical Point(50+i) = {10+i+1}; //noeuds Sigma
Physical Curve(20+i) = {30+i}; //interface Sigma
EndFor
Physical Point(54) = {15};
// Due to TC_main.geo particular case, Mesh mush not be automatically computed
Printf(OnelabAction);
Solver.AutoMesh = -1; // the .geo file creates the mesh
If(StrCmp(OnelabAction, "compute") == 0)
Mesh 2;
Save StrCat[WorkingDir, "main.msh"];
CreateDir StrCat[WorkingDir, "aux"];
Plugin(BoundaryAngles).Remove = 1;
Plugin(BoundaryAngles).Save = 1;
Plugin(BoundaryAngles).Dir = "aux";
Plugin(BoundaryAngles).Filename = "TC_Theta";
Plugin(BoundaryAngles).Run;
EndIf
Include "TC_main.data";
If(!StrCmp(TC_PREFIX,""))
TC_PREFIX_ = "";
Else
TC_PREFIX_ = StrCat[TC_PREFIX, "_"];
EndIf
Function {
I[] = Complex[0.,1.]; //i complexe
alpha_vect[] = Vector[Cos[alpha],Sin[alpha],0]; //vecteur angle alpha
If(!StrCmp(INCIDENT_WAVE,"Plane Wave"))
//one plane incidente vect alpha
uInc[] = Complex[Cos[kap*(XYZ[]*alpha_vect[])],Sin[kap*(XYZ[]*alpha_vect[])]];
dx_uInc[] = I[]*kap*Cos[alpha]*uInc[];
ElseIf(!StrCmp(INCIDENT_WAVE,"Eigenmode"))
uInc[] = Cos[m_mode*Pi*(Y[]+l/2)/l];
dx_uInc[] = I[]*kap*uInc[];
EndIf
}
// 2 subdomains (DDM) case)
Group {
D() = {};
For i In {0:(NDom-1)}
Omega~{i} = Region[{i}];
D() += i;
j = (i+1)%NDom; // neighbor
D~{i} = {j};
If (i==0)
GammaIncoming~{i} = Region[{31}];
GammaOutgoing~{i} = Region[{}];
Else
GammaOutgoing~{i} = Region[{33}];
GammaIncoming~{i} = Region[{}];
EndIf
GammaTop~{i} = Region[{(100+i)}];
GammaBottom~{i} = Region[{(200+i)}];
//Sigma_ij = dOmega1 cap dOmega2
Sigma~{i}~{j} = Region[{20,21,22,23}];
Delta~{i}~{j} = Region[{20,21,22,23}];
Sigma~{i} = Region[{20,21,22,23}];
Delta~{i} = Region[{20,21,22,23}];
If(CORNER_CONDITION == 0) //Dirichlet
// Segments are considered as "only one"
r = 0; // only one "segment"
NDelta~{i}~{j} = 1;
Delta~{i}~{j}~{r} = Region[{20,21,22,23}];
dDelta~{i}~{j}~{r} = Region[{}];
ListdDelta~{i}~{j}~{r}() = {}; // No corner point (empty list)
ListdDelta~{i}~{j}() = {}; // No corner point (empty list)
Else
// Sigma_{ij} is split into segments ("broken" H^1(Sigma_ij))
NDelta~{i}~{j} = 4;
// if i == 0 then r = {3,2,1,0}
// else i == 1 then r = {0,1,2,3}
For rr In {0:(NDelta~{i}~{j}-1)}
r = (i==0)?(NDelta~{i}~{j}-1-rr):rr;
s = (i==0)?(rr):(NDelta~{i}~{j}-1-rr); // indice of this segment from the other side
ListdDelta~{i}~{j}~{r}() = {}; // contains the indices r' of the adjacent segments (possibly unique)
Delta~{i}~{j}~{r} = Region[{(20+rr)}];
//Nods Q^ij_r = Delta^ij_r \cap Delta^ij_r'
If (r==0) // First segment
next = r+1;
Q~{i}~{j}~{r}~{next} = Region[{((i==0)?53-r:51+r)}];
dDelta~{i}~{j}~{r}~{0} = Region[{((i==0)?53-r:51+r)}];
ListdDelta~{i}~{j}~{r}() += next; // next segment
dDelta~{i}~{j}~{r} = Region[{((i==0)?53-r:51+r)}];
ElseIf (r==1) // second segment
next = r+1;
Q~{i}~{j}~{r}~{next} = Region[{((i==0)?53-r:51+r)}];
dDelta~{i}~{j}~{r}~{0} = Region[{((i==0)?53-r:51+r)}];
ListdDelta~{i}~{j}~{r}() += next; // next segment
prev = r-1;
Q~{i}~{j}~{r}~{prev} = Region[{((i==0)?53-prev:51+prev)}];
dDelta~{i}~{j}~{r}~{1} = Region[{((i==0)?53-prev:51+prev)}];
ListdDelta~{i}~{j}~{r}() += prev; // next segment
dDelta~{i}~{j}~{r} = Region[{((i==0)?53-r:51+r), ((i==0)?53-prev:51+prev)}];
ElseIf (r==2) // third segment
next = r+1;
Q~{i}~{j}~{r}~{next} = Region[{((i==0)?53-r:51+r)}];
dDelta~{i}~{j}~{r}~{0} = Region[{((i==0)?53-r:51+r)}];
ListdDelta~{i}~{j}~{r}() += next; // next segment
prev = r-1;
Q~{i}~{j}~{r}~{prev} = Region[{((i==0)?53-prev:51+prev)}];
dDelta~{i}~{j}~{r}~{1} = Region[{((i==0)?53-prev:51+prev)}];
ListdDelta~{i}~{j}~{r}() += prev; // next segment
dDelta~{i}~{j}~{r} = Region[{((i==0)?53-r:51+r), ((i==0)?53-prev:51+prev)}];
Else // last segment
prev = r-1;
Q~{i}~{j}~{r}~{prev} = Region[{((i==0)?53-prev:51+prev)}];
dDelta~{i}~{j}~{r}~{0} = Region[{((i==0)?53-prev:51+prev)}];
ListdDelta~{i}~{j}~{r}() += prev; // next segment
dDelta~{i}~{j}~{r} = Region[{((i==0)?53-prev:51+prev)}];
EndIf
EndFor
EndIf
EndFor
}
Group{
Omega = Region[{(0:(NDom-1))}];
Sigma = Region[{20,21,22,23}];
Delta = Region[{20,21,22,23}];
GammaIncoming = Region[{31}];
GammaOutgoing = Region[{33}];
GammaTop = Region[{100, 101}];
GammaBottom = Region[{200, 201}];
}
Include "TC_Decomposition.pro";
Include "TC_WaveGuide.pro";
Include "TC_WaveGuide_DDM.pro";
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment