Skip to content
Snippets Groups Projects
Commit dcfd65cc authored by Axel Modave's avatar Axel Modave
Browse files

AcademicWaves: update + more onelab (EM not checked)

parent b395ad20
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 476 additions and 664 deletions
AcademicWaves/formulations_scalarWaves.pro
+ 155
− 159
View file @ dcfd65cc
......@@ -9,258 +9,254 @@
//==============================================================================
Function {
DefineFunction[ fNeumann, alphaBT, betaBT, pmlScal, pmlTens ] ;
DefineVariable[ DIM, ORDER ] ;
DefineVariable[ dt, beta, gamma, theta, tf ] ;
DefineFunction[ fNeumann, alphaBT, betaBT, pmlScal, pmlTens ];
DefineVariable[ dt, beta, gamma, theta, tf, FREQ, k ];
}
Group{
DefineGroup[ Omega, OmegaPml, GammaD, GammaN, GammaR, GammaBT ] ;
SurAll = Region[{GammaD, GammaN, GammaR, GammaBT}] ;
VolAll = Region[{Omega, OmegaPml}] ;
TotAll = Region[{VolAll, SurAll}] ;
DefineGroup[ Omega, OmegaPml, GammaD, GammaN, GammaR, GammaBT ];
SurAll = Region[{GammaD, GammaN, GammaR, GammaBT}];
VolAll = Region[{Omega, OmegaPml}];
TotAll = Region[{VolAll, SurAll}];
}
Jacobian {
{ Name JVol ;
Case {
{ Region All ; Jacobian Vol ; }
}
}
{ Name JSur ;
Case {
{ Region All ; Jacobian Sur ; }
}
}
}
Integration {
{ Name I1 ;
Case {
{ Type Gauss ;
Case {
{ GeoElement Point ; NumberOfPoints 1 ; }
{ GeoElement Line ; NumberOfPoints 4 ; }
{ GeoElement Triangle ; NumberOfPoints 6 ; }
{ GeoElement Quadrangle ; NumberOfPoints 7 ; }
{ GeoElement Tetrahedron ; NumberOfPoints 15 ; }
{ GeoElement Hexahedron ; NumberOfPoints 34 ; }
}
}
}
}
}
//==============================================================================
FunctionSpace {
{ Name pSpace ; Type Form0 ;
{ Name pSpace; Type Form0;
BasisFunction {
{ Name sn ; NameOfCoef un ; Function BF_Node ;
Support TotAll ; Entity NodesOf[All] ; }
If (ORDER==2)
{ Name sn2 ; NameOfCoef un2 ; Function BF_Node_2E ;
Support TotAll ; Entity EdgesOf[All] ; }
EndIf
{ Name sn; NameOfCoef un; Function BF_Node;
Support TotAll; Entity NodesOf[All]; }
}
Constraint {
{ NameOfCoef un ; EntityType NodesOf ; NameOfConstraint pConstraint ; }
If (ORDER==2)
{ NameOfCoef un2 ; EntityType EdgesOf ; NameOfConstraint pConstraint2 ; }
EndIf
{ NameOfCoef un; EntityType NodesOf; NameOfConstraint pConstraint; }
}
}
{ Name vSpace ; Type Form0 ;
{ Name vSpace; Type Form0;
BasisFunction {
{ Name sn ; NameOfCoef vnA ; Function BF_Node ;
Support TotAll ; Entity NodesOf[All] ; }
{ Name sn; NameOfCoef vnA; Function BF_Node;
Support TotAll; Entity NodesOf[All]; }
}
Constraint {
{ NameOfCoef vnA ; EntityType NodesOf ; NameOfConstraint uConstraint ; }
{ NameOfCoef vnA; EntityType NodesOf; NameOfConstraint uConstraint; }
}
}
If (DIM==2)
{ Name uSpace ; Type Form2P ;
If (FLAG_DIM==2)
{ Name uSpace; Type Form2P;
BasisFunction {
{ Name se ; NameOfCoef he ; Function BF_PerpendicularFacet ;
Support TotAll ; Entity EdgesOf[All] ; }
{ Name se; NameOfCoef he; Function BF_PerpendicularFacet;
Support TotAll; Entity EdgesOf[All]; }
}
Constraint {
{ NameOfCoef he ; EntityType EdgesOf ; NameOfConstraint uConstraint ; }
{ NameOfCoef he; EntityType EdgesOf; NameOfConstraint uConstraint; }
}
}
EndIf
If (DIM==3)
{ Name uSpace ; Type Form2 ;
EndIf
If (FLAG_DIM==3)
{ Name uSpace; Type Form2;
BasisFunction {
{ Name se ; NameOfCoef he ; Function BF_Facet ;
Support TotAll ; Entity FacetsOf[All] ; }
{ Name se; NameOfCoef he; Function BF_Facet;
Support TotAll; Entity FacetsOf[All]; }
}
Constraint {
{ NameOfCoef he ; EntityType FacetsOf ; NameOfConstraint uConstraint ; }
{ NameOfCoef he; EntityType FacetsOf; NameOfConstraint uConstraint; }
}
}
EndIf
EndIf
}
//==============================================================================
Formulation {
{ Name WaveEquation_Form ; Type FemEquation ;
{ Name WaveEqn_Form; Type FemEquation;
Quantity {
{ Name p ; Type Local ; NameOfSpace pSpace ; }
{ Name p; Type Local; NameOfSpace pSpace; }
}
Equation {
Galerkin { DtDtDof [ Dof{p} , {p} ] ;
In VolAll ; Integration I1 ; Jacobian JVol ; }
Galerkin { [ Dof{d p} , {d p} ] ;
In Omega ; Jacobian JVol ; Integration I1 ; }
Galerkin { DtDtDof [ Dof{p} , {p} ];
In Omega; Jacobian JVol; Integration I1; }
Galerkin { [ Dof{d p} , {d p} ];
In Omega; Jacobian JVol; Integration I1; }
// BC: Neumann
Galerkin { [ -fNeumann[] , {p} ] ;
In GammaN ; Jacobian JSur ; Integration I1 ; }
Galerkin { [ -fNeumann[] , {p} ];
In GammaN; Jacobian JSur; Integration I1; }
// BC: Radiation
// Galerkin { [ -I[]*k * Dof{p} , {p} ] ;
Galerkin { DtDof [ -Dof{p} , {p} ] ;
In GammaR ; Jacobian JSur ; Integration I1 ; }
Galerkin { DtDof [ Dof{p} , {p} ];
In GammaR; Jacobian JSur; Integration I1; }
// BC: Bayliss-Turkel
// Galerkin { [ -I[]*k * Dof{p} , {p} ] ;
Galerkin { DtDof [ Dof{p} , {p} ] ;
In GammaBT ; Jacobian JSur ; Integration I1 ; }
Galerkin { [ alphaBT[] * Dof{p} , {p} ] ;
In GammaBT ; Jacobian JSur ; Integration I1 ; }
Galerkin { [ betaBT[] * Dof{d p} , {d p} ] ;
In GammaBT ; Jacobian JSur ; Integration I1 ; }
Galerkin { DtDof [ Dof{p} , {p} ];
In GammaBT; Jacobian JSur; Integration I1; }
Galerkin { [ alphaBT[] * Dof{p} , {p} ];
In GammaBT; Jacobian JSur; Integration I1; }
Galerkin { [ betaBT[] * Dof{d p} , {d p} ];
In GammaBT; Jacobian JSur; Integration I1; }
// PML (only for resolutions in the frequency domain)
Galerkin { [ pmlTens[] * Dof{d p} , {d p} ] ;
In OmegaPml ; Jacobian JVol ; Integration I1 ; }
Galerkin { [ -k^2 * (1/pmlScal[]) * Dof{p} , {p} ] ;
In OmegaPml ; Jacobian JVol ; Integration I1 ; }
Galerkin { [ pmlTens[] * Dof{d p} , {d p} ];
In OmegaPml; Jacobian JVol; Integration I1; }
Galerkin { [ -k^2 * (1/pmlScal[]) * Dof{p} , {p} ];
In OmegaPml; Jacobian JVol; Integration I1; }
}
}
{ Name MembraneEquations_Form ; Type FemEquation ;
{ Name MembraneSys_Form; Type FemEquation;
Quantity {
{ Name p ; Type Local ; NameOfSpace pSpace ; }
{ Name v ; Type Local ; NameOfSpace vSpace ; }
{ Name p; Type Local; NameOfSpace pSpace; }
{ Name v; Type Local; NameOfSpace vSpace; }
}
Equation {
Galerkin { DtDof [ Dof{p} , {p} ] ;
In VolAll ; Integration I1 ; Jacobian JVol ; }
Galerkin { [ - Dof{v} , {p} ] ;
In VolAll ; Integration I1 ; Jacobian JVol ; }
Galerkin { DtDof [ Dof{p} , {p} ];
In VolAll; Integration I1; Jacobian JVol; }
Galerkin { [ -Dof{v} , {p} ];
In VolAll; Integration I1; Jacobian JVol; }
Galerkin { DtDof [ Dof{v} , {v} ] ;
In VolAll ; Integration I1 ; Jacobian JVol ; }
Galerkin { [ Dof{d p} , {d v} ] ;
In VolAll ; Integration I1 ; Jacobian JVol ; }
Galerkin { DtDof [ Dof{v} , {v} ];
In VolAll; Integration I1; Jacobian JVol; }
Galerkin { [ Dof{d p} , {d v} ];
In VolAll; Integration I1; Jacobian JVol; }
}
}
{ Name EulerEquations_Form ; Type FemEquation ;
{ Name EulerSys_Form; Type FemEquation;
Quantity {
{ Name p ; Type Local ; NameOfSpace pSpace ; }
{ Name u ; Type Local ; NameOfSpace uSpace ; }
{ Name p; Type Local; NameOfSpace pSpace; }
{ Name u; Type Local; NameOfSpace uSpace; }
}
Equation {
Galerkin { DtDof [ Dof{p} , {p} ] ;
In VolAll ; Integration I1 ; Jacobian JVol ; }
Galerkin { [ - Dof{u} , {d p} ] ;
In VolAll ; Integration I1 ; Jacobian JVol ; }
Galerkin { DtDof [ Dof{p} , {p} ];
In VolAll; Integration I1; Jacobian JVol; }
Galerkin { [ - Dof{u} , {d p} ];
In VolAll; Integration I1; Jacobian JVol; }
Galerkin { DtDof [ Dof{u} , {u} ] ;
In VolAll ; Integration I1 ; Jacobian JVol ; }
Galerkin { [ - Dof{p} , {d u} ] ;
In VolAll ; Integration I1 ; Jacobian JVol ; }
Galerkin { DtDof [ Dof{u} , {u} ];
In VolAll; Integration I1; Jacobian JVol; }
Galerkin { [ - Dof{p} , {d u} ];
In VolAll; Integration I1; Jacobian JVol; }
// BC: Radiation
Galerkin { [ Dof{p} , {p} ] ;
In GammaR ; Integration I1 ; Jacobian JSur ; }
Galerkin { [ (Dof{u} * Normal[]) * Normal[] , {u} ] ;
In GammaR ; Integration I1 ; Jacobian JSur ; }
Galerkin { [ Dof{p} , {p} ];
In GammaR; Integration I1; Jacobian JSur; }
Galerkin { [ (Dof{u} * Normal[]) * Normal[] , {u} ];
In GammaR; Integration I1; Jacobian JSur; }
}
}
}
//==============================================================================
Resolution {
{ Name WaveEquation_FreqReso ;
If(FLAG_RES == RES_FREQ)
{ Name WaveEqn_Reso;
System {
{ Name A ; NameOfFormulation WaveEquation_Form ; Frequency Freq ; }
{ Name A; NameOfFormulation WaveEqn_Form; Frequency FREQ; }
}
Operation {
CreateDir["res/"] ;
Generate[A] ;
Solve[A] ;
SaveSolution[A] ;
CreateDir["res/"];
Generate[A]; Solve[A]; SaveSolution[A];
}
}
{ Name MembraneEquations_FreqReso ;
{ Name MembraneSys_Reso;
System {
{ Name A ; NameOfFormulation MembraneEquations_Form ; Frequency Freq ; }
{ Name A; NameOfFormulation MembraneSys_Form; Frequency FREQ; }
}
Operation {
CreateDir["res/"] ;
Generate[A] ;
Solve[A] ;
SaveSolution[A] ;
CreateDir["res/"];
Generate[A]; Solve[A]; SaveSolution[A];
}
}
{ Name EulerEquations_FreqReso ;
{ Name EulerSys_Reso;
System {
{ Name A ; NameOfFormulation EulerEquations_Form ; Frequency Freq ; }
{ Name A; NameOfFormulation EulerSys_Form; Frequency FREQ; }
}
Operation {
CreateDir["res/"] ;
Generate[A] ;
Solve[A] ;
SaveSolution[A] ;
CreateDir["res/"];
Generate[A]; Solve[A]; SaveSolution[A];
}
}
{ Name WaveEquation_TimeReso ;
EndIf
If(FLAG_RES == RES_TIME)
{ Name WaveEqn_Reso;
System {
{ Name A ; NameOfFormulation WaveEquation_Form ; }
{ Name A; NameOfFormulation WaveEqn_Form; }
}
Operation {
CreateDir["res/"] ;
InitSolution[A] ;
InitSolution[A] ;
CreateDir["res/"];
InitSolution[A];
InitSolution[A];
TimeLoopNewmark[0., tf, dt, beta, gamma] {
Generate[A] ;
Solve[A] ;
SaveSolution[A] ;
}
Generate[A]; Solve[A]; SaveSolution[A]; }
}
}
{ Name MembraneEquations_TimeReso ;
{ Name MembraneSys_Reso;
System {
{ Name A ; NameOfFormulation MembraneEquations_Form ; }
{ Name A; NameOfFormulation MembraneSys_Form; }
}
Operation {
CreateDir["res/"] ;
InitSolution[A] ;
CreateDir["res/"];
InitSolution[A];
TimeLoopTheta{
Time0 0 ; DTime dt ; Theta theta ; TimeMax tf ;
Operation {
Generate[A] ;
Solve[A] ;
SaveSolution[A] ;
}
Time0 0; DTime dt; Theta theta; TimeMax tf;
Operation { Generate[A]; Solve[A]; SaveSolution[A]; }
}
}
}
{ Name EulerEquations_TimeReso ;
{ Name EulerSys_Reso;
System {
{ Name A ; NameOfFormulation EulerEquations_Form ; }
{ Name A; NameOfFormulation EulerSys_Form; }
}
Operation {
CreateDir["res/"] ;
InitSolution[A] ;
CreateDir["res/"];
InitSolution[A];
TimeLoopTheta{
Time0 0 ; DTime dt ; Theta theta ; TimeMax tf ;
Operation {
Generate[A] ;
Solve[A] ;
SaveSolution[A] ;
}
Time0 0; DTime dt; Theta theta; TimeMax tf;
Operation { Generate[A]; Solve[A]; SaveSolution[A]; }
}
}
}
EndIf
}
//==============================================================================
PostProcessing {
{ Name WaveEqn_PostPro; NameOfFormulation WaveEqn_Form; NameOfSystem A;
Quantity {
{ Name p; Value{ Local{ [ {p} ]; In VolAll; Jacobian JVol; }}}
}
}
{ Name MembraneSys_PostPro; NameOfFormulation MembraneSys_Form; NameOfSystem A;
Quantity {
{ Name p; Value{ Local{ [ {p} ]; In VolAll; Jacobian JVol; }}}
{ Name v; Value{ Local{ [ {v} ]; In VolAll; Jacobian JVol; }}}
}
}
{ Name EulerSys_PostPro; NameOfFormulation EulerSys_Form; NameOfSystem A;
Quantity {
{ Name p; Value{ Local{ [ {p} ]; In VolAll; Jacobian JVol; }}}
{ Name u; Value{ Local{ [ {u} ]; In VolAll; Jacobian JVol; }}}
}
}
}
//==============================================================================
PostOperation {
{ Name WaveEqn_PostOp; NameOfPostProcessing WaveEqn_PostPro;
Operation {
Print[ p, OnElementsOf VolAll, File "res/p.pos"];
}
}
{ Name MembraneSys_PostOp; NameOfPostProcessing MembraneSys_PostPro;
Operation {
Print[ p, OnElementsOf VolAll, File "res/p.pos"];
Print[ v, OnElementsOf VolAll, File "res/u.pos"];
}
}
{ Name EulerSys_PostOp; NameOfPostProcessing EulerSys_PostPro;
Operation {
Print[ p, OnElementsOf VolAll, File "res/p.pos"];
Print[ u, OnElementsOf VolAll, File "res/u.pos"];
}
}
}
AcademicWaves/formulations_vectorWaves.pro
+ 22
− 52
View file @ dcfd65cc
......@@ -5,51 +5,21 @@
//========================================================
Function{
DefineFunction[ fNeumann, fDirichlet ] ;
DefineFunction[ fNeumann, fDirichlet ];
}
Group{
DefineGroup[ Omega, GammaD, GammaN, GammaR ] ;
SurAll = Region[{GammaD, GammaN, GammaR }] ;
VolAll = Region[{Omega}] ;
TotAll = Region[{VolAll,SurAll}] ;
}
Jacobian {
{ Name JVol ;
Case {
{ Region All ; Jacobian Vol ; }
}
}
{ Name JSur ;
Case {
{ Region All ; Jacobian Sur ; }
}
}
}
Integration {
{ Name I1 ;
Case {
{ Type Gauss ;
Case {
{ GeoElement Line ; NumberOfPoints 4 ; }
{ GeoElement Triangle ; NumberOfPoints 4 ; }
{ GeoElement Quadrangle ; NumberOfPoints 4 ; }
{ GeoElement Tetrahedron ; NumberOfPoints 4 ; }
{ GeoElement Hexahedron ; NumberOfPoints 6 ; }
{ GeoElement Prism ; NumberOfPoints 9 ; }
}
}
}
}
DefineGroup[ Omega, GammaD, GammaN, GammaR ];
SurAll = Region[{GammaD, GammaN, GammaR }];
VolAll = Region[{Omega}];
TotAll = Region[{VolAll,SurAll}];
}
FunctionSpace {
{ Name eSpace ; Type Form1 ;
{ Name eSpace; Type Form1;
BasisFunction {
{ Name se ; NameOfCoef ee ; Function BF_Edge ;
Support TotAll ; Entity EdgesOf[All] ; }
{ Name se; NameOfCoef ee; Function BF_Edge;
Support TotAll; Entity EdgesOf[All]; }
}
}
{ Name eSpace_DirichletBC; Type Form1;
......@@ -61,16 +31,16 @@ FunctionSpace {
}
Formulation {
{ Name Form ; Type FemEquation ;
{ Name Form; Type FemEquation;
Quantity {
{ Name e ; Type Local ; NameOfSpace eSpace ; }
{ Name lambda ; Type Local ; NameOfSpace eSpace_DirichletBC ; }
{ Name e; Type Local; NameOfSpace eSpace; }
{ Name lambda; Type Local; NameOfSpace eSpace_DirichletBC; }
}
Equation {
Galerkin { DtDtDof [ epsilon[] * Dof{e} , {e} ] ;
In VolAll ; Integration I1 ; Jacobian JVol ; }
Galerkin { [ nu[] * Dof{d e} , {d e} ] ;
In VolAll ; Integration I1 ; Jacobian JVol ; }
Galerkin { DtDtDof [ epsilon[] * Dof{e} , {e} ];
In VolAll; Integration I1; Jacobian JVol; }
Galerkin { [ nu[] * Dof{d e} , {d e} ];
In VolAll; Integration I1; Jacobian JVol; }
// BC: Dirichlet
Galerkin { [ Dof{lambda} , {e} ];
......@@ -81,21 +51,21 @@ Formulation {
In GammaD; Integration I1; Jacobian JSur; }
// BC: Radiation
Galerkin { [ I[]*k * nu[] * ( extNormal[] /\ Dof{e} ) /\ extNormal[] , {e} ] ;
In GammaR ; Integration I1 ; Jacobian JSur ; }
Galerkin { [ I[]*k * nu[] * ( extNormal[] /\ Dof{e} ) /\ extNormal[] , {e} ];
In GammaR; Integration I1; Jacobian JSur; }
}
}
}
Resolution {
{ Name Reso ;
{ Name Reso;
System {
{ Name A ; NameOfFormulation Form ; Type Complex ; Frequency Freq ; }
{ Name A; NameOfFormulation Form; Type Complex; Frequency FREQ; }
}
Operation {
Generate[A] ;
Solve[A] ;
SaveSolution[A] ;
Generate[A];
Solve[A];
SaveSolution[A];
}
}
}
AcademicWaves/gaussInCuboid.dat
+ 8
− 9
View file @ dcfd65cc
......@@ -3,6 +3,9 @@
// File: parameters
//========================================================
FLAG_DIM = 3;
FLAG_RES = RES_TIME;
//---------------------
// Physical parameters
//---------------------
......@@ -15,16 +18,12 @@ tf = 1;
// Numerical parameters
//---------------------
lc = 0.1;
Flag_SecondOrder = 0;
dt = lc/4;
// For the theta-scheme (if wave system)
theta = 0.5;
LC = 0.1;
dt = LC/4;
// For the Newmark scheme (if wave equation)
beta = 0.25;
gamma = 0.5;
theta = 0.5; // for theta-scheme (if wave system)
beta = 0.25; // for Newmark scheme (if wave equation)
gamma = 0.5; // for Newmark scheme (if wave equation)
//---------------------
// Geometrical tags
......
AcademicWaves/gaussInCuboid.geo
+ 3
− 1
View file @ dcfd65cc
......@@ -5,7 +5,9 @@
Include "gaussInCuboid.dat";
Mesh.CharacteristicLengthMax = lc;
Mesh.CharacteristicLengthMax = LC;
Mesh.CharacteristicLengthFactor = 1;
Mesh.Optimize = 1;
// Points
......
AcademicWaves/gaussInCuboid.pro
+ 19
− 29
View file @ dcfd65cc
......@@ -3,49 +3,39 @@
// File: GetDP simulation
//========================================================
Include "gaussInCuboid.dat" ;
Include "gaussInCuboid.dat";
Group {
GammaD = Region[{}] ;
GammaS = Region[{BOUNDARY}] ;
Omega = Region[{DOMAIN}] ;
GammaD = Region[{}];
GammaN = Region[{}];
GammaR = Region[{BOUNDARY}];
Omega = Region[{DOMAIN}];
}
Function {
InitValue[] = Exp[- ((X[])^2 + (Y[])^2 + (Z[])^2) / (r*r)] ;
InitValue[] = Exp[-(X[]^2+Y[]^2+Z[]^2)/(r*r)];
}
Constraint {
{ Name pConstraint ;
{ Name pConstraint;
Case {
{ Region Omega ; Type Init ; Value InitValue[] ; }
{ Region Omega; Type Init; Value InitValue[]; }
}
}
{ Name uConstraint ;
{ Name uConstraint;
Case {
{ Region GammaD ; Type Assign ; Value 0. ; }
{ Region Omega ; Type Init ; Value 0. ; }
{ Region GammaD; Type Assign; Value 0.; }
{ Region Omega; Type Init; Value 0.; }
}
}
}
Include "form_scalarWaveSyst_time3D.pro" ;
//Include "form_scalarWaveEqn_time3D.pro" ;
Include "formulations_scalarWaves.pro";
PostProcessing {
{ Name PostPro ; NameOfFormulation Form ; NameOfSystem A ;
Quantity {
{ Name p ; Value{ Local{ [ {p} ] ; In VolAll ; } } }
{ Name u ; Value{ Local{ [ {u} ] ; In VolAll ; } } }
}
}
}
PostOperation {
{ Name PostOp ; NameOfPostProcessing PostPro ;
Operation {
Print[ p, OnElementsOf VolAll, File "res/p.pos"] ;
Print[ u, OnElementsOf VolAll, File "res/u.pos"] ;
}
}
}
DefineConstant[
C_ = {"-solve -pos -v2", Name "GetDP/9ComputeCommand", Visible 0},
R_ = {"WaveEqn_Reso", Name "GetDP/1ResolutionChoices", Visible 1,
Choices {"WaveEqn_Reso", "MembraneSys_Reso", "EulerSys_Reso"}},
P_ = {"WaveEqn_PostOp", Name "GetDP/2PostOperationChoices", Visible 1,
Choices {"WaveEqn_PostOp", "MembraneSys_PostOp", "EulerSys_PostOp"}}
];
AcademicWaves/gaussInSquare.dat
+ 9
− 10
View file @ dcfd65cc
......@@ -3,28 +3,27 @@
// File: parameters
//========================================================
FLAG_DIM = 2;
FLAG_RES = RES_TIME;
//---------------------
// Physical parameters
//---------------------
ldom = 1;
r = ldom/5 ;
tf = 1;
tf = 1.;
//---------------------
// Numerical parameters
//---------------------
lc = 0.025;
Flag_SecondOrder = 0;
dt = lc/4;
// For the theta-scheme (if wave system)
theta = 0.5;
LC = 0.025;
dt = LC/4;
// For the Newmark scheme (if wave equation)
beta = 0.25;
gamma = 0.5;
theta = 0.5; // for theta-scheme (if wave system)
beta = 0.25; // for Newmark scheme (if wave equation)
gamma = 0.5; // for Newmark scheme (if wave equation)
//---------------------
// Geometrical tags
......
AcademicWaves/gaussInSquare.geo
+ 8
− 4
View file @ dcfd65cc
......@@ -5,10 +5,14 @@
Include "gaussInSquare.dat";
Point(1) = {-ldom/2, -ldom/2, 0, lc};
Point(2) = { ldom/2, -ldom/2, 0, lc};
Point(3) = { ldom/2, ldom/2, 0, lc};
Point(4) = {-ldom/2, ldom/2, 0, lc};
Mesh.CharacteristicLengthMax = LC;
Mesh.CharacteristicLengthFactor = 1;
Mesh.Optimize = 1;
Point(1) = {-ldom/2, -ldom/2, 0};
Point(2) = { ldom/2, -ldom/2, 0};
Point(3) = { ldom/2, ldom/2, 0};
Point(4) = {-ldom/2, ldom/2, 0};
Line(1) = {1, 2};
Line(2) = {2, 3};
......
AcademicWaves/gaussInSquare.pro
+ 19
− 32
View file @ dcfd65cc
......@@ -3,52 +3,39 @@
// File: GetDP simulation
//========================================================
Include "gaussInSquare.dat" ;
Include "gaussInSquare.dat";
Group {
GammaD = Region[{}] ;
GammaN = Region[{BOUNDARY}] ;
GammaS = Region[{}] ;
Omega = Region[{DOMAIN}] ;
GammaD = Region[{}];
GammaN = Region[{}];
GammaR = Region[{BOUNDARY}];
Omega = Region[{DOMAIN}];
}
Function {
InitValue[] = Exp[- ((X[])^2 + (Y[])^2) / (r*r)] ;
fNeumann[] = 0.;
InitValue[] = Exp[-(X[]^2+Y[]^2)/(r*r)];
}
Constraint {
{ Name pConstraint ;
{ Name pConstraint;
Case {
{ Region Omega ; Type Init ; Value InitValue[] ; }
{ Region Omega; Type Init; Value InitValue[]; }
}
}
{ Name uConstraint ;
{ Name uConstraint;
Case {
{ Region GammaD ; Type Assign ; Value 0. ; }
{ Region Omega ; Type Init ; Value 0. ; }
{ Region GammaD; Type Assign; Value 0.; }
{ Region Omega; Type Init; Value 0.; }
}
}
}
//Include "form_scalarWaveSyst_time2D.pro" ;
//Include "form_scalarWaveSyst_time2D_membrane.pro" ;
Include "form_scalarWaveEqn_time.pro" ;
Include "formulations_scalarWaves.pro";
PostProcessing {
{ Name PostPro ; NameOfFormulation Form ; NameOfSystem A ;
Quantity {
{ Name p ; Value{ Local{ [ {p} ] ; In VolAll ; } } }
{ Name u ; Value{ Local{ [ {u} ] ; In VolAll ; } } }
}
}
}
PostOperation {
{ Name PostOp ; NameOfPostProcessing PostPro ;
Operation {
Print[ p, OnElementsOf VolAll, File "res/p.pos"] ;
Print[ u, OnElementsOf VolAll, File "res/u.pos"] ;
}
}
}
DefineConstant[
C_ = {"-solve -pos -v2", Name "GetDP/9ComputeCommand", Visible 0},
R_ = {"WaveEqn_Reso", Name "GetDP/1ResolutionChoices", Visible 1,
Choices {"WaveEqn_Reso", "MembraneSys_Reso", "EulerSys_Reso"}},
P_ = {"WaveEqn_PostOp", Name "GetDP/2PostOperationChoices", Visible 1,
Choices {"WaveEqn_PostOp", "MembraneSys_PostOp", "EulerSys_PostOp"}}
];
AcademicWaves/loudspeaker.dat
+ 7
− 26
View file @ dcfd65cc
......@@ -3,23 +3,13 @@
// File: parameters
//========================================================
FLAG_DIM = 2;
FLAG_RES = RES_FREQ;
//---------------------
// Model choices
// Geometrical parameters
//---------------------
DefineConstant[
Flag_ModelEquations = {1, Choices{1="Wave equation (p)", 2="Wave system (p and u)"},
Name "Input/01Model equations", Highlight "Black"}
];
//---------------------
// Physical parameters
//---------------------
// Geometry
feed = 0.1;
app = 0.5;
a = 0.1;
......@@ -28,24 +18,15 @@ thick = 0.03;
x0 = thick-(a+b+thick)/2;
Ro = 1;
//---------------------
// Incident signal
//---------------------
k = 20;
lambda = 2*Pi/k;
nlambda = 10;
FREQ = k/(2*Pi);
LC = 0.01;
theta = 0;
//---------------------
// Numerical parameters
//---------------------
lc0 = 0.02;
lc1 = 0.02;
ORDER = 1;
//---------------------
// Geometrical tags
//---------------------
......
AcademicWaves/loudspeaker.geo
+ 21
− 19
View file @ dcfd65cc
......@@ -5,24 +5,27 @@
Include "loudspeaker.dat";
Mesh.CharacteristicLengthMax = LC;
Mesh.CharacteristicLengthFactor = 1;
Mesh.Optimize = 1;
//--------------------------------------------------------
// The loudspeaker
//--------------------------------------------------------
cen0=newp ; Point(newp) = {0,0,0,lc0} ;
pObs[]+=newp ; Point(newp) = {x0, -feed/2, 0, lc0} ;
pObs[]+=newp ; Point(newp) = {x0, feed/2, 0, lc0} ;
pObs[]+=newp ; Point(newp) = {a+x0, feed/2, 0, lc0} ;
pObs[]+=newp ; Point(newp) = {a+b+x0, app/2, 0, lc0} ;
pObs[]+=newp ; Point(newp) = {a+b+x0, app/2+thick, 0, lc0} ;
pObs[]+=newp ; Point(newp) = {a+x0, feed/2+thick, 0, lc0} ;
pObs[]+=newp ; Point(newp) = {-thick+x0, feed/2+thick, 0, lc0} ;
pObs[]+=newp ; Point(newp) = {-thick+x0, -feed/2-thick, 0, lc0} ;
pObs[]+=newp ; Point(newp) = {a+x0, -feed/2-thick, 0, lc0} ;
pObs[]+=newp ; Point(newp) = {a+b+x0, -app/2-thick, 0, lc0} ;
pObs[]+=newp ; Point(newp) = {a+b+x0, -app/2, 0, lc0} ;
pObs[]+=newp ; Point(newp) = {a+x0, -feed/2, 0, lc0} ;
cen0=newp ; Point(newp) = {0,0,0} ;
pObs[]+=newp ; Point(newp) = {x0, -feed/2, 0} ;
pObs[]+=newp ; Point(newp) = {x0, feed/2, 0} ;
pObs[]+=newp ; Point(newp) = {a+x0, feed/2, 0} ;
pObs[]+=newp ; Point(newp) = {a+b+x0, app/2, 0} ;
pObs[]+=newp ; Point(newp) = {a+b+x0, app/2+thick, 0} ;
pObs[]+=newp ; Point(newp) = {a+x0, feed/2+thick, 0} ;
pObs[]+=newp ; Point(newp) = {-thick+x0, feed/2+thick, 0} ;
pObs[]+=newp ; Point(newp) = {-thick+x0, -feed/2-thick, 0} ;
pObs[]+=newp ; Point(newp) = {a+x0, -feed/2-thick, 0} ;
pObs[]+=newp ; Point(newp) = {a+b+x0, -app/2-thick, 0} ;
pObs[]+=newp ; Point(newp) = {a+b+x0, -app/2, 0} ;
pObs[]+=newp ; Point(newp) = {a+x0, -feed/2, 0} ;
For i In {0:#pObs[]-1}
lObs[]+=newl ; Line(newl) = {pObs[i],pObs[(i+1)%#pObs[]]} ;
EndFor
......@@ -32,16 +35,15 @@ llObs=newll ; Line Loop(newll) = lObs[] ;
Physical Line(GAMMA_SOURCE) = {-lObs[{0}]} ;
Physical Line(GAMMA_WALL) = {-lObs[{1:11}]} ;
//--------------------------------------------------------
// Main domain with its boundary
//--------------------------------------------------------
cen0 = newp ; Point(newp) = {0,0,0,lc0} ;
pi[]+=newp ; Point(newp) = {Ro,0,0,lc1} ;
pi[]+=newp ; Point(newp) = {0,Ro,0,lc1} ;
pi[]+=newp ; Point(newp) = {-Ro,0,0,lc1} ;
pi[]+=newp ; Point(newp) = {0,-Ro,0,lc1} ;
cen0 = newp ; Point(newp) = {0,0,0} ;
pi[]+=newp ; Point(newp) = {Ro,0,0} ;
pi[]+=newp ; Point(newp) = {0,Ro,0} ;
pi[]+=newp ; Point(newp) = {-Ro,0,0} ;
pi[]+=newp ; Point(newp) = {0,-Ro,0} ;
li[]+=newl ; Circle(newl) = {pi[0], cen0, pi[1]} ;
li[]+=newl ; Circle(newl) = {pi[1], cen0, pi[2]} ;
li[]+=newl ; Circle(newl) = {pi[2], cen0, pi[3]} ;
......
AcademicWaves/loudspeaker.pro
+ 16
− 69
View file @ dcfd65cc
......@@ -3,92 +3,39 @@
// File: GetDP simulation
//========================================================
Include "loudspeaker.dat" ;
Include "loudspeaker.dat";
Group {
Omega = Region[{DOMAIN}] ;
If (Flag_ModelEquations==1)
GammaD = Region[{GAMMA_SOURCE}] ;
GammaN = Region[{GAMMA_WALL}] ;
EndIf
If (Flag_ModelEquations==2)
GammaD_p = Region[{GAMMA_SOURCE}] ;
GammaD_u = Region[{GAMMA_WALL}] ;
GammaD = Region[{GammaD_p,GammaD_u}] ;
EndIf
GammaR = Region[{GAMMA_DOM}] ;
Omega = Region[{DOMAIN}];
GammaD = Region[{GAMMA_SOURCE}];
GammaN = Region[{GAMMA_WALL}];
GammaR = Region[{GAMMA_DOM}];
}
Function {
I[] = Complex[0, 1] ;
I[] = Complex[0, 1];
// Incident signal
kVec[] = Vector[ k*Cos[theta], k*Sin[theta], 0] ;
pInc[] = Complex[ Cos[kVec[]*XYZ[]], Sin[kVec[]*XYZ[]] ] ;
kVec[] = Vector[ k*Cos[theta], k*Sin[theta], 0];
pInc[] = Complex[ Cos[kVec[]*XYZ[]], Sin[kVec[]*XYZ[]] ];
// BC: Neumann
If (Flag_ModelEquations==2)
fNeumann[] = 0. ;
EndIf
fNeumann[] = 0.;
}
Constraint {
{ Name pConstraint ;
Case {
{ Region Region[{GAMMA_SOURCE}] ; Value pInc[] ; }
}
{ Name pConstraint;
Case {{ Region GammaD; Value pInc[]; }}
}
If(Flag_SecondOrder)
{ Name pConstraint2 ;
Case {
{ Region Region[{GAMMA_SOURCE}] ; Value 0. ; }
}
{ Name uConstraint;
Case {{ Region GammaN; Value 0.; }}
}
EndIf
{ Name uConstraint ;
Case {
{ Region Region[{GAMMA_WALL}] ; Value 0. ; }
}
}
If(Flag_SecondOrder)
{ Name uConstraint2;
Case {
{ Region Region[{GAMMA_WALL}] ; Value 0. ; }
}
}
EndIf
}
Include "formulations_scalarWaves.pro" ;
PostProcessing {
If (Flag_ModelEquations==1)
{ Name PostPro ; NameOfFormulation WaveEquation_Form ;
Quantity {
{ Name p ; Value{ Local{ [ {p} ] ; In VolAll ; Jacobian JVol ; } } }
}
}
EndIf
If (Flag_ModelEquations==2)
{ Name PostPro ; NameOfFormulation EulerEquations_Form ;
Quantity {
{ Name p ; Value{ Local{ [ {p} ] ; In VolAll ; Jacobian JVol ; } } }
}
}
EndIf
}
PostOperation {
{ Name Get_Fields ; NameOfPostProcessing PostPro ;
Operation {
Print[ p, OnElementsOf VolAll, File "res/p.pos"] ;
}
}
}
Include "formulations_scalarWaves.pro";
DefineConstant[
C_ = {"-solve -pos -v2", Name "GetDP/9ComputeCommand", Visible 0},
P_ = {"Get_Fields", Name "GetDP/2PostOperationChoices", Visible 0, ReadOnly 1}
R_ = {"WaveEqn_Reso", Name "GetDP/1ResolutionChoices", Visible 1, Choices {"WaveEqn_Reso", "MembraneSys_Reso", "EulerSys_Reso"} },
P_ = {"WaveEqn_PostOp", Name "GetDP/2PostOperationChoices", Visible 1, Choices {"WaveEqn_PostOp", "MembraneSys_PostOp", "EulerSys_PostOp"}}
];
AcademicWaves/scattKite.dat
+ 4
− 16
View file @ dcfd65cc
......@@ -3,14 +3,14 @@
// File: parameters
//========================================================
FLAG_DIM = 2;
FLAG_RES = RES_FREQ;
//---------------------
// Model choices
//---------------------
DefineConstant[
Flag_ModelEquations = {1, Choices{1="Wave equation (p)", 2="Wave system (p and u)"},
Name "Input/01Model equations", Highlight "Black"}
Flag_DomTruncShape = {1, Choices{1="Circular domain", 2="Squared domain"},
Name "Input/1Domain truncation/01Shape of the domain", Highlight "Black"},
Flag_DomTruncMethod_Shape1 = {1, Choices{1="Sommerfeld BC", 2="Bayliss-Turkel BC"},
......@@ -25,13 +25,10 @@ If (Flag_DomTruncShape==2)
Flag_DomTruncMethod = Flag_DomTruncMethod_Shape2;
EndIf
//---------------------
// Physical parameters
//---------------------
DIM = 2;
// Incident signal
DefineConstant[
......@@ -42,7 +39,7 @@ DefineConstant[
theta = { 0., Min 0., Max 2*Pi, Step 0.01,
Name "Input/2Incident signal/12Theta", Highlight "Yellow"}
];
FREQ = k/(2*Pi);
// Geometry -- Fig. 7, page 26 (Ecevit, Reitich, "Multiple scattering iterations for high frequency")
......@@ -67,15 +64,7 @@ DefineConstant[
// Numerical parameters
//---------------------
DefineConstant[
lc = { 0.05, Min 0.001, Max 1, Step 0.001,
Name "Input/2Numerical parameters/10Size of a mesh cell", Highlight "Yellow"}
];
lc0 = lc;
lc1 = lc;
ORDER = 1;
LC = lambda/10;
//---------------------
// Geometrical tags
......@@ -88,4 +77,3 @@ GAMMA_OBS2 = 202;
DOMAIN = 301;
LAYER = 302;
AcademicWaves/scattKite.geo
+ 72
− 70
View file @ dcfd65cc
......@@ -3,33 +3,35 @@
// File: GMSH geometry
//========================================================
Include "scattKite.dat" ;
Include "scattKite.dat";
Mesh.CharacteristicLengthMax = LC;
Mesh.CharacteristicLengthFactor = 1;
Mesh.Optimize = 1;
//--------------------------------------------------------
// Obstacle 1
//--------------------------------------------------------
cen0=newp ; Point(newp) = {0,0,0,lc0} ;
pObs[]+=newp ; Point(newp) = {R0,0,0,lc0} ;
pObs[]+=newp ; Point(newp) = {0,R1,0,lc0} ;
pObs[]+=newp ; Point(newp) = {-R0,0,0,lc0} ;
pObs[]+=newp ; Point(newp) = {0,-R1,0,lc0} ;
lObs[]+=newl ; Ellipsis(newl) = {pObs[0], cen0, pObs[0], pObs[1]} ;
lObs[]+=newl ; Ellipsis(newl) = {pObs[1], cen0, pObs[1], pObs[2]} ;
lObs[]+=newl ; Ellipsis(newl) = {pObs[2], cen0, pObs[2], pObs[3]} ;
lObs[]+=newl ; Ellipsis(newl) = {pObs[3], cen0, pObs[3], pObs[0]} ;
cen0=newp; Point(newp) = {0,0,0};
pObs[]+=newp; Point(newp) = {R0,0,0};
pObs[]+=newp; Point(newp) = {0,R1,0};
pObs[]+=newp; Point(newp) = {-R0,0,0};
pObs[]+=newp; Point(newp) = {0,-R1,0};
lObs[]+=newl; Ellipsis(newl) = {pObs[0], cen0, pObs[0], pObs[1]};
lObs[]+=newl; Ellipsis(newl) = {pObs[1], cen0, pObs[1], pObs[2]};
lObs[]+=newl; Ellipsis(newl) = {pObs[2], cen0, pObs[2], pObs[3]};
lObs[]+=newl; Ellipsis(newl) = {pObs[3], cen0, pObs[3], pObs[0]};
llellipsis=newll ; Line Loop(newll) = {lObs[]} ;
surfellipsis=news ; Plane Surface(news) = {llellipsis} ;
llellipsis=newll; Line Loop(newll) = {lObs[]};
surfellipsis=news; Plane Surface(news) = {llellipsis};
ll[] += llellipsis ;
ll[] += llellipsis;
Rotate {{0, 0, 1}, {0, 0, 0}, phi} { Surface{surfellipsis}; }
Translate {cenx, ceny, 0} { Surface{surfellipsis}; }
Physical Line(GAMMA_OBS1) = {-lObs[]} ;
Physical Line(GAMMA_OBS1) = {-lObs[]};
//--------------------------------------------------------
// Obstacle 2 : Kite shaped domain
......@@ -37,19 +39,19 @@ Physical Line(GAMMA_OBS1) = {-lObs[]} ;
// .. 0 <= t <= 2*Pi
//--------------------------------------------------------
t[] = {0:2*Pi:(2*Pi/50)} ;
t[] = {0:2*Pi:(2*Pi/50)};
For kk In {0:#t[]-1}
pk[]+=newp ; Point(newp) = { factor*(Cos(t[kk])+0.65*Cos(2*t[kk])-0.65), factor*(1.5*Sin(t[kk])), 0., lc0 } ;
pk[]+=newp; Point(newp) = { factor*(Cos(t[kk])+0.65*Cos(2*t[kk])-0.65), factor*(1.5*Sin(t[kk])), 0.};
EndFor
lkite=newl ; Spline(newl) = {pk[],pk[0]} ;
llkite=newll ; Line Loop(newll) = {lkite} ;
skite=news ; Plane Surface(news) = {newll-1} ;
lkite=newl; Spline(newl) = {pk[],pk[0]};
ll[] += llkite ;
llkite=newll; Line Loop(newll) = {lkite};
skite=news; Plane Surface(news) = {newll-1};
Physical Line(GAMMA_OBS2) = {-lkite} ;
ll[] += llkite;
Physical Line(GAMMA_OBS2) = {-lkite};
//--------------------------------------------------------
// Main domain and Layers with the boundaries (CIRCLE)
......@@ -57,38 +59,38 @@ Physical Line(GAMMA_OBS2) = {-lkite} ;
If (Flag_DomTruncShape==1)
cen0 = newp ; Point(newp) = {0,0,0,lc0} ;
pi[]+=newp ; Point(newp) = {Rdom,0,0,lc1} ;
pi[]+=newp ; Point(newp) = {0,Rdom,0,lc1} ;
pi[]+=newp ; Point(newp) = {-Rdom,0,0,lc1} ;
pi[]+=newp ; Point(newp) = {0,-Rdom,0,lc1} ;
li[]+=newl ; Circle(newl) = {pi[0], cen0, pi[1]} ;
li[]+=newl ; Circle(newl) = {pi[1], cen0, pi[2]} ;
li[]+=newl ; Circle(newl) = {pi[2], cen0, pi[3]} ;
li[]+=newl ; Circle(newl) = {pi[3], cen0, pi[0]} ;
cen0 = newp; Point(newp) = {0,0,0};
pi[]+=newp; Point(newp) = {Rdom,0,0};
pi[]+=newp; Point(newp) = {0,Rdom,0};
pi[]+=newp; Point(newp) = {-Rdom,0,0};
pi[]+=newp; Point(newp) = {0,-Rdom,0};
li[]+=newl; Circle(newl) = {pi[0], cen0, pi[1]};
li[]+=newl; Circle(newl) = {pi[1], cen0, pi[2]};
li[]+=newl; Circle(newl) = {pi[2], cen0, pi[3]};
li[]+=newl; Circle(newl) = {pi[3], cen0, pi[0]};
llOmega=newll ; Line Loop(llOmega) = li[];
surfOmega=news ; Plane Surface(news) = {llOmega, ll[]};
llOmega=newll; Line Loop(llOmega) = li[];
surfOmega=news; Plane Surface(news) = {llOmega, ll[]};
Physical Line(GAMMA_DOM) = li[] ;
Physical Line(GAMMA_DOM) = li[];
Physical Surface(DOMAIN) = {surfOmega};
If (Flag_DomTruncMethod==9)
cen0 = newp ; Point(newp) = {0,0,0,lc0} ;
pe[]+=newp ; Point(newp) = {Rdom+Lpml,0,0,lc1} ;
pe[]+=newp ; Point(newp) = {0,Rdom+Lpml,0,lc1} ;
pe[]+=newp ; Point(newp) = {-Rdom-Lpml,0,0,lc1} ;
pe[]+=newp ; Point(newp) = {0,-Rdom-Lpml,0,lc1} ;
le[]+=newl ; Circle(newl) = {pe[0], cen0, pe[1]} ;
le[]+=newl ; Circle(newl) = {pe[1], cen0, pe[2]} ;
le[]+=newl ; Circle(newl) = {pe[2], cen0, pe[3]} ;
le[]+=newl ; Circle(newl) = {pe[3], cen0, pe[0]} ;
cen0 = newp; Point(newp) = {0,0,0};
pe[]+=newp; Point(newp) = {Rdom+Lpml,0,0};
pe[]+=newp; Point(newp) = {0,Rdom+Lpml,0};
pe[]+=newp; Point(newp) = {-Rdom-Lpml,0,0};
pe[]+=newp; Point(newp) = {0,-Rdom-Lpml,0};
le[]+=newl; Circle(newl) = {pe[0], cen0, pe[1]};
le[]+=newl; Circle(newl) = {pe[1], cen0, pe[2]};
le[]+=newl; Circle(newl) = {pe[2], cen0, pe[3]};
le[]+=newl; Circle(newl) = {pe[3], cen0, pe[0]};
llLayer=newll ; Line Loop(llLayer) = le[];
surfLayer=news ; Plane Surface(news) = {llOmega, llLayer};
llLayer=newll; Line Loop(llLayer) = le[];
surfLayer=news; Plane Surface(news) = {llOmega, llLayer};
Physical Line(GAMMA_LAY) = le[] ;
Physical Line(GAMMA_LAY) = le[];
Physical Surface(LAYER) = {surfLayer};
EndIf
......@@ -102,37 +104,37 @@ EndIf
If (Flag_DomTruncShape==2)
cen0=newp ; Point(newp) = {0,0,0,lc0} ;
pi[]+=newp ; Point(newp) = { Ldom/2, Ldom/2,0,lc1} ;
pi[]+=newp ; Point(newp) = {-Ldom/2, Ldom/2,0,lc1} ;
pi[]+=newp ; Point(newp) = {-Ldom/2,-Ldom/2,0,lc1} ;
pi[]+=newp ; Point(newp) = { Ldom/2,-Ldom/2,0,lc1} ;
li[]+=newl ; Line(newl) = {pi[0], pi[1]} ;
li[]+=newl ; Line(newl) = {pi[1], pi[2]} ;
li[]+=newl ; Line(newl) = {pi[2], pi[3]} ;
li[]+=newl ; Line(newl) = {pi[3], pi[0]} ;
cen0=newp; Point(newp) = {0,0,0};
pi[]+=newp; Point(newp) = { Ldom/2, Ldom/2,0};
pi[]+=newp; Point(newp) = {-Ldom/2, Ldom/2,0};
pi[]+=newp; Point(newp) = {-Ldom/2,-Ldom/2,0};
pi[]+=newp; Point(newp) = { Ldom/2,-Ldom/2,0};
li[]+=newl; Line(newl) = {pi[0], pi[1]};
li[]+=newl; Line(newl) = {pi[1], pi[2]};
li[]+=newl; Line(newl) = {pi[2], pi[3]};
li[]+=newl; Line(newl) = {pi[3], pi[0]};
llOmega=newll ; Line Loop(llOmega) = li[];
surfOmega=news ; Plane Surface(news) = {llOmega, ll[]};
llOmega=newll; Line Loop(llOmega) = li[];
surfOmega=news; Plane Surface(news) = {llOmega, ll[]};
Physical Line(GAMMA_DOM) = li[] ;
Physical Line(GAMMA_DOM) = li[];
Physical Surface(DOMAIN) = {surfOmega};
If (Flag_DomTruncMethod==9)
pe[]+=newp ; Point(newp) = { Ldom/2+Lpml, Ldom/2+Lpml,0,lc1} ;
pe[]+=newp ; Point(newp) = {-Ldom/2-Lpml, Ldom/2+Lpml,0,lc1} ;
pe[]+=newp ; Point(newp) = {-Ldom/2-Lpml,-Ldom/2-Lpml,0,lc1} ;
pe[]+=newp ; Point(newp) = { Ldom/2+Lpml,-Ldom/2-Lpml,0,lc1} ;
le[]+=newl ; Line(newl) = {pe[0], pe[1]} ;
le[]+=newl ; Line(newl) = {pe[1], pe[2]} ;
le[]+=newl ; Line(newl) = {pe[2], pe[3]} ;
le[]+=newl ; Line(newl) = {pe[3], pe[0]} ;
pe[]+=newp; Point(newp) = { Ldom/2+Lpml, Ldom/2+Lpml,0};
pe[]+=newp; Point(newp) = {-Ldom/2-Lpml, Ldom/2+Lpml,0};
pe[]+=newp; Point(newp) = {-Ldom/2-Lpml,-Ldom/2-Lpml,0};
pe[]+=newp; Point(newp) = { Ldom/2+Lpml,-Ldom/2-Lpml,0};
le[]+=newl; Line(newl) = {pe[0], pe[1]};
le[]+=newl; Line(newl) = {pe[1], pe[2]};
le[]+=newl; Line(newl) = {pe[2], pe[3]};
le[]+=newl; Line(newl) = {pe[3], pe[0]};
llLayer=newll ; Line Loop(llLayer) = le[];
surfLayer=news ; Plane Surface(news) = {llOmega, llLayer};
llLayer=newll; Line Loop(llLayer) = le[];
surfLayer=news; Plane Surface(news) = {llOmega, llLayer};
Physical Line(GAMMA_LAY) = le[] ;
Physical Line(GAMMA_LAY) = le[];
Physical Surface(LAYER) = {surfLayer};
EndIf
......
AcademicWaves/scattKite.pro
+ 34
− 61
View file @ dcfd65cc
......@@ -3,113 +3,86 @@
// File: GetDP simulation
//========================================================
Include "scattKite.dat" ;
Include "scattKite.dat";
Group {
Omega = Region[{DOMAIN}] ;
GammaD = Region[{GAMMA_OBS1, GAMMA_OBS2}] ;
Omega = Region[{DOMAIN}];
GammaD = Region[{GAMMA_OBS1, GAMMA_OBS2}];
If(Flag_DomTruncMethod==1)
GammaR = Region[{GAMMA_DOM}] ;
GammaR = Region[{GAMMA_DOM}];
EndIf
If(Flag_DomTruncMethod==2)
GammaBT = Region[{GAMMA_DOM}] ;
GammaBT = Region[{GAMMA_DOM}];
EndIf
If(Flag_DomTruncMethod==9)
OmegaPml = Region[{LAYER}] ;
GammaD += Region[{GAMMA_LAY}] ;
OmegaPml = Region[{LAYER}];
GammaD += Region[{GAMMA_LAY}];
EndIf
}
Function {
I[] = Complex[0, 1] ;
Freq = 1./lambda;
kVec[] = Vector[ k*Cos[theta], k*Sin[theta], 0] ;
pInc[] = Complex[ Cos[kVec[]*XYZ[]], Sin[kVec[]*XYZ[]] ] ;
uInc[] = -pInc[] * kVec[]/k ;
I[] = Complex[0, 1];
// Incident signal
kVec[] = Vector[ k*Cos[theta], k*Sin[theta], 0];
pInc[] = Complex[ Cos[kVec[]*XYZ[]], Sin[kVec[]*XYZ[]] ];
// BC: Bayliss-Turkel
If(Flag_DomTruncMethod==2)
alphaBT[] = 1/(2*Ldom) - I[]/(8*k*Ldom^2*(1+I[]/(k*Ldom))) ;
betaBT[] = - 1/(2*I[]*k*(1+I[]/(k*Ldom))) ;
alphaBT[] = 1/(2*Ldom) - I[]/(8*k*Ldom^2*(1+I[]/(k*Ldom)));
betaBT[] = - 1/(2*I[]*k*(1+I[]/(k*Ldom)));
EndIf
// PML
If(Flag_DomTruncMethod==9)
xLoc[] = Fabs[X[]]-Ldom/2;
yLoc[] = Fabs[Y[]]-Ldom/2;
absFuncX[] = (xLoc[]>=0) ? 1/(Lpml-xLoc[]) : 0 ;
absFuncY[] = (yLoc[]>=0) ? 1/(Lpml-yLoc[]) : 0 ;
absFuncX[] = (xLoc[]>=0) ? Log[Lpml-xLoc[]] : 0;
absFuncY[] = (yLoc[]>=0) ? Log[Lpml-yLoc[]] : 0;
hx[] = Complex[1, absFuncX[]/k];
hy[] = Complex[1, absFuncY[]/k];
layerScal[OmegaPml] = 1/(hx[]*hy[]);
layerTens[OmegaPml] = TensorDiag[hy[]/hx[], hx[]/hy[], 1.];
pmlScal[OmegaPml] = 1/(hx[]*hy[]);
pmlTens[OmegaPml] = TensorDiag[hy[]/hx[], hx[]/hy[], 1.];
EndIf
}
Constraint {
{ Name pConstraint ;
{ Name pConstraint;
Case {
{ Region Region[{GAMMA_OBS1, GAMMA_OBS2}] ; Value -pInc[] ; }
{ Region Region[{GAMMA_OBS1, GAMMA_OBS2}]; Value -pInc[]; }
If(Flag_DomTruncMethod==9)
{ Region Region[{GAMMA_LAY}] ; Value 0. ; }
{ Region Region[{GAMMA_LAY}]; Value 0.; }
EndIf
}
}
{ Name uConstraint ;
Case {
{ }
}
}
{ Name uConstraint ; Case {} }
}
Include "formulations_scalarWaves.pro" ;
Include "formulations_scalarWaves.pro";
PostProcessing {
If (Flag_ModelEquations==1)
{ Name PostPro ; NameOfFormulation WaveEquation_Form ;
{ Name WaveEqn2_PostPro; NameOfFormulation WaveEqn_Form;
Quantity {
{ Name pSct ; Value{ Local{ [ {p} ] ; In VolAll ; Jacobian JVol ; } } }
{ Name pInc ; Value{ Local{ [ pInc[] ] ; In VolAll ; Jacobian JVol ; } } }
{ Name pTot ; Value{ Local{ [ pInc[]+{p} ] ; In VolAll ; Jacobian JVol ; } } }
{ Name pSct; Value{ Local{ [ {p} ]; In VolAll; Jacobian JVol; }}}
{ Name pInc; Value{ Local{ [ pInc[] ]; In VolAll; Jacobian JVol; }}}
{ Name pTot; Value{ Local{ [ pInc[]+{p} ]; In VolAll; Jacobian JVol; }}}
}
}
EndIf
If (Flag_ModelEquations==2)
{ Name PostPro ; NameOfFormulation EulerEquations_Form ;
Quantity {
{ Name pSct ; Value{ Local{ [ {p} ] ; In VolAll ; Jacobian JVol ; } } }
{ Name pInc ; Value{ Local{ [ pInc[] ] ; In VolAll ; Jacobian JVol ; } } }
{ Name pTot ; Value{ Local{ [ pInc[]+{p} ] ; In VolAll ; Jacobian JVol ; } } }
{ Name uSct ; Value{ Local{ [ {u} ] ; In VolAll ; Jacobian JVol ; } } }
{ Name uInc ; Value{ Local{ [ uInc[] ] ; In VolAll ; Jacobian JVol ; } } }
{ Name uTot ; Value{ Local{ [ uInc[]+{u} ] ; In VolAll ; Jacobian JVol ; } } }
}
}
EndIf
}
PostOperation {
{ Name Get_Fields ; NameOfPostProcessing PostPro ;
{ Name WaveEqn2_PostOp; NameOfPostProcessing WaveEqn2_PostPro;
Operation {
Print[ pSct, OnElementsOf VolAll, File "res/pSct.pos"] ;
//Print[ pInc, OnElementsOf VolAll, File "res/pInc.pos"] ;
//Print[ pTot, OnElementsOf VolAll, File "res/pTot.pos"] ;
If (Flag_ModelEquations==2)
//Print[ uSct, OnElementsOf VolAll, File "res/uSct.pos"] ;
//Print[ uInc, OnElementsOf VolAll, File "res/uInc.pos"] ;
//Print[ uTot, OnElementsOf VolAll, File "res/uTot.pos"] ;
EndIf
Print[ pSct, OnElementsOf VolAll, File "res/pSct.pos"];
Print[ pInc, OnElementsOf VolAll, File "res/pInc.pos"];
Print[ pTot, OnElementsOf VolAll, File "res/pTot.pos"];
}
}
}
DefineConstant[
C_ = {"-solve -pos -v2", Name "GetDP/9ComputeCommand", Visible 0},
P_ = {"Get_Fields", Name "GetDP/2PostOperationChoices", Visible 0, ReadOnly 1}
R_ = {"WaveEqn_Reso", Name "GetDP/1ResolutionChoices", Visible 0, Choices {"WaveEqn_Reso"} },
P_ = {"WaveEqn2_PostOp", Name "GetDP/2PostOperationChoices", Visible 0, Choices {"WaveEqn2_PostOp"}}
];
AcademicWaves/scattPlates.dat
+ 4
− 12
View file @ dcfd65cc
......@@ -3,14 +3,14 @@
// File: parameters
//========================================================
FLAG_DIM = 2;
FLAG_RES = RES_FREQ;
//---------------------
// Model choices
//---------------------
DefineConstant[
Flag_ModelEquations = {1, Choices{1="Wave equation (p)", 2="Wave system (p and u)"},
Name "Input/01Model equations", Highlight "Black"}
Flag_DomTruncShape = {1, Choices{1="Circular domain", 2="Squared domain"},
Name "Input/1Domain truncation/01Shape of the domain", Highlight "Black"},
Flag_DomTruncMethod_Shape1 = {1, Choices{1="Sommerfeld BC", 2="Bayliss-Turkel BC"},
......@@ -40,6 +40,7 @@ DefineConstant[
theta = { Pi/4, Min 0., Max 2*Pi, Step 0.01,
Name "Input/2Incident signal/12Theta", Highlight "Yellow"}
];
FREQ = k/(2*Pi);
// Geometry
......@@ -57,20 +58,11 @@ DefineConstant[
Name "Input/1Domain truncation/12Thickness of the layer", Highlight "Yellow", Visible (Flag_DomTruncMethod==9)}
];
//---------------------
// Numerical parameters
//---------------------
DefineConstant[
lc = { 0.01, Min 0.001, Max 1, Step 0.001,
Name "Input/2Numerical parameters/10Size of a mesh cell", Highlight "Yellow"}
];
lc0 = lc;
lc1 = lc;
Flag_SecondOrder = 0;
LC = lambda/10;
//---------------------
// Geometrical tags
......
AcademicWaves/scattPlates.geo
+ 30
− 27
View file @ dcfd65cc
......@@ -5,15 +5,18 @@
Include "scattPlates.dat";
Mesh.CharacteristicLengthMax = LC;
Mesh.CharacteristicLengthFactor = 1;
Mesh.Optimize = 1;
//--------------------------------------------------------
// Obstacle 1
//--------------------------------------------------------
pObs1[]+=newp ; Point(newp) = { dx/2,-dy/2+dist, 0, lc0} ;
pObs1[]+=newp ; Point(newp) = { dx/2, dy/2+dist, 0, lc0} ;
pObs1[]+=newp ; Point(newp) = {-dx/2, dy/2+dist, 0, lc0} ;
pObs1[]+=newp ; Point(newp) = {-dx/2,-dy/2+dist, 0, lc0} ;
pObs1[]+=newp ; Point(newp) = { dx/2,-dy/2+dist, 0} ;
pObs1[]+=newp ; Point(newp) = { dx/2, dy/2+dist, 0} ;
pObs1[]+=newp ; Point(newp) = {-dx/2, dy/2+dist, 0} ;
pObs1[]+=newp ; Point(newp) = {-dx/2,-dy/2+dist, 0} ;
lObs1[]+=newl ; Line(newl) = {pObs1[0],pObs1[1]} ;
lObs1[]+=newl ; Line(newl) = {pObs1[1],pObs1[2]} ;
......@@ -30,10 +33,10 @@ Physical Line(GAMMA_OBS1) = {-lObs1[]} ;
// Obstacle 2
//--------------------------------------------------------
pObs2[]+=newp ; Point(newp) = { dx/2,-dy/2-dist, 0, lc0} ;
pObs2[]+=newp ; Point(newp) = { dx/2, dy/2-dist, 0, lc0} ;
pObs2[]+=newp ; Point(newp) = {-dx/2, dy/2-dist, 0, lc0} ;
pObs2[]+=newp ; Point(newp) = {-dx/2,-dy/2-dist, 0, lc0} ;
pObs2[]+=newp ; Point(newp) = { dx/2,-dy/2-dist, 0} ;
pObs2[]+=newp ; Point(newp) = { dx/2, dy/2-dist, 0} ;
pObs2[]+=newp ; Point(newp) = {-dx/2, dy/2-dist, 0} ;
pObs2[]+=newp ; Point(newp) = {-dx/2,-dy/2-dist, 0} ;
lObs2[]+=newl ; Line(newl) = {pObs2[0],pObs2[1]} ;
lObs2[]+=newl ; Line(newl) = {pObs2[1],pObs2[2]} ;
......@@ -52,11 +55,11 @@ Physical Line(GAMMA_OBS2) = {-lObs2[]} ;
If (Flag_DomTruncShape==1)
cen0 = newp ; Point(newp) = {0,0,0,lc0} ;
pi[]+=newp ; Point(newp) = {Rdom,0,0,lc1} ;
pi[]+=newp ; Point(newp) = {0,Rdom,0,lc1} ;
pi[]+=newp ; Point(newp) = {-Rdom,0,0,lc1} ;
pi[]+=newp ; Point(newp) = {0,-Rdom,0,lc1} ;
cen0 = newp ; Point(newp) = {0,0,0} ;
pi[]+=newp ; Point(newp) = {Rdom,0,0} ;
pi[]+=newp ; Point(newp) = {0,Rdom,0} ;
pi[]+=newp ; Point(newp) = {-Rdom,0,0} ;
pi[]+=newp ; Point(newp) = {0,-Rdom,0} ;
li[]+=newl ; Circle(newl) = {pi[0], cen0, pi[1]} ;
li[]+=newl ; Circle(newl) = {pi[1], cen0, pi[2]} ;
li[]+=newl ; Circle(newl) = {pi[2], cen0, pi[3]} ;
......@@ -70,11 +73,11 @@ If (Flag_DomTruncShape==1)
If (Flag_DomTruncMethod==9)
cen0 = newp ; Point(newp) = {0,0,0,lc0} ;
pe[]+=newp ; Point(newp) = {Rdom+Lpml,0,0,lc1} ;
pe[]+=newp ; Point(newp) = {0,Rdom+Lpml,0,lc1} ;
pe[]+=newp ; Point(newp) = {-Rdom-Lpml,0,0,lc1} ;
pe[]+=newp ; Point(newp) = {0,-Rdom-Lpml,0,lc1} ;
cen0 = newp ; Point(newp) = {0,0,0} ;
pe[]+=newp ; Point(newp) = {Rdom+Lpml,0,0} ;
pe[]+=newp ; Point(newp) = {0,Rdom+Lpml,0} ;
pe[]+=newp ; Point(newp) = {-Rdom-Lpml,0,0} ;
pe[]+=newp ; Point(newp) = {0,-Rdom-Lpml,0} ;
le[]+=newl ; Circle(newl) = {pe[0], cen0, pe[1]} ;
le[]+=newl ; Circle(newl) = {pe[1], cen0, pe[2]} ;
le[]+=newl ; Circle(newl) = {pe[2], cen0, pe[3]} ;
......@@ -97,11 +100,11 @@ EndIf
If (Flag_DomTruncShape==2)
cen0=newp ; Point(newp) = {0,0,0,lc0} ;
pi[]+=newp ; Point(newp) = { Ldom/2, Ldom/2,0,lc1} ;
pi[]+=newp ; Point(newp) = {-Ldom/2, Ldom/2,0,lc1} ;
pi[]+=newp ; Point(newp) = {-Ldom/2,-Ldom/2,0,lc1} ;
pi[]+=newp ; Point(newp) = { Ldom/2,-Ldom/2,0,lc1} ;
cen0=newp ; Point(newp) = {0,0,0} ;
pi[]+=newp ; Point(newp) = { Ldom/2, Ldom/2,0} ;
pi[]+=newp ; Point(newp) = {-Ldom/2, Ldom/2,0} ;
pi[]+=newp ; Point(newp) = {-Ldom/2,-Ldom/2,0} ;
pi[]+=newp ; Point(newp) = { Ldom/2,-Ldom/2,0} ;
li[]+=newl ; Line(newl) = {pi[0], pi[1]} ;
li[]+=newl ; Line(newl) = {pi[1], pi[2]} ;
li[]+=newl ; Line(newl) = {pi[2], pi[3]} ;
......@@ -115,10 +118,10 @@ If (Flag_DomTruncShape==2)
If (Flag_DomTruncMethod==9)
pe[]+=newp ; Point(newp) = { Ldom/2+Lpml, Ldom/2+Lpml,0,lc1} ;
pe[]+=newp ; Point(newp) = {-Ldom/2-Lpml, Ldom/2+Lpml,0,lc1} ;
pe[]+=newp ; Point(newp) = {-Ldom/2-Lpml,-Ldom/2-Lpml,0,lc1} ;
pe[]+=newp ; Point(newp) = { Ldom/2+Lpml,-Ldom/2-Lpml,0,lc1} ;
pe[]+=newp ; Point(newp) = { Ldom/2+Lpml, Ldom/2+Lpml,0} ;
pe[]+=newp ; Point(newp) = {-Ldom/2-Lpml, Ldom/2+Lpml,0} ;
pe[]+=newp ; Point(newp) = {-Ldom/2-Lpml,-Ldom/2-Lpml,0} ;
pe[]+=newp ; Point(newp) = { Ldom/2+Lpml,-Ldom/2-Lpml,0} ;
le[]+=newl ; Line(newl) = {pe[0], pe[1]} ;
le[]+=newl ; Line(newl) = {pe[1], pe[2]} ;
le[]+=newl ; Line(newl) = {pe[2], pe[3]} ;
......
AcademicWaves/scattPlates.pro
+ 35
− 57
View file @ dcfd65cc
......@@ -3,108 +3,86 @@
// File: GetDP simulation
//========================================================
Include "scattPlates.dat" ;
Include "scattPlates.dat";
Group {
Omega = Region[{DOMAIN}] ;
GammaD = Region[{GAMMA_OBS1, GAMMA_OBS2}] ;
Omega = Region[{DOMAIN}];
GammaD = Region[{GAMMA_OBS1, GAMMA_OBS2}];
If(Flag_DomTruncMethod==1)
GammaR = Region[{GAMMA_DOM}] ;
GammaR = Region[{GAMMA_DOM}];
EndIf
If(Flag_DomTruncMethod==2)
GammaBT = Region[{GAMMA_DOM}] ;
GammaBT = Region[{GAMMA_DOM}];
EndIf
If(Flag_DomTruncMethod==9)
OmegaPml = Region[{LAYER}] ;
GammaD += Region[{GAMMA_LAY}] ;
OmegaPml = Region[{LAYER}];
GammaD += Region[{GAMMA_LAY}];
EndIf
}
Function {
I[] = Complex[0, 1] ;
I[] = Complex[0, 1];
// Incident signal
kVec[] = Vector[ k*Cos[theta], k*Sin[theta], 0] ;
pInc[] = Complex[ Cos[kVec[]*XYZ[]], Sin[kVec[]*XYZ[]] ] ;
kVec[] = Vector[ k*Cos[theta], k*Sin[theta], 0];
pInc[] = Complex[ Cos[kVec[]*XYZ[]], Sin[kVec[]*XYZ[]] ];
// BC: Bayliss-Turkel
If(Flag_DomTruncMethod==2)
alphaBT[] = 1/(2*Ldom) - I[]/(8*k*Ldom^2*(1+I[]/(k*Ldom))) ;
betaBT[] = - 1/(2*I[]*k*(1+I[]/(k*Ldom))) ;
alphaBT[] = 1/(2*Ldom) - I[]/(8*k*Ldom^2*(1+I[]/(k*Ldom)));
betaBT[] = - 1/(2*I[]*k*(1+I[]/(k*Ldom)));
EndIf
// PML
If(Flag_DomTruncMethod==9)
xLoc[] = Fabs[X[]]-Ldom/2;
yLoc[] = Fabs[Y[]]-Ldom/2;
absFuncX[] = (xLoc[]>=0) ? 1/(Lpml-xLoc[]) : 0 ;
absFuncY[] = (yLoc[]>=0) ? 1/(Lpml-yLoc[]) : 0 ;
absFuncX[] = (xLoc[]>=0) ? Log[Lpml-xLoc[]] : 0;
absFuncY[] = (yLoc[]>=0) ? Log[Lpml-yLoc[]] : 0;
hx[] = Complex[1, absFuncX[]/k];
hy[] = Complex[1, absFuncY[]/k];
layerScal[OmegaPml] = 1/(hx[]*hy[]);
layerTens[OmegaPml] = TensorDiag[hy[]/hx[], hx[]/hy[], 1.];
pmlScal[OmegaPml] = 1/(hx[]*hy[]);
pmlTens[OmegaPml] = TensorDiag[hy[]/hx[], hx[]/hy[], 1.];
EndIf
}
Constraint {
{ Name pConstraint ;
{ Name pConstraint;
Case {
{ Region Region[{GAMMA_OBS1, GAMMA_OBS2}] ; Value -pInc[] ; }
{ Region Region[{GAMMA_OBS1, GAMMA_OBS2}]; Value -pInc[]; }
If(Flag_DomTruncMethod==9)
{ Region Region[{GAMMA_LAY}] ; Value 0. ; }
{ Region Region[{GAMMA_LAY}]; Value 0.; }
EndIf
}
}
{ Name uConstraint ;
Case {
}
}
{ Name uConstraint ; Case {} }
}
If (Flag_ModelEquations==1)
Include "form_scalarWaveEqn_frequency.pro" ;
EndIf
If (Flag_ModelEquations==2)
Include "form_scalarWaveSyst2D_frequency.pro" ;
EndIf
Include "formulations_scalarWaves.pro";
PostProcessing {
{ Name PostPro ; NameOfFormulation Form ;
{ Name WaveEqn2_PostPro; NameOfFormulation WaveEqn_Form;
Quantity {
{ Name pSct ; Value{ Local{ [ {p} ] ; In VolAll ; Jacobian JVol ; } } }
{ Name pInc ; Value{ Local{ [ pInc[] ] ; In VolAll ; Jacobian JVol ; } } }
{ Name pTot ; Value{ Local{ [ pInc[]+{p} ] ; In VolAll ; Jacobian JVol ; } } }
If (Flag_ModelEquations==2)
{ Name uSct ; Value{ Local{ [ {u} ] ; In VolAll ; Jacobian JVol ; } } }
{ Name uInc ; Value{ Local{ [ uInc[] ] ; In VolAll ; Jacobian JVol ; } } }
{ Name uTot ; Value{ Local{ [ uInc[]+{u} ] ; In VolAll ; Jacobian JVol ; } } }
EndIf
{ Name pSct; Value{ Local{ [ {p} ]; In VolAll; Jacobian JVol; }}}
{ Name pInc; Value{ Local{ [ pInc[] ]; In VolAll; Jacobian JVol; }}}
{ Name pTot; Value{ Local{ [ pInc[]+{p} ]; In VolAll; Jacobian JVol; }}}
}
}
}
PostOperation {
{ Name Get_Fields ; NameOfPostProcessing PostPro ;
{ Name WaveEqn2_PostOp; NameOfPostProcessing WaveEqn2_PostPro;
Operation {
Print[ pSct, OnElementsOf VolAll, File "res/pSct.pos"] ;
Print[ pInc, OnElementsOf VolAll, File "res/pInc.pos"] ;
Print[ pTot, OnElementsOf VolAll, File "res/pTot.pos"] ;
If (Flag_ModelEquations==2)
Print[ uSct, OnElementsOf VolAll, File "res/uSct.pos"] ;
Print[ uInc, OnElementsOf VolAll, File "res/uInc.pos"] ;
Print[ uTot, OnElementsOf VolAll, File "res/uTot.pos"] ;
EndIf
Print[ pSct, OnElementsOf VolAll, File "res/pSct.pos"];
Print[ pInc, OnElementsOf VolAll, File "res/pInc.pos"];
Print[ pTot, OnElementsOf VolAll, File "res/pTot.pos"];
}
}
}
DefineConstant[
R_ = {"Reso", Name "GetDP/1ResolutionChoices", Visible 0},
C_ = {"-solve -pos -v2", Name "GetDP/9ComputeCommand", Visible 0},
P_ = {"Get_Fields", Name "GetDP/2PostOperationChoices", Visible 0, ReadOnly 1}
R_ = {"WaveEqn_Reso", Name "GetDP/1ResolutionChoices", Visible 0, Choices {"WaveEqn_Reso"} },
P_ = {"WaveEqn2_PostOp", Name "GetDP/2PostOperationChoices", Visible 0, Choices {"WaveEqn2_PostOp"}}
];
AcademicWaves/scattSphere.dat
+ 5
− 5
View file @ dcfd65cc
......@@ -3,6 +3,8 @@
// File: parameters
//========================================================
FLAG_DIM = 3;
FLAG_RES = RES_FREQ;
//---------------------
// Physical parameters
......@@ -26,8 +28,8 @@ z0 = Sqrt[mu0/ep0];
// Incident field
Freq = 1e8;
k = 2*Pi*Freq/c0;
FREQ = 1e8;
k = 2*Pi*FREQ/c0;
// Domain truncation
......@@ -35,13 +37,11 @@ Flag_DomTruncMethod = 9;
// 1: Sommerfeld
// 9: PML
//---------------------
// Numerical parameters
//---------------------
lc = 0.3;
LC = 0.3;
//---------------------
// Geometrical tags
......
AcademicWaves/scattSphere.geo
+ 1
− 2
View file @ dcfd65cc
......@@ -5,13 +5,12 @@
Include "scattSphere.dat";
Mesh.CharacteristicLengthMax = lc;
Mesh.CharacteristicLengthMax = LC;
Mesh.CharacteristicLengthFactor = 1;
Mesh.Optimize = 1;
c0 = newp ; Point(newp) = {0,0,0};
// SPHERE
ps0[]+=newp ; Point(newp) = { R0, 0, 0};
......
AcademicWaves/scattSphere.pro
+ 4
− 4
View file @ dcfd65cc
......@@ -63,9 +63,9 @@ Function {
DampingProfileY[] = (Fabs[Y[]]>=ay) ? c0 * 1/(lpml-(Fabs[Y[]]-ay)) : 0 ;
DampingProfileZ[] = (Fabs[Z[]]>=az) ? c0 * 1/(lpml-(Fabs[Z[]]-az)) : 0 ;
cX[] = Complex[1, -1/(2*Pi*Freq) * DampingProfileX[]] ;
cY[] = Complex[1, -1/(2*Pi*Freq) * DampingProfileY[]] ;
cZ[] = Complex[1, -1/(2*Pi*Freq) * DampingProfileZ[]] ;
cX[] = Complex[1, -1/(2*Pi*FREQ) * DampingProfileX[]] ;
cY[] = Complex[1, -1/(2*Pi*FREQ) * DampingProfileY[]] ;
cZ[] = Complex[1, -1/(2*Pi*FREQ) * DampingProfileZ[]] ;
epsilon[Omega_Pml] = ep0 * TensorDiag[(cY[]*cZ[])/cX[], (cX[]*cZ[])/cY[], (cX[]*cY[])/cZ[]] ;
nu[Omega_Pml] = nu0 * TensorDiag[cX[]/(cY[]*cZ[]), cY[]/(cX[]*cZ[]), cZ[]/(cX[]*cY[])] ;
......@@ -76,7 +76,7 @@ Function {
}
Include "form_vectorWaveEqn_frequency.pro" ;
Include "formulations_vectorWaves.pro" ;
PostProcessing {
{ Name PostPro ; NameOfFormulation Form ;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment