Skip to content
Snippets Groups Projects
Commit 68ffc945 authored by Guillaume Demesy's avatar Guillaume Demesy
Browse files

Revert "Merge branch 'master' of https://gitlab.onelab.info/doc/models"

This reverts commit 7fb2fd83, reversing
changes made to 83c0a56a.
parent 7fb2fd83
No related branches found
No related tags found
No related merge requests found
Pipeline #5805 passed
Showing
with 0 additions and 3871 deletions
Include "domDisk.dat";
//==================================================================================================
// OPTIONS and PARAMETERS
//==================================================================================================
DefineConstant[
nPade = {4, Choices {0, 1, 2, 3, 4, 5, 6},
Name "Input/5Model/05Pade: Number of fields",
Visible (BND_TYPE == BND_Pade)},
thetaPadeInput = {3, Choices {0, 1, 2, 3, 4},
Name "Input/5Model/06Pade: Rotation of branch cut",
Visible (BND_TYPE == BND_Pade)}
];
If(BND_TYPE == BND_Pade)
If(thetaPadeInput == 0)
thetaPade = 0;
ElseIf(thetaPadeInput == 1)
thetaPade = Pi/8;
ElseIf(thetaPadeInput == 2)
thetaPade = Pi/4;
ElseIf(thetaPadeInput == 3)
thetaPade = Pi/3;
ElseIf(thetaPadeInput == 4)
thetaPade = Pi/2;
EndIf
mPade = 2*nPade+1;
For j In{1:nPade}
cPade~{j} = Tan[j*Pi/mPade]^2;
EndFor
Else
nPade = 0;
thetaPadeInput = 0;
EndIf
If(BND_TYPE == BND_PML)
nPade = Npml;
EndIf
Group {
Dom = Region[{DOM}];
BndSca = Region[{BND_Scatt}];
If(BND_TYPE == BND_PML)
DomPml = Region[{DOM_PML}];
BndExt = Region[{BND_PmlExt}];
Else
DomPml = Region[{}];
BndExt = Region[{BND_Dom}];
EndIf
DomAll = Region[{Dom,DomPml,BndSca,BndExt}];
}
Function {
If(BND_TYPE == BND_Pade)
kEps[] = WAVENUMBER + I[] * 0.39 * WAVENUMBER^(1/3) * (1/R_DOM)^(2/3);
EndIf
If(BND_TYPE == BND_Pade)
ExpPTheta[] = Complex[Cos[ thetaPade],Sin[ thetaPade]];
ExpMTheta[] = Complex[Cos[-thetaPade],Sin[-thetaPade]];
ExpPTheta2[] = Complex[Cos[thetaPade/2.],Sin[thetaPade/2.]];
For i In{1:nPade}
For j In{1:nPade}
coefA~{i}~{j}[] = 2./mPade * cPade~{j} * (cPade~{i}-1+ExpMTheta[]) / (cPade~{i}+cPade~{j}+ExpMTheta[]);
coefB~{i}~{j}[] = 2./mPade * cPade~{j} * (-1-cPade~{i}) / (cPade~{i}+cPade~{j}+ExpMTheta[]);
EndFor
EndFor
EndIf
If(BND_TYPE == BND_PML)
rLoc[DomPml] = R[]-R_DOM;
absFuncS[DomPml] = 1/(Lpml-rLoc[]);
absFuncF[DomPml] = -Log[1-rLoc[]/Lpml];
//absFuncS[DomPml] = 1/(Lpml-rLoc[]) - 1/Lpml;
//absFuncF[DomPml] = -Log[1-rLoc[]/Lpml] - rLoc[]/Lpml;
If(rotPml < 91)
rot[DomPml] = Complex[Sin[rotPml*Pi/180.], Cos[rotPml*Pi/180.]]; // I (rotPml=0, prop) - 1 (rotPml=Pi/2, evan)
Else
rot[DomPml] = Complex[1., 1.];
EndIf
s1[DomPml] = 1 + rot[] * absFuncS[]/k[];
s2[DomPml] = 1 + rot[] * (1/R[]) * absFuncF[]/k[];
nVec[DomPml] = XYZ[]/R[];
tVec[DomPml] = nVec[] /\ Vector[0,0,1];
nTen[DomPml] = Tensor[ CompX[nVec[]]*CompX[nVec[]], CompX[nVec[]]*CompY[nVec[]], CompX[nVec[]]*CompZ[nVec[]],
CompY[nVec[]]*CompX[nVec[]], CompY[nVec[]]*CompY[nVec[]], CompY[nVec[]]*CompZ[nVec[]],
CompZ[nVec[]]*CompX[nVec[]], CompZ[nVec[]]*CompY[nVec[]], CompZ[nVec[]]*CompZ[nVec[]] ];
tTen[DomPml] = Tensor[ CompX[tVec[]]*CompX[tVec[]], CompX[tVec[]]*CompY[tVec[]], CompX[tVec[]]*CompZ[tVec[]],
CompY[tVec[]]*CompX[tVec[]], CompY[tVec[]]*CompY[tVec[]], CompY[tVec[]]*CompZ[tVec[]],
CompZ[tVec[]]*CompX[tVec[]], CompZ[tVec[]]*CompY[tVec[]], CompZ[tVec[]]*CompZ[tVec[]] ];
pmlScal[DomPml] = s1[]*s2[];
pmlTens[DomPml] = (s2[]/s1[]) * nTen[] + (s1[]/s2[]) * tTen[];
EndIf
}
//==================================================================================================
// FONCTION SPACES with CONSTRAINTS
//==================================================================================================
If(FLAG_SIGNAL_BC == SIGNAL_Dirichlet)
Constraint {
{ Name DirichletBC; Case {
If(FLAG_SIGNAL_BC == SIGNAL_Dirichlet)
{ Region BndSca; Value -f_inc[]; }
EndIf
}}
}
EndIf
FunctionSpace {
{ Name H_ref; Type Form0; BasisFunction {{ Name si; NameOfCoef pi; Function BF_Node; Support Region[{DomAll}]; Entity NodesOf[All]; }}
If(FLAG_SIGNAL_BC == SIGNAL_Dirichlet)
Constraint {{ NameOfCoef pi; EntityType NodesOf; NameOfConstraint DirichletBC; }}
EndIf
}
{ Name H_num; Type Form0; BasisFunction {{ Name si; NameOfCoef pi; Function BF_Node; Support Region[{DomAll}]; Entity NodesOf[All]; }}
If((FLAG_SIGNAL_BC == SIGNAL_Dirichlet) || (BND_TYPE == BND_PML))
Constraint {{ NameOfCoef pi; EntityType NodesOf; NameOfConstraint DirichletBC; }}
EndIf
}
If(BND_TYPE == BND_Pade)
For i In {1:nPade}
{ Name H~{i}; Type Form0; BasisFunction {{ Name sn; NameOfCoef un; Function BF_Node; Support Region[{BndExt}]; Entity NodesOf[All]; }}}
EndFor
EndIf
}
//==================================================================================================
// FORMULATIONS
//==================================================================================================
Formulation {
{ Name NumSol; Type FemEquation;
Quantity {
{ Name u_num; Type Local; NameOfSpace H_num; }
If(BND_TYPE == BND_Pade)
For i In {1:nPade}
{ Name u~{i}; Type Local; NameOfSpace H~{i}; }
EndFor
EndIf
}
Equation {
// Helmholtz
Galerkin{ [ Dof{d u_num} , {d u_num} ]; In Dom; Jacobian JVol; Integration I1; }
Galerkin{ [ -k[]^2*Dof{u_num} , {u_num} ]; In Dom; Jacobian JVol; Integration I1; }
If(FLAG_SIGNAL_BC == SIGNAL_Neumann)
Galerkin{ [ -df_inc[] , {u_num} ]; In BndSca; Jacobian JSur; Integration I1; }
EndIf
If(BND_TYPE == BND_PML)
Galerkin{ [ pmlTens[] * Dof{d u_num} , {d u_num} ]; In DomPml; Jacobian JVol; Integration I1; }
Galerkin{ [ -k[]^2 * pmlScal[] * Dof{u_num} , {u_num} ]; In DomPml; Jacobian JVol; Integration I1; }
EndIf
// Sommerfeld ABC
If(BND_TYPE == BND_Sommerfeld)
Galerkin{ [ -I[]*k[]*Dof{u_num} , {u_num} ]; In BndExt; Jacobian JSur; Integration I1; }
EndIf
// Second-order ABC
If(BND_TYPE == BND_Second)
Galerkin { [ - I[]*k[] * Dof{u_num} , {u_num} ]; In BndExt; Jacobian JSur; Integration I1; }
Galerkin { [ - 1/(2*I[]*k[]) * Dof{d u_num} , {d u_num} ]; In BndExt; Jacobian JSur; Integration I1; }
EndIf
// HABC
If(BND_TYPE == BND_Pade)
Galerkin { [ -I[]*k[]*ExpPTheta2[] * Dof{u_num} , {u_num} ]; In BndExt; Jacobian JSur; Integration I1; }
For i In{1:nPade}
Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{i} * Dof{u~{i}} , {u_num} ]; In BndExt; Jacobian JSur; Integration I1; }
Galerkin { [ -I[]*k[]*ExpPTheta2[] * 2./mPade * cPade~{i} * Dof{u_num} , {u_num} ]; In BndExt; Jacobian JSur; Integration I1; }
Galerkin { [ Dof{d u~{i}} , {d u~{i}} ]; In BndExt; Jacobian JSur; Integration I1; }
Galerkin { [ -kEps[]^2*(ExpPTheta[]*cPade~{i}+1) * Dof{u~{i}} , {u~{i}} ]; In BndExt; Jacobian JSur; Integration I1; }
Galerkin { [ -kEps[]^2*ExpPTheta[]*(cPade~{i}+1) * Dof{u_num} , {u~{i}} ]; In BndExt; Jacobian JSur; Integration I1; }
EndFor
EndIf
}
}
{ Name ProjSol; Type FemEquation;
Quantity {
{ Name u_refProj; Type Local; NameOfSpace H_ref; }
}
Equation {
Galerkin{ [ Dof{u_refProj} , {u_refProj} ]; In Dom; Jacobian JVol; Integration I1; }
Galerkin{ [ -f_ref[] , {u_refProj} ]; In Dom; Jacobian JVol; Integration I1; }
}
}
}
//==================================================================================================
// RESOLUTION
//==================================================================================================
Resolution {
{ Name NumSol;
System {
{ Name A; NameOfFormulation NumSol; Type Complex; }
}
Operation {
Generate[A]; Solve[A]; SaveSolution[A];
}
}
{ Name ProjSol;
System {
{ Name A; NameOfFormulation ProjSol; Type Complex; }
}
Operation {
Generate[A]; Solve[A]; SaveSolution[A];
}
}
}
//==================================================================================================
// POSTPRO / POSTOP
//==================================================================================================
PostProcessing {
{ Name NumSol; NameOfFormulation NumSol;
Quantity {
// { Name param1; Value { Local { [ rLoc[] ]; In DomPml; Jacobian JVol; }}}
// { Name param2; Value { Local { [ nVec[] ]; In DomPml; Jacobian JVol; }}}
// { Name param3; Value { Local { [ tVec[] ]; In DomPml; Jacobian JVol; }}}
{ Name u_ref; Value { Local { [ f_ref[] ]; In Region[{Dom,BndSca}]; Jacobian JVol; }}}
{ Name u_num~{BND_TYPE}~{nPade}~{thetaPadeInput}; Value { Local { [ {u_num} ]; In Region[{Dom,DomPml}]; Jacobian JVol; }}}
{ Name u_err~{BND_TYPE}~{nPade}~{thetaPadeInput}; Value { Local { [ f_ref[]-{u_num} ]; In Region[{Dom,BndSca}]; Jacobian JVol; }}}
}
}
{ Name Errors; NameOfFormulation NumSol;
Quantity {
{ Name energy; Value { Integral { [ Abs[f_ref[]]^2 ]; In Dom; Jacobian JVol; Integration I1; }}}
{ Name error2; Value { Integral { [ Abs[f_ref[]-{u_num}]^2 ]; In Dom; Jacobian JVol; Integration I1; }}}
{ Name errorAbs; Value { Term { Type Global; [Sqrt[$error2Var]] ; In Dom; Jacobian JVol; }}}
{ Name errorRel; Value { Term { Type Global; [Sqrt[$error2Var/$energyVar]] ; In Dom; Jacobian JVol; }}}
}
}
{ Name ProjError; NameOfFormulation ProjSol;
Quantity {
{ Name energy; Value { Integral { [ Abs[f_ref[]]^2 ]; In Dom; Jacobian JVol; Integration I1; }}}
{ Name errorProj2; Value { Integral { [ Abs[f_ref[]-{u_refProj}]^2 ]; In Dom; Jacobian JVol; Integration I1; }}}
{ Name errorProjAbs; Value { Term { Type Global; [Sqrt[$errorProj2Var]] ; In Dom; Jacobian JVol; }}}
{ Name errorProjRel; Value { Term { Type Global; [Sqrt[$errorProj2Var/$energyVar]] ; In Dom; Jacobian JVol; }}}
}
}
{ Name DtNError; NameOfFormulation NumSol;
Quantity {
{ Name energyDtN; Value { Integral { [ Abs[f_ref[]]^2 ]; In BndSca; Jacobian JLin; Integration I1; }}}
{ Name errorDtN2; Value { Integral { [ Abs[f_ref[]-{u_num}]^2 ]; In BndSca; Jacobian JLin; Integration I1; }}}
{ Name normDtN; Value { Term { Type Global; [Sqrt[$energyDtNVar]] ; In BndSca; Jacobian JLin; }}}
{ Name errorDtNAbs; Value { Term { Type Global; [Sqrt[$errorDtN2Var]] ; In BndSca; Jacobian JLin; }}}
{ Name errorDtNRel; Value { Term { Type Global; [Sqrt[$errorDtN2Var/$energyDtNVar]] ; In BndSca; Jacobian JLin; }}}
}
}
}
PostOperation{
{ Name NumSol; NameOfPostProcessing NumSol;
Operation {
// Print [param1, OnElementsOf DomPml, File "out/param1.pos"];
// Print [param2, OnElementsOf DomPml, File "out/param2.pos"];
// Print [param3, OnElementsOf DomPml, File "out/param3.pos"];
tmp1 = Sprintf("out/solNum_%g_%g_%g_%g.pos", FLAG_SIGNAL_BC, BND_TYPE, nPade, thetaPadeInput);
tmp2 = Sprintf("out/solErr_%g_%g_%g_%g.pos", FLAG_SIGNAL_BC, BND_TYPE, nPade, thetaPadeInput);
Print [u_ref, OnElementsOf Region[{DOM,BND_Scatt}], File "out/solRef.pos"];
Print [u_num~{BND_TYPE}~{nPade}~{thetaPadeInput}, OnElementsOf Region[{Dom,DomPml}], File tmp1];
Print [u_err~{BND_TYPE}~{nPade}~{thetaPadeInput}, OnElementsOf Region[{Dom,BndSca}], File tmp2];
}
}
{ Name Errors; NameOfPostProcessing Errors;
Operation {
tmp3 = Sprintf("out/errorAbs_%g_%g_%g_%g.dat", FLAG_SIGNAL_BC, BND_TYPE, nPade, thetaPadeInput);
tmp4 = Sprintf("out/errorRel_%g_%g_%g_%g.dat", FLAG_SIGNAL_BC, BND_TYPE, nPade, thetaPadeInput);
Print [energy[Dom], OnRegion Dom, Format Table, StoreInVariable $energyVar];
Print [error2[Dom], OnRegion Dom, Format Table, StoreInVariable $error2Var];
Print [errorAbs, OnRegion Dom, Format Table, SendToServer "Output/1L2-Error (absolute)", File > tmp3];
Print [errorRel, OnRegion Dom, Format Table, SendToServer "Output/2L2-Error (relative)", File > tmp4];
}
}
{ Name ProjError; NameOfPostProcessing ProjError;
Operation {
tmp5 = Sprintf("out/errorProjAbs_%g_%g_%g_%g.dat", FLAG_SIGNAL_BC, BND_TYPE, nPade, thetaPadeInput);
tmp6 = Sprintf("out/errorProjRel_%g_%g_%g_%g.dat", FLAG_SIGNAL_BC, BND_TYPE, nPade, thetaPadeInput);
Print [energy[Dom], OnRegion Dom, Format Table, StoreInVariable $energyVar];
Print [errorProj2[Dom], OnRegion Dom, Format Table, StoreInVariable $errorProj2Var];
Print [errorProjAbs, OnRegion Dom, Format Table, SendToServer "Output/3L2-ErrorProj (absolute)", File > tmp5];
Print [errorProjRel, OnRegion Dom, Format Table, SendToServer "Output/4L2-ErrorProj (relative)", File > tmp6];
}
}
{ Name DtNError; NameOfPostProcessing DtNError;
Operation {
tmp7 = Sprintf("out/normDtNAbs_%g_%g_%g_%g.dat", FLAG_SIGNAL_BC, BND_TYPE, nPade, thetaPadeInput);
tmp8 = Sprintf("out/errorDtNAbs_%g_%g_%g_%g.dat", FLAG_SIGNAL_BC, BND_TYPE, nPade, thetaPadeInput);
tmp9 = Sprintf("out/errorDtNRel_%g_%g_%g_%g.dat", FLAG_SIGNAL_BC, BND_TYPE, nPade, thetaPadeInput);
Print [energyDtN[BndSca], OnRegion BndSca, Format Table, StoreInVariable $energyDtNVar];
Print [errorDtN2[BndSca], OnRegion BndSca, Format Table, StoreInVariable $errorDtN2Var];
Print [normDtN, OnRegion BndSca, Format Table, SendToServer "Output/4L2-DtN-Norm", File > tmp7];
Print [errorDtNAbs, OnRegion BndSca, Format Table, SendToServer "Output/5L2-DtN-Error (absolute)", File > tmp8];
Print [errorDtNRel, OnRegion BndSca, Format Table, SendToServer "Output/6L2-DtN-Error (relative)", File > tmp9];
}
}
}
DefineConstant[
R_ = {"NumSol", Name "GetDP/1ResolutionChoices", Visible 1, Choices {"NumSol", "ProjSol"}},
P_ = {"NumSol", Name "GetDP/2PostOperationChoices", Visible 1, Choices {"NumSol", "Errors", "ProjError", "DtNError"}}
];
REF_Ana = 0;
REF_Num = 1;
DefineConstant[
dimL = { 3.3, Min 0.1, Step 0.1, Max 20, Name "Input/1Geometry/1Wedge length"},
alphaDom = { 90, Min 60, Step 10 , Max 180, Name "Input/1Geometry/2Wedge angle"},
distSca = { 2.2, Min 0.1, Step 0.1, Max 10, Name "Input/1Geometry/3Scatterer distance"},
// X_SCA = { 1.1, Min 1, Step 0.01, Max 10, Name "Input/1Geometry/2Scatterer position (x0)"},
// Y_SCA = { 1.1, Min 1, Step 0.01, Max 10, Name "Input/1Geometry/3Scatterer position (y0)"}
FLAG_REF = {REF_Num,
Name "Input/1Reference solution",
Choices {REF_Ana = "Analytic solution", REF_Num = "Numeric solution"}}
];
alphaDomSave = alphaDom;
alphaDom = alphaDom*Pi/180;
alphaSca = alphaDom/2;
//dimL = (R_SCA*1.1)*(1.+1./Sin[alphaSca]);
//distSca = (R_SCA*1.1)/Sin[alphaSca];
X_SCA = distSca * Sin[alphaSca];
Y_SCA = distSca * Cos[alphaSca];
X~{1} = - X_SCA;
Y~{1} = dimL - Y_SCA;
X~{100} = Sin[alphaDom/2]*dimL - X_SCA;
Y~{100} = Cos[alphaDom/2]*dimL - Y_SCA;
X~{2} = Sin[alphaDom ]*dimL - X_SCA;
Y~{2} = Cos[alphaDom ]*dimL - Y_SCA;
X~{3} = - X_SCA;
Y~{3} = - Y_SCA;
X~{200} =-Sin[alphaDom/2]*dimL - X_SCA;
Y~{200} =-Cos[alphaDom/2]*dimL - Y_SCA;
CRN_1 = 101;
CRN_2 = 102;
CRN_3 = 103;
BND_1 = 201;
BND_2 = 202;
BND_3 = 203;
BND_Scatt = 205;
BND_EXT = 206;
DOM = 301;
DOM_EXT = 302;
Include "domPieWedge.dat";
Point(0) = {0, 0, 0};
Point(1) = {X~{1}, Y~{1}, 0};
Point(2) = {X~{2}, Y~{2}, 0};
Point(3) = {X~{3}, Y~{3}, 0};
Point(100) = {X~{100}, Y~{100}, 0};
Point(5) = {-R_SCA, 0, 0};
Point(6) = { 0,-R_SCA, 0};
Point(7) = { R_SCA, 0, 0};
Point(8) = { 0, R_SCA, 0};
Circle(0) = {2, 3, 100};
Circle(1) = {100, 3, 1};
Line(2) = {3, 2};
Line(3) = {1, 3};
Circle(5) = {5, 0, 6};
Circle(6) = {6, 0, 7};
Circle(7) = {7, 0, 8};
Circle(8) = {8, 0, 5};
Line Loop(1) = {0, 1, 2, 3, -5, -6, -7, -8};
Plane Surface(1) = {1};
If(FLAG_REF == REF_Num)
Point(200) = {X~{200}, Y~{200}, 0};
Circle(10) = {1, 3, 200};
Circle(11) = {200, 3, 2};
Line Loop(2) = {-2, -3, 10, 11};
Plane Surface(2) = {2};
EndIf
Physical Point(CRN_1) = {1}; // Hybrid1
Physical Point(CRN_2) = {2}; // Hybrid2
Physical Point(CRN_3) = {3}; // Pade
Physical Line(BND_1) = {0, 1}; // BGT
Physical Line(BND_2) = {2}; // Pade-Left
Physical Line(BND_3) = {3}; // Pade-Down
Physical Line(BND_Scatt) = {5,6,7,8};
SetOrder 1;
Mesh.ElementOrder = 1;
Mesh 1;
Save "mainCurv.msh";
SetOrder ORDER;
Mesh.ElementOrder = ORDER;
Physical Surface(DOM) = {1};
If(FLAG_REF == REF_Num)
Physical Line(BND_EXT) = {10, 11};
Physical Surface(DOM_EXT) = {2};
EndIf
This diff is collapsed.
REF_Ana = 1;
REF_Num = 2;
DefineConstant[
DOM_midR = { 1.65, Min 1, Step 0.1, Max 10, Name "Input/1Geometry/1Polygon: mid-radius"},
DOM_Num = { 3, Choices {3, 4, 5, 6, 8, 12, 18, 1}, Name "Input/1Geometry/2Polygon: number edges"},
DOM_AngleEdges = { 180.-360./DOM_Num, Name "Input/1Geometry/32Polygon: angle edges", ReadOnly 1},
FLAG_REF = {REF_Ana,
Name "Input/1Reference solution",
Choices {REF_Ana = "Analytic solution", REF_Num = "Numeric solution"}},
L_EXT = { 5, Min 1, Step 0.1, Max 10, Name "Input/2Geometry/0Reference length",
Visible (FLAG_REF == REF_Num)}
];
angleInterior = 360/DOM_Num;
angleEdges = 180-angleInterior;
DOM_R = 2./(1.+Cos[Pi/DOM_Num])*DOM_midR;
If(DOM_Num == 1)
DOM_R = DOM_midR;
EndIf
For i In{1:DOM_Num}
X~{i} = DOM_R*Cos[(i-0.5)*angleInterior*Pi/180-Pi/2];
Y~{i} = DOM_R*Sin[(i-0.5)*angleInterior*Pi/180-Pi/2];
CRN~{i} = 100+i;
BND~{i} = 200+i;
EndFor
BND_SCA = 299;
DOM = 301;
For i In{1:4}
CRN_EXT~{i} = 400+i;
BND_EXT~{i} = 500+i;
EndFor
BND_EXT = 599;
DOM_EXT = 601;
Include "domPolygon.dat";
Point(0) = {0, 0, 0};
Point(1) = {-R_SCA, 0, 0};
Point(2) = { 0,-R_SCA, 0};
Point(3) = { R_SCA, 0, 0};
Point(4) = { 0, R_SCA, 0};
Circle(1) = {1, 0, 2};
Circle(2) = {2, 0, 3};
Circle(3) = {3, 0, 4};
Circle(4) = {4, 0, 1};
llScat[] = {-1, -2, -3, -4};
If(DOM_Num == 1)
Point(11) = {-DOM_R, 0, 0};
Point(12) = { 0,-DOM_R, 0};
Point(13) = { DOM_R, 0, 0};
Point(14) = { 0, DOM_R, 0};
Circle(11) = {11, 0, 12};
Circle(12) = {12, 0, 13};
Circle(13) = {13, 0, 14};
Circle(14) = {14, 0, 11};
llBnd[] = {-11, -12, -13, -14};
Else
For iCrn In{1:DOM_Num}
pBnd~{iCrn} = newp;
Point(pBnd~{iCrn}) = {X~{iCrn}, Y~{iCrn}, 0};
EndFor
For iEdg In{1:DOM_Num}
iCrn1 = (iEdg == 1) ? DOM_Num : iEdg-1;
iCrn2 = iEdg;
lBnd~{iEdg} = newl;
Line(lBnd~{iEdg}) = {pBnd~{iCrn1}, pBnd~{iCrn2}};
llBnd[] += lBnd~{iEdg};
EndFor
EndIf
Line Loop(1) = {llBnd[], llScat[]};
Plane Surface(1) = {1};
If(DOM_Num == 1)
Physical Line(BND~{1}) = {11, 12, 13, 14};
Else
For i In{1:DOM_Num}
Physical Point(CRN~{i}) = {pBnd~{i}};
Physical Line(BND~{i}) = {lBnd~{i}};
EndFor
EndIf
SetOrder 1;
Mesh.ElementOrder = 1;
Mesh 1;
Save "mainCurv.msh";
SetOrder ORDER;
Mesh.ElementOrder = ORDER;
Physical Line(BND_SCA) = {1,2,3,4};
Physical Surface(DOM) = {1};
If(FLAG_REF == REF_Num)
pExt~{1} = newp; Point(pExt~{1}) = {-L_EXT/2,-L_EXT/2, 0};
pExt~{2} = newp; Point(pExt~{2}) = { L_EXT/2,-L_EXT/2, 0};
pExt~{3} = newp; Point(pExt~{3}) = { L_EXT/2, L_EXT/2, 0};
pExt~{4} = newp; Point(pExt~{4}) = {-L_EXT/2, L_EXT/2, 0};
lExt~{1} = newl; Line(lExt~{1}) = {pExt~{4}, pExt~{1}};
lExt~{2} = newl; Line(lExt~{2}) = {pExt~{1}, pExt~{2}};
lExt~{3} = newl; Line(lExt~{3}) = {pExt~{2}, pExt~{3}};
lExt~{4} = newl; Line(lExt~{4}) = {pExt~{3}, pExt~{4}};
For i In{1:4}
Physical Point(CRN_EXT~{i}) = {pExt~{i}};
Physical Line(BND_EXT~{i}) = {lExt~{i}};
llExt[] += lExt~{i};
EndFor
Line Loop(2) = {llBnd[], llExt[]};
Plane Surface(2) = {2};
Physical Surface(DOM_EXT) = {2};
EndIf
This diff is collapsed.
DefineConstant[
dimRmean = {2, Min 1, Step 0.1, Max 5, Name "Input/1Geometry/3Polyhedron: midradius"},
FAC_Num = {6,
Name "Input/1Geometry/2Polyhedron: number faces", Highlight "Blue",
Choices {4 = "Tetrahedron",
6 = "Cube",
8 = "Octahedron",
12 = "Dodecahedron",
20 = "Icosahedron",
1 = "Sphere"} }
];
GoldRatio = 1.61803398875;
If(FAC_Num == 4)
dimRext = dimRmean * Sin[Pi/4]/Cos[Pi/3] * Tan[Pi/3]/Sqrt[2];
dimRint = dimRext/(Tan[Pi/3]*Tan[Pi/3]);
ElseIf(FAC_Num == 6)
dimRext = dimRmean * Sin[Pi/6]/Cos[Pi/4] * Tan[Pi/3];
dimRint = dimRext/(Tan[Pi/4]*Tan[Pi/3]);
ElseIf(FAC_Num == 8)
dimRext = dimRmean * Sin[Pi/6]/Cos[Pi/3] * Tan[Pi/4]*Sqrt[2];
dimRint = dimRext/(Tan[Pi/3]*Tan[Pi/4]);
ElseIf(FAC_Num == 12)
dimRext = dimRmean * Sin[Pi/10]/Cos[Pi/5] * Tan[Pi/3] * GoldRatio;
dimRint = dimRext/(Tan[Pi/5]*Tan[Pi/3]);
ElseIf(FAC_Num == 20)
dimRext = dimRmean * Sin[Pi/10]/Cos[Pi/3] * Tan[Pi/5] * GoldRatio*GoldRatio;
dimRint = dimRext/(Tan[Pi/3]*Tan[Pi/5]);
ElseIf(FAC_Num == 1)
dimRext = dimRmean;
dimRint = dimRmean;
EndIf
Printf("... Polyhedron circumradius =", dimRext);
Printf("... Polyhedron inradius =", dimRint);
// ADJACENCY
If(FAC_Num == 4)
EDG_Num = 6;
EdgNeigh~{1}~{1} = 1; EdgNeigh~{1}~{2} = 2; // Neighboring faces
EdgNeigh~{2}~{1} = 1; EdgNeigh~{2}~{2} = 3;
EdgNeigh~{3}~{1} = 2; EdgNeigh~{3}~{2} = 3;
EdgNeigh~{4}~{1} = 1; EdgNeigh~{4}~{2} = 4;
EdgNeigh~{5}~{1} = 2; EdgNeigh~{5}~{2} = 4;
EdgNeigh~{6}~{1} = 3; EdgNeigh~{6}~{2} = 4;
CRN_Num = 4;
CRN_NumEdgPerCrn = 3;
CrnNeigh~{1}~{1} = 1; CrnNeigh~{1}~{2} = 2; CrnNeigh~{1}~{3} = 3; // Neighboring edges
CrnNeigh~{2}~{1} = 1; CrnNeigh~{2}~{2} = 4; CrnNeigh~{2}~{3} = 5;
CrnNeigh~{3}~{1} = 2; CrnNeigh~{3}~{2} = 4; CrnNeigh~{3}~{3} = 6;
CrnNeigh~{4}~{1} = 3; CrnNeigh~{4}~{2} = 5; CrnNeigh~{4}~{3} = 6;
EndIf
If(FAC_Num == 6)
EDG_Num = 12;
EdgNeigh~{1}~{1} = 3; EdgNeigh~{1}~{2} = 5;
EdgNeigh~{2}~{1} = 3; EdgNeigh~{2}~{2} = 6;
EdgNeigh~{3}~{1} = 4; EdgNeigh~{3}~{2} = 6;
EdgNeigh~{4}~{1} = 4; EdgNeigh~{4}~{2} = 5;
EdgNeigh~{5}~{1} = 1; EdgNeigh~{5}~{2} = 5;
EdgNeigh~{6}~{1} = 2; EdgNeigh~{6}~{2} = 5;
EdgNeigh~{7}~{1} = 2; EdgNeigh~{7}~{2} = 6;
EdgNeigh~{8}~{1} = 1; EdgNeigh~{8}~{2} = 6;
EdgNeigh~{9}~{1} = 1; EdgNeigh~{9}~{2} = 3;
EdgNeigh~{10}~{1} = 1; EdgNeigh~{10}~{2} = 4;
EdgNeigh~{11}~{1} = 2; EdgNeigh~{11}~{2} = 4;
EdgNeigh~{12}~{1} = 2; EdgNeigh~{12}~{2} = 3;
CRN_Num = 8;
CRN_NumEdgPerCrn = 3;
CrnNeigh~{1}~{1} = 9; CrnNeigh~{1}~{2} = 5; CrnNeigh~{1}~{3} = 1;
CrnNeigh~{2}~{1} = 12; CrnNeigh~{2}~{2} = 6; CrnNeigh~{2}~{3} = 1;
CrnNeigh~{3}~{1} = 12; CrnNeigh~{3}~{2} = 7; CrnNeigh~{3}~{3} = 2;
CrnNeigh~{4}~{1} = 9; CrnNeigh~{4}~{2} = 8; CrnNeigh~{4}~{3} = 2;
CrnNeigh~{5}~{1} = 10; CrnNeigh~{5}~{2} = 8; CrnNeigh~{5}~{3} = 3;
CrnNeigh~{6}~{1} = 11; CrnNeigh~{6}~{2} = 7; CrnNeigh~{6}~{3} = 3;
CrnNeigh~{7}~{1} = 10; CrnNeigh~{7}~{2} = 5; CrnNeigh~{7}~{3} = 4;
CrnNeigh~{8}~{1} = 11; CrnNeigh~{8}~{2} = 6; CrnNeigh~{8}~{3} = 4;
EndIf
If(FAC_Num == 8)
EDG_Num = 12;
EdgNeigh~{1}~{1} = 1; EdgNeigh~{1}~{2} = 4;
EdgNeigh~{2}~{1} = 1; EdgNeigh~{2}~{2} = 2;
EdgNeigh~{3}~{1} = 2; EdgNeigh~{3}~{2} = 3;
EdgNeigh~{4}~{1} = 3; EdgNeigh~{4}~{2} = 4;
EdgNeigh~{5}~{1} = 1; EdgNeigh~{5}~{2} = 5;
EdgNeigh~{6}~{1} = 2; EdgNeigh~{6}~{2} = 6;
EdgNeigh~{7}~{1} = 3; EdgNeigh~{7}~{2} = 7;
EdgNeigh~{8}~{1} = 4; EdgNeigh~{8}~{2} = 8;
EdgNeigh~{9}~{1} = 5; EdgNeigh~{9}~{2} = 8;
EdgNeigh~{10}~{1} = 5; EdgNeigh~{10}~{2} = 6;
EdgNeigh~{11}~{1} = 6; EdgNeigh~{11}~{2} = 7;
EdgNeigh~{12}~{1} = 7; EdgNeigh~{12}~{2} = 8;
CRN_Num = 6;
CRN_NumEdgPerCrn = 4;
CrnNeigh~{1}~{1} = 2; CrnNeigh~{1}~{2} = 1; CrnNeigh~{1}~{3} = 3; CrnNeigh~{1}~{4} = 4;
CrnNeigh~{2}~{1} = 1; CrnNeigh~{2}~{2} = 5; CrnNeigh~{2}~{3} = 8; CrnNeigh~{2}~{4} = 9;
CrnNeigh~{3}~{1} = 2; CrnNeigh~{3}~{2} = 5; CrnNeigh~{3}~{3} = 6; CrnNeigh~{3}~{4} = 10;
CrnNeigh~{4}~{1} = 3; CrnNeigh~{4}~{2} = 6; CrnNeigh~{4}~{3} = 7; CrnNeigh~{4}~{4} = 11;
CrnNeigh~{5}~{1} = 4; CrnNeigh~{5}~{2} = 7; CrnNeigh~{5}~{3} = 8; CrnNeigh~{5}~{4} = 12;
CrnNeigh~{6}~{1} = 10; CrnNeigh~{6}~{2} = 9; CrnNeigh~{6}~{3} = 11; CrnNeigh~{6}~{4} = 12;
EndIf
If(FAC_Num == 12)
EDG_Num = 30;
EdgNeigh~{1}~{1} = 1; EdgNeigh~{1}~{2} = 2;
EdgNeigh~{2}~{1} = 1; EdgNeigh~{2}~{2} = 3;
EdgNeigh~{3}~{1} = 1; EdgNeigh~{3}~{2} = 4;
EdgNeigh~{4}~{1} = 1; EdgNeigh~{4}~{2} = 5;
EdgNeigh~{5}~{1} = 1; EdgNeigh~{5}~{2} = 6;
EdgNeigh~{6}~{1} = 2; EdgNeigh~{6}~{2} = 6;
EdgNeigh~{7}~{1} = 2; EdgNeigh~{7}~{2} = 3;
EdgNeigh~{8}~{1} = 3; EdgNeigh~{8}~{2} = 4;
EdgNeigh~{9}~{1} = 4; EdgNeigh~{9}~{2} = 5;
EdgNeigh~{10}~{1} = 5; EdgNeigh~{10}~{2} = 6;
EdgNeigh~{11}~{1} = 6; EdgNeigh~{11}~{2} = 7;
EdgNeigh~{12}~{1} = 2; EdgNeigh~{12}~{2} = 7;
EdgNeigh~{13}~{1} = 2; EdgNeigh~{13}~{2} = 8;
EdgNeigh~{14}~{1} = 3; EdgNeigh~{14}~{2} = 8;
EdgNeigh~{15}~{1} = 3; EdgNeigh~{15}~{2} = 9;
EdgNeigh~{16}~{1} = 4; EdgNeigh~{16}~{2} = 9;
EdgNeigh~{17}~{1} = 4; EdgNeigh~{17}~{2} = 10;
EdgNeigh~{18}~{1} = 5; EdgNeigh~{18}~{2} = 10;
EdgNeigh~{19}~{1} = 5; EdgNeigh~{19}~{2} = 11;
EdgNeigh~{20}~{1} = 6; EdgNeigh~{20}~{2} = 11;
EdgNeigh~{21}~{1} = 7; EdgNeigh~{21}~{2} = 8;
EdgNeigh~{22}~{1} = 8; EdgNeigh~{22}~{2} = 9;
EdgNeigh~{23}~{1} = 9; EdgNeigh~{23}~{2} = 10;
EdgNeigh~{24}~{1} = 10; EdgNeigh~{24}~{2} = 11;
EdgNeigh~{25}~{1} = 7; EdgNeigh~{25}~{2} = 11;
EdgNeigh~{26}~{1} = 8; EdgNeigh~{26}~{2} = 12;
EdgNeigh~{27}~{1} = 9; EdgNeigh~{27}~{2} = 12;
EdgNeigh~{28}~{1} = 10; EdgNeigh~{28}~{2} = 12;
EdgNeigh~{29}~{1} = 11; EdgNeigh~{29}~{2} = 12;
EdgNeigh~{30}~{1} = 7; EdgNeigh~{30}~{2} = 12;
CRN_Num = 20;
CRN_NumEdgPerCrn = 3;
CrnNeigh~{1}~{1} = 1; CrnNeigh~{1}~{2} = 5; CrnNeigh~{1}~{3} = 6;
CrnNeigh~{2}~{1} = 1; CrnNeigh~{2}~{2} = 2; CrnNeigh~{2}~{3} = 7;
CrnNeigh~{3}~{1} = 2; CrnNeigh~{3}~{2} = 3; CrnNeigh~{3}~{3} = 8;
CrnNeigh~{4}~{1} = 3; CrnNeigh~{4}~{2} = 4; CrnNeigh~{4}~{3} = 9;
CrnNeigh~{5}~{1} = 4; CrnNeigh~{5}~{2} = 5; CrnNeigh~{5}~{3} = 10;
CrnNeigh~{6}~{1} = 6; CrnNeigh~{6}~{2} = 12; CrnNeigh~{6}~{3} = 11;
CrnNeigh~{7}~{1} = 7; CrnNeigh~{7}~{2} = 13; CrnNeigh~{7}~{3} = 14;
CrnNeigh~{8}~{1} = 8; CrnNeigh~{8}~{2} = 15; CrnNeigh~{8}~{3} = 16;
CrnNeigh~{9}~{1} = 9; CrnNeigh~{9}~{2} = 17; CrnNeigh~{9}~{3} = 18;
CrnNeigh~{10}~{1} = 10; CrnNeigh~{10}~{2} = 19; CrnNeigh~{10}~{3} = 20;
CrnNeigh~{11}~{1} = 12; CrnNeigh~{11}~{2} = 13; CrnNeigh~{11}~{3} = 21;
CrnNeigh~{12}~{1} = 14; CrnNeigh~{12}~{2} = 15; CrnNeigh~{12}~{3} = 22;
CrnNeigh~{13}~{1} = 16; CrnNeigh~{13}~{2} = 17; CrnNeigh~{13}~{3} = 23;
CrnNeigh~{14}~{1} = 18; CrnNeigh~{14}~{2} = 19; CrnNeigh~{14}~{3} = 24;
CrnNeigh~{15}~{1} = 11; CrnNeigh~{15}~{2} = 20; CrnNeigh~{15}~{3} = 25;
CrnNeigh~{16}~{1} = 21; CrnNeigh~{16}~{2} = 30; CrnNeigh~{16}~{3} = 26;
CrnNeigh~{17}~{1} = 22; CrnNeigh~{17}~{2} = 26; CrnNeigh~{17}~{3} = 27;
CrnNeigh~{18}~{1} = 23; CrnNeigh~{18}~{2} = 27; CrnNeigh~{18}~{3} = 28;
CrnNeigh~{19}~{1} = 24; CrnNeigh~{19}~{2} = 28; CrnNeigh~{19}~{3} = 29;
CrnNeigh~{20}~{1} = 25; CrnNeigh~{20}~{2} = 30; CrnNeigh~{20}~{3} = 29;
EndIf
If(FAC_Num == 20)
EDG_Num = 30;
EdgNeigh~{1}~{1} = 1; EdgNeigh~{1}~{2} = 5;
EdgNeigh~{2}~{1} = 1; EdgNeigh~{2}~{2} = 2;
EdgNeigh~{3}~{1} = 2; EdgNeigh~{3}~{2} = 3;
EdgNeigh~{4}~{1} = 3; EdgNeigh~{4}~{2} = 4;
EdgNeigh~{5}~{1} = 4; EdgNeigh~{5}~{2} = 5;
EdgNeigh~{6}~{1} = 1; EdgNeigh~{6}~{2} = 6;
EdgNeigh~{7}~{1} = 2; EdgNeigh~{7}~{2} = 7;
EdgNeigh~{8}~{1} = 3; EdgNeigh~{8}~{2} = 8;
EdgNeigh~{9}~{1} = 4; EdgNeigh~{9}~{2} = 9;
EdgNeigh~{10}~{1} = 5; EdgNeigh~{10}~{2} = 10;
EdgNeigh~{11}~{1} = 6; EdgNeigh~{11}~{2} = 15;
EdgNeigh~{12}~{1} = 6; EdgNeigh~{12}~{2} = 11;
EdgNeigh~{13}~{1} = 7; EdgNeigh~{13}~{2} = 11;
EdgNeigh~{14}~{1} = 7; EdgNeigh~{14}~{2} = 12;
EdgNeigh~{15}~{1} = 8; EdgNeigh~{15}~{2} = 12;
EdgNeigh~{16}~{1} = 8; EdgNeigh~{16}~{2} = 13;
EdgNeigh~{17}~{1} = 9; EdgNeigh~{17}~{2} = 13;
EdgNeigh~{18}~{1} = 9; EdgNeigh~{18}~{2} = 14;
EdgNeigh~{19}~{1} = 10; EdgNeigh~{19}~{2} = 14;
EdgNeigh~{20}~{1} = 10; EdgNeigh~{20}~{2} = 15;
EdgNeigh~{21}~{1} = 11; EdgNeigh~{21}~{2} = 16;
EdgNeigh~{22}~{1} = 12; EdgNeigh~{22}~{2} = 17;
EdgNeigh~{23}~{1} = 13; EdgNeigh~{23}~{2} = 18;
EdgNeigh~{24}~{1} = 14; EdgNeigh~{24}~{2} = 19;
EdgNeigh~{25}~{1} = 15; EdgNeigh~{25}~{2} = 20;
EdgNeigh~{26}~{1} = 16; EdgNeigh~{26}~{2} = 20;
EdgNeigh~{27}~{1} = 16; EdgNeigh~{27}~{2} = 17;
EdgNeigh~{28}~{1} = 17; EdgNeigh~{28}~{2} = 18;
EdgNeigh~{29}~{1} = 18; EdgNeigh~{29}~{2} = 19;
EdgNeigh~{30}~{1} = 19; EdgNeigh~{30}~{2} = 20;
CRN_Num = 12;
CRN_NumEdgPerCrn = 5;
CrnNeigh~{1}~{1} = 1; CrnNeigh~{1}~{2} = 2; CrnNeigh~{1}~{3} = 3; CrnNeigh~{1}~{4} = 4; CrnNeigh~{1}~{5} = 5;
CrnNeigh~{2}~{1} = 1; CrnNeigh~{2}~{2} = 6; CrnNeigh~{2}~{3} = 10; CrnNeigh~{2}~{4} = 11; CrnNeigh~{2}~{5} = 20;
CrnNeigh~{3}~{1} = 2; CrnNeigh~{3}~{2} = 6; CrnNeigh~{3}~{3} = 7; CrnNeigh~{3}~{4} = 12; CrnNeigh~{3}~{5} = 13;
CrnNeigh~{4}~{1} = 3; CrnNeigh~{4}~{2} = 7; CrnNeigh~{4}~{3} = 8; CrnNeigh~{4}~{4} = 14; CrnNeigh~{4}~{5} = 15;
CrnNeigh~{5}~{1} = 4; CrnNeigh~{5}~{2} = 8; CrnNeigh~{5}~{3} = 9; CrnNeigh~{5}~{4} = 16; CrnNeigh~{5}~{5} = 17;
CrnNeigh~{6}~{1} = 5; CrnNeigh~{6}~{2} = 9; CrnNeigh~{6}~{3} = 10; CrnNeigh~{6}~{4} = 18; CrnNeigh~{6}~{5} = 19;
CrnNeigh~{7}~{1} = 11; CrnNeigh~{7}~{2} = 12; CrnNeigh~{7}~{3} = 21; CrnNeigh~{7}~{4} = 25; CrnNeigh~{7}~{5} = 26;
CrnNeigh~{8}~{1} = 13; CrnNeigh~{8}~{2} = 14; CrnNeigh~{8}~{3} = 21; CrnNeigh~{8}~{4} = 22; CrnNeigh~{8}~{5} = 27;
CrnNeigh~{9}~{1} = 15; CrnNeigh~{9}~{2} = 16; CrnNeigh~{9}~{3} = 22; CrnNeigh~{9}~{4} = 23; CrnNeigh~{9}~{5} = 28;
CrnNeigh~{10}~{1} = 17; CrnNeigh~{10}~{2} = 18; CrnNeigh~{10}~{3} = 23; CrnNeigh~{10}~{4} = 24; CrnNeigh~{10}~{5} = 29;
CrnNeigh~{11}~{1} = 19; CrnNeigh~{11}~{2} = 20; CrnNeigh~{11}~{3} = 24; CrnNeigh~{11}~{4} = 25; CrnNeigh~{11}~{5} = 30;
CrnNeigh~{12}~{1} = 26; CrnNeigh~{12}~{2} = 27; CrnNeigh~{12}~{3} = 28; CrnNeigh~{12}~{4} = 29; CrnNeigh~{12}~{5} = 30;
EndIf
If(FAC_Num == 1)
EDG_Num = 0;
EDG_NumFacPerEdg = 0;
CRN_Num = 0;
CRN_NumEdgPerCrn = 0;
Else
EDG_NumFacPerEdg = 2;
EndIf
For i In {1:CRN_Num}
CRN~{i} = i;
EndFor
For i In {1:EDG_Num}
EDG~{i} = 100+i;
EndFor
For i In {1:FAC_Num}
SUR~{i} = 200+i;
EndFor
SUR_Scatt = 200;
VOL = 301;
This diff is collapsed.
This diff is collapsed.
BND_Neumann = 0;
BND_Sommerfeld = 1;
BND_Second = 2;
BND_Pade = 3;
BND_CRBC = 4;
BND_PML = 5;
DefineConstant[
BND_TYPE = {BND_Pade,
Name "Input/5Model/02Boundary condition (edges)", Highlight "Blue",
Choices {BND_Neumann = "Homogeneous Neumann",
BND_Sommerfeld = "Sommerfeld ABC",
BND_Second = "Second-order ABC",
BND_Pade = "Pade ABC",
BND_CRBC = "CRBC",
BND_PML = "PML"}}
];
DefineConstant[
dimL = { 2.2, Min 1, Step 0.1, Max 10, Name "Input/1Geometry/1Domain length"},
X_SCA = { 1.1, Min 1, Step 0.01, Max 10, Name "Input/1Geometry/2Scatterer position (x0)"},
Y_SCA = { 1.1, Min 1, Step 0.01, Max 10, Name "Input/1Geometry/3Scatterer position (y0)"},
Npml = { 3, Min 1, Step 1, Max 5, Name "Input/5Model/03PML: Thickness (N*Lc)", Visible (BND_TYPE == BND_PML)},
//sigmaPml = { 2e4, Choices {2e4}, Name "Input/5Model/03PML: Sigma Mult", Visible (BND_TYPE == BND_PML)},
rotPml = { 0, Min -1.57079632679, Max 1.57079632679, Step 0.0001, Name "Input/5Model/03PML: Rotation", Visible (BND_TYPE == BND_PML)}
]; // 1.57079632679
Lpml = LC*Npml;
X~{1} = -X_SCA;
Y~{1} = -Y_SCA;
X~{2} = -X_SCA+dimL;
Y~{2} = -Y_SCA;
X~{3} = -X_SCA+dimL;
Y~{3} = -Y_SCA+dimL;
X~{4} = -X_SCA;
Y~{4} = -Y_SCA+dimL;
CRN_1 = 101;
CRN_2 = 102;
CRN_3 = 103;
CRN_4 = 104;
BND_1 = 201;
BND_2 = 202;
BND_3 = 203;
BND_4 = 204;
BND_Scatt = 205;
BND_PmlExt = 206;
DOM = 301;
DOM_PML = 302;
Include "domSquare.dat";
Point(0) = {0, 0, 0};
Point(1) = {X~{1}, Y~{1}, 0};
Point(2) = {X~{2}, Y~{2}, 0};
Point(3) = {X~{3}, Y~{3}, 0};
Point(4) = {X~{4}, Y~{4}, 0};
Point(5) = {-R_SCA, 0, 0};
Point(6) = { 0,-R_SCA, 0};
Point(7) = { R_SCA, 0, 0};
Point(8) = { 0, R_SCA, 0};
Line(1) = {4, 1};
Line(2) = {1, 2};
Line(3) = {2, 3};
Line(4) = {3, 4};
Circle(5) = {5, 0, 6};
Circle(6) = {6, 0, 7};
Circle(7) = {7, 0, 8};
Circle(8) = {8, 0, 5};
Line Loop(1) = {1, 2, 3, 4, -5, -6, -7, -8};
Plane Surface(1) = {1};
If(BND_TYPE == BND_PML)
out0[] = Extrude {-Lpml,0,0} { Line{1}; Layers{Npml}; Recombine; };
Printf("==== %g ====\n", out0[0]);
Printf("==== %g ====\n", out0[1]);
Printf("==== %g ====\n", out0[2]);
Printf("==== %g ====\n", out0[3]);
out1[] = Extrude { Lpml,0,0} { Line{3}; Layers{Npml}; Recombine; };
Printf("==== %g ====\n", out1[0]);
Printf("==== %g ====\n", out1[1]);
Printf("==== %g ====\n", out1[2]);
Printf("==== %g ====\n", out1[3]);
out2[] = Extrude {0,-Lpml,0} { Line{11,2,14}; Layers{Npml}; Recombine; };
Printf("==== %g ====\n", out2[0]);
Printf("==== %g ====\n", out2[1]);
Printf("==== %g ====\n", out2[2]);
Printf("==== %g ====\n", out2[3]);
out3[] = Extrude {0, Lpml,0} { Line{10,4,15}; Layers{Npml}; Recombine; };
Printf("==== %g ====\n", out3[0]);
Printf("==== %g ====\n", out3[1]);
Printf("==== %g ====\n", out3[2]);
Printf("==== %g ====\n", out3[3]);
EndIf
Physical Point(CRN_1) = {1};
Physical Point(CRN_2) = {2};
Physical Point(CRN_3) = {3};
Physical Point(CRN_4) = {4};
Physical Line(BND_1) = {1}; // Left
Physical Line(BND_2) = {2}; // Down
Physical Line(BND_3) = {3}; // Right
Physical Line(BND_4) = {4}; // Top
Physical Line(BND_Scatt) = {5,6,7,8};
Physical Surface(DOM) = {1};
If(BND_TYPE == BND_PML)
Physical Surface(DOM_PML) = {12,16,20,24,28,32,36,40};
If(Npml == 0)
Physical Line(BND_PmlExt) = {5,6,7,8};
Else
Physical Line(BND_PmlExt) = {17,21,25,27,13,39,37,33,29,31,9,19};
EndIf
EndIf
This diff is collapsed.
DIR = "out/";
GEO_SQUARE = 1;
GEO_DISK = 2;
GEO_POLYGON = 3;
GEO_SLICE = 4;
GEO_CUBE = 5;
GEO_POLYHEDRON = 6;
SIGNAL_Harmonic = 1;
SIGNAL_Scatt = 2;
SIGNAL_Dirichlet = 1;
SIGNAL_Neumann = 2;
DefineConstant[
FLAG_GEO = {GEO_SQUARE,
Name "Input/1Geometry/0Type of domain", Highlight "Blue",
GmshOption "Reset", Autocheck 0,
Choices {GEO_SQUARE = "Squared domain (2D)",
GEO_DISK = "Circular domain (2D)",
GEO_POLYGON = "Polygonal domain (2D)",
GEO_SLICE = "Slice domain (2D)",
GEO_CUBE = "Cuboidal domain (3D)",
GEO_POLYHEDRON = "Polyhedral domain (3D)"}}
];
If ((FLAG_GEO==GEO_SQUARE) || (FLAG_GEO==GEO_DISK) || (FLAG_GEO==GEO_POLYGON) || (FLAG_GEO==GEO_SLICE))
FLAG_DIM = 2;
defWaveNum = 25;
defNLambda = 10;
ElseIf ((FLAG_GEO==GEO_CUBE) || (FLAG_GEO==GEO_POLYHEDRON))
FLAG_DIM = 3;
defWaveNum = 10;
defNLambda = 10;
EndIf
DefineConstant[
FLAG_SIGNAL = {SIGNAL_Scatt,
Name "Input/2Signal/1Type of solution",
Choices {SIGNAL_Harmonic = "Harmonic (cylindrical or spherical)",
SIGNAL_Scatt = "Scattering of plane wave"}},
FLAG_SIGNAL_BC = {SIGNAL_Neumann,
Name "Input/2Signal/2Type of condition",
Choices {SIGNAL_Dirichlet = "Sound-soft (Dirichlet BC)",
SIGNAL_Neumann = "Sound-hard (Neumann BC)"}},
SIGNAL_MODE = {0, Min 0, Step 1, Max 50, Name "Input/2Signal/3Mode number", Visible (FLAG_SIGNAL == SIGNAL_Harmonic)},
WAVENUMBER = {defWaveNum, Min 0.1, Step 0.1, Max 300, Name "Input/2Signal/4Wavenumber"},
LAMBDA = {2*Pi/WAVENUMBER, Name "Input/2Signal/5Wavelength", ReadOnly 1},
R_SCA = {1, Min 0.1, Step 0.1, Max 10, Name "Input/1Geometry/8Scatterer radius"},
N_LAMBDA = {defNLambda, Choices {5, 7.5, 10, 12.5, 15}, Name "Input/3Discretization/1Points per wavelength"},
ORDER = {2, Choices {1 = "First-order", 2 = "Second-order"}, Name "Input/3Discretization/2Polynomial order"}
];
// LAMBDA = 2.*Pi/WAVENUMBER;
// WAVENUMBER = 2.*Pi/LAMBDA;
LC = LAMBDA/N_LAMBDA;
LinkGeo = 0;
LinkPro = 0;
If (FLAG_GEO==GEO_SQUARE)
LinkGeo = "domSquare.geo";
LinkPro = "domSquare.pro";
ElseIf (FLAG_GEO==GEO_DISK)
LinkGeo = "domDisk.geo";
LinkPro = "domDisk.pro";
ElseIf (FLAG_GEO==GEO_POLYGON)
LinkGeo = "domPolygon.geo";
LinkPro = "domPolygon.pro";
ElseIf (FLAG_GEO==GEO_SLICE)
LinkGeo = "domPieWedge.geo";
LinkPro = "domPieWedge.pro";
ElseIf (FLAG_GEO==GEO_CUBE)
LinkGeo = "domCube.geo";
LinkPro = "domCube.pro";
ElseIf (FLAG_GEO==GEO_POLYHEDRON)
LinkGeo = "domPolyhedron.geo";
LinkPro = "domPolyhedron.pro";
EndIf
Include "main.dat" ;
CreateDir Str(DIR);
SetOrder ORDER;
Mesh.ElementOrder = ORDER;
Mesh.SecondOrderLinear = 0;
Mesh.CharacteristicLengthMax = LC;
Mesh.CharacteristicLengthFactor = 1;
Mesh.Optimize = 1;
// Solver.AutoMesh = 1;
Include Str[LinkGeo];
This diff is collapsed.
File deleted
This diff is collapsed.
Models developed by [[User:Ruth V. Sabariego|R. Sabariego]], J. Gyselinck. Copyright (c) 2005-2019 KU Leuven-ULg-ULB.
This diff is collapsed.
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment