Skip to content
Snippets Groups Projects
Commit b43da713 authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

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

parents b99f89f0 0e5b0d2d
No related branches found
No related tags found
No related merge requests found
Pipeline #4868 passed
Stage: test
Show whitespace changes
Inline Side-by-side
DiffractionGratings/grating2D_data.geo
+ 1
− 1
View file @ b43da713
...@@ -80,7 +80,7 @@ DefineConstant[ ...@@ -80,7 +80,7 @@ DefineConstant[
h_pmltop = {lambda_max, Name StrCat[pp3, "0top PML size [nm]"] , ReadOnly 1, Highlight Str[colorro], ServerAction "Reset GetDP/T0, GetDP/R0, GetDP/Lambda_step, GetDP/total absorption"}, h_pmltop = {lambda_max, Name StrCat[pp3, "0top PML size [nm]"] , ReadOnly 1, Highlight Str[colorro], ServerAction "Reset GetDP/T0, GetDP/R0, GetDP/Lambda_step, GetDP/total absorption"},
h_pmlbot = {lambda_max, Name StrCat[pp3, "1bottom PML size [nm]"] , ReadOnly 1, Highlight Str[colorro], ServerAction "Reset GetDP/T0, GetDP/R0, GetDP/Lambda_step, GetDP/total absorption"}, h_pmlbot = {lambda_max, Name StrCat[pp3, "1bottom PML size [nm]"] , ReadOnly 1, Highlight Str[colorro], ServerAction "Reset GetDP/T0, GetDP/R0, GetDP/Lambda_step, GetDP/total absorption"},
Flag_o2_interp = {1 , Name StrCat[pp3, "2Second order interpolation?"] , Choices {0,1}, ServerAction "Reset GetDP/T0, GetDP/R0, GetDP/Lambda_step, GetDP/total absorption"}, Flag_o2_interp = {1 , Name StrCat[pp3, "2Second order interpolation?"] , Choices {0,1}, ServerAction "Reset GetDP/T0, GetDP/R0, GetDP/Lambda_step, GetDP/total absorption"},
Flag_o2_geom = {0 , Name StrCat[pp3, "3Second order mesh?"] , Choices {0,1}, Visible 0, ServerAction "Reset GetDP/T0, GetDP/R0, GetDP/Lambda_step, GetDP/total absorption"}, Flag_o2_geom = {1 , Name StrCat[pp3, "3Second order mesh?"] , Choices {0,1}, Visible 1, ServerAction "Reset GetDP/T0, GetDP/R0, GetDP/Lambda_step, GetDP/total absorption"},
paramaille = {10 , Name StrCat[pp3, "4nb of mesh elements per wavelength [-]" ] , ReadOnly 0, Highlight Str[colorpp], ServerAction "Reset GetDP/T0, GetDP/R0, GetDP/Lambda_step, GetDP/total absorption"}, paramaille = {10 , Name StrCat[pp3, "4nb of mesh elements per wavelength [-]" ] , ReadOnly 0, Highlight Str[colorpp], ServerAction "Reset GetDP/T0, GetDP/R0, GetDP/Lambda_step, GetDP/total absorption"},
paramaille_scale_sub = {1 , Name StrCat[pp3, "6Custom mesh parameters/refine substrate [-]" ] , ReadOnly 0, Highlight Str[colorpp], Closed 1, ServerAction "Reset GetDP/T0, GetDP/R0, GetDP/Lambda_step, GetDP/total absorption"}, paramaille_scale_sub = {1 , Name StrCat[pp3, "6Custom mesh parameters/refine substrate [-]" ] , ReadOnly 0, Highlight Str[colorpp], Closed 1, ServerAction "Reset GetDP/T0, GetDP/R0, GetDP/Lambda_step, GetDP/total absorption"},
paramaille_scale_layer_dep = {1 , Name StrCat[pp3, "6Custom mesh parameters/refine deposit layer [-]"] , ReadOnly 0, Highlight Str[colorpp], Closed 1, ServerAction "Reset GetDP/T0, GetDP/R0, GetDP/Lambda_step, GetDP/total absorption"}, paramaille_scale_layer_dep = {1 , Name StrCat[pp3, "6Custom mesh parameters/refine deposit layer [-]"] , ReadOnly 0, Highlight Str[colorpp], Closed 1, ServerAction "Reset GetDP/T0, GetDP/R0, GetDP/Lambda_step, GetDP/total absorption"},
......
ElectromagneticScattering/README.md
+ 32
− 7
View file @ b43da713
...@@ -2,22 +2,47 @@ A Onelab model for 3D scattering problems in nanophotonics. ...@@ -2,22 +2,47 @@ A Onelab model for 3D scattering problems in nanophotonics.
## Synopsis ## Synopsis
This project contains a [Onelab](http://onelab.info/wiki/ONELAB) model for solving various 3D electromagnetic scattering problems on an isolated object : This project contains a [Onelab](http://onelab.info/wiki/ONELAB) model for solving various 3D electromagnetic problems on an isolated arbitrary object:
* T-matrix computation: See https://arxiv.org/abs/1802.00596 for details * T-matrix computation [1]
* Quasi-normal modes * Quasi-normal modes
* Plane wave response * Plane wave response
* Green's function (LDOS) * Green's function (e.g. to compute the Local Density Of States)
It contains various usefull features in electromag:
* Vector partial waves (AKA vector spherical harmonics)
* Total field formuation with a point source illumination ("an oriented delta")
* Linear eigenvalue problems
* Scattered field formulation
* Spherical PMLs
## Installation ## Installation
This model requires the following (open-source) programs: This model requires the following programs:
* [onelab](http://onelab.info/wiki/ONELAB), which bundles both [gmsh](http://www.gmsh.info/) and [getdp](http://www.getdp.info/) * [gmsh](http://www.gmsh.info/)
* python (v2.7.x or v3.6.x) with numpy, scipy and matplotlib * [getdp](http://www.getdp.info/) compiled with python support (see below)
* python (>3.5.x) with numpy, scipy and matplotlib
## Running the model ## Running the model
Open `scattering.pro` with Gmsh. Open `scattering.pro` with Gmsh.
The default parameters are set to compute the T-matrix of a sphere. It retrieves the results from [1].
## Authors ## Authors
Guillaume Demésy and Brian Stout Guillaume Demésy and Brian Stout
## References
[1] See "[Scattering matrix of arbitrarily shaped objects: Combining Finite Elements and Vector Partial Waves](https://arxiv.org/abs/1802.00596)" for details about T-matrices
## Installation notes
To enable python support (Python[] function) in GetDP, follow [these instructions (with complex arithmetic)](https://gitlab.onelab.info/getdp/getdp/wikis/GetDP-compilation) and add to the final cmake line:
`-DENABLE_PYTHON=ON -DPYTHON_LIBRARY=/path/to/pythonlib -DPYTHON_INCLUDE_DIR=/path/to/pythoninclude`
* On Debian/Ubuntu systems, for python3.6 installed with apt-get,
* `/path/to/pythonlib` is `/usr/lib/x86_64-linux-gnu/libpython3.6m.so`
* `/path/to/python/include` is `/usr/include/python3.6m`
* For python versions installed through anaconda in some environment (e.g. env py36 below), a common location is:
* `/somepath/anaconda3/envs/py36/lib/libpython3.6m.so`
* `/somepath/anaconda3/envs/py36/include/python3.6m`
NonLinearEVP/NonLinearEVP.geo
+ 132
− 52
View file @ b43da713
...@@ -20,7 +20,6 @@ lc_cell = a_lat/paramaille; ...@@ -20,7 +20,6 @@ lc_cell = a_lat/paramaille;
lc_sq = lc_cell/4; lc_sq = lc_cell/4;
lc_pml = lc_cell*1.2; lc_pml = lc_cell*1.2;
lc_sqa = lc_sq; lc_sqa = lc_sq;
lc_sq_inside = lc_sq*1.4; lc_sq_inside = lc_sq*1.4;
epsc = lc_sq*0.8; epsc = lc_sq*0.8;
...@@ -30,43 +29,37 @@ If (flag_Tmesh==0) ...@@ -30,43 +29,37 @@ If (flag_Tmesh==0)
Point(2) = { a_lat/2. ,-d_sq/2.-space2pml, 0. , lc_cell}; Point(2) = { a_lat/2. ,-d_sq/2.-space2pml, 0. , lc_cell};
Point(3) = { a_lat/2. , d_sq/2.+space2pml, 0. , lc_cell}; Point(3) = { a_lat/2. , d_sq/2.+space2pml, 0. , lc_cell};
Point(4) = {-a_lat/2. , d_sq/2.+space2pml, 0. , lc_cell}; Point(4) = {-a_lat/2. , d_sq/2.+space2pml, 0. , lc_cell};
Point(5) = {-d_sq/2. ,-d_sq/2., 0. , lc_sq};
Point(6) = { d_sq/2. ,-d_sq/2., 0. , lc_sq};
Point(7) = { d_sq/2. , d_sq/2., 0. , lc_sq};
Point(8) = {-d_sq/2. , d_sq/2., 0. , lc_sq};
Point(9) = {-a_lat/2. ,-d_sq/2.-space2pml-pmlsize, 0. , lc_pml}; Point(9) = {-a_lat/2. ,-d_sq/2.-space2pml-pmlsize, 0. , lc_pml};
Point(10) = { a_lat/2. ,-d_sq/2.-space2pml-pmlsize, 0. , lc_pml}; Point(10) = { a_lat/2. ,-d_sq/2.-space2pml-pmlsize, 0. , lc_pml};
Point(11) = { a_lat/2. , d_sq/2.+space2pml+pmlsize, 0. , lc_pml}; Point(11) = { a_lat/2. , d_sq/2.+space2pml+pmlsize, 0. , lc_pml};
Point(12) = {-a_lat/2. , d_sq/2.+space2pml+pmlsize, 0. , lc_pml}; Point(12) = {-a_lat/2. , d_sq/2.+space2pml+pmlsize, 0. , lc_pml};
Point(13) = {-a_lat/2. , a_lat/2., 0. , lc_sqa}; Point(13) = {-a_lat/2. , a_lat/2., 0. , lc_sqa};
Point(14) = {-a_lat/2. ,-a_lat/2., 0. , lc_sqa}; Point(14) = {-a_lat/2. ,-a_lat/2., 0. , lc_sqa};
Point(15) = { a_lat/2. , a_lat/2., 0. , lc_sqa}; Point(15) = { a_lat/2. , a_lat/2., 0. , lc_sqa};
Point(16) = { a_lat/2. ,-a_lat/2., 0. , lc_sqa}; Point(16) = { a_lat/2. ,-a_lat/2., 0. , lc_sqa};
If (flag_rounding==0)
Point(5) = {-d_sq/2. ,-d_sq/2., 0. , lc_sq};
Point(6) = { d_sq/2. ,-d_sq/2., 0. , lc_sq};
Point(7) = { d_sq/2. , d_sq/2., 0. , lc_sq};
Point(8) = {-d_sq/2. , d_sq/2., 0. , lc_sq};
Line(1) = {1, 2}; Line(1) = {1, 2};
Line(3) = {3, 4}; Line(3) = {3, 4};
Line(5) = {5, 6}; Line(5) = {5, 6};
Line(6) = {6, 7}; Line(6) = {6, 7};
Line(7) = {7, 8}; Line(7) = {7, 8};
Line(8) = {8, 5}; Line(8) = {8, 5};
Line(15) = {1, 14}; Line(15) = {1, 14};
Line(16) = {14, 13}; Line(16) = {14, 13};
Line(17) = {13, 4}; Line(17) = {13, 4};
Line(18) = {2, 16}; Line(18) = {2, 16};
Line(19) = {16, 15}; Line(19) = {16, 15};
Line(20) = {15, 3}; Line(20) = {15, 3};
Line(9) = {4, 12}; Line(9) = {4, 12};
Line(10) = {12, 11}; Line(10) = {12, 11};
Line(11) = {11, 3}; Line(11) = {11, 3};
Line(12) = {1, 9}; Line(12) = {1, 9};
Line(13) = {9, 10}; Line(13) = {9, 10};
Line(14) = {10, 2}; Line(14) = {10, 2};
Line Loop(1) = {12, 13, 14, -1}; Line Loop(1) = {12, 13, 14, -1};
Plane Surface(1) = {1}; Plane Surface(1) = {1};
Line Loop(2) = {5, 6, 7, 8}; Line Loop(2) = {5, 6, 7, 8};
...@@ -75,22 +68,104 @@ If (flag_Tmesh==0) ...@@ -75,22 +68,104 @@ If (flag_Tmesh==0)
Plane Surface(3) = {2, -3}; Plane Surface(3) = {2, -3};
Line Loop(4) = {9, 10, 11, 3}; Line Loop(4) = {9, 10, 11, 3};
Plane Surface(4) = {-4}; Plane Surface(4) = {-4};
Periodic Line { 14,18,19,20,11 } = {12,15,16,17,9 } Translate {a_lat,0,0} ; Periodic Line { 14,18,19,20,11 } = {12,15,16,17,9 } Translate {a_lat,0,0} ;
Physical Line("SCATBOUND",1005) = {5,6,7,8}; // bound for lag
Else
// // with a circle
// Point(20) = {0.,0., 0. , lc_sq};
// Point(21) = {-d_sq/2. ,0, 0. , lc_sq};
// Point(22) = {0 ,-d_sq/2., 0. , lc_sq};
// Point(23) = {d_sq/2. ,0, 0. , lc_sq};
// Point(24) = {0 ,d_sq/2., 0. , lc_sq};
// Line(1) = {1, 2};
// Line(3) = {3, 4};
// Line(15) = {1, 14};
// Line(16) = {14, 13};
// Line(17) = {13, 4};
// Line(18) = {2, 16};
// Line(19) = {16, 15};
// Line(20) = {15, 3};
// Line(9) = {4, 12};
// Line(10) = {12, 11};
// Line(11) = {11, 3};
// Line(12) = {1, 9};
// Line(13) = {9, 10};
// Line(14) = {10, 2};
// Circle(21) = {21, 20, 22};
// Circle(22) = {22, 20, 23};
// Circle(23) = {23, 20, 24};
// Circle(24) = {24, 20, 21};
// Line Loop(1) = {12, 13, 14, -1};
// Plane Surface(1) = {1};
// Curve Loop(2) = {21,22,23,24};
// Plane Surface(2) = {2};
// Curve Loop(3) = {20, 3, -17, -16, -15, 1, 18, 19};
// Plane Surface(3) = {2, 3};
// Line Loop(4) = {9, 10, 11, 3};
// Plane Surface(4) = {-4};
// // Rotate {{0, 0, 1}, {0, 0, 0}, 2.*Pi/180.} { Surface{ 2 } ; }
// Periodic Line { 14,18,19,20,11 } = {12,15,16,17,9 } Translate {a_lat,0,0} ;
// Physical Line("SCATBOUND",1005) = {21,22,23,24}; // bound for lag
Point(20) = {-d_sq/2.+corner_rad ,-d_sq/2., 0. , lc_sq};
Point(21) = {-d_sq/2. ,-d_sq/2.+corner_rad, 0. , lc_sq};
Point(22) = {-d_sq/2.+corner_rad ,-d_sq/2.+corner_rad, 0. , lc_sq};
Point(23) = { d_sq/2.-corner_rad ,-d_sq/2., 0. , lc_sq};
Point(24) = { d_sq/2. ,-d_sq/2.+corner_rad, 0. , lc_sq};
Point(25) = { d_sq/2.-corner_rad ,-d_sq/2.+corner_rad, 0. , lc_sq};
Point(26) = { d_sq/2.-corner_rad , d_sq/2., 0. , lc_sq};
Point(27) = { d_sq/2. , d_sq/2.-corner_rad, 0. , lc_sq};
Point(28) = { d_sq/2.-corner_rad , d_sq/2.-corner_rad, 0. , lc_sq};
Point(29) = {-d_sq/2.+corner_rad , d_sq/2., 0. , lc_sq};
Point(30) = {-d_sq/2. , d_sq/2.-corner_rad, 0. , lc_sq};
Point(31) = {-d_sq/2.+corner_rad , d_sq/2.-corner_rad, 0. , lc_sq};
Physical Surface(100) = {2}; // 1 dom in Line(1) = {1, 2};
Physical Surface(101) = {3}; // 2 dom out Line(3) = {3, 4};
Physical Surface(102) = {1}; // PML bot Line(15) = {1, 14};
Physical Surface(103) = {4}; // PML top Line(16) = {14, 13};
Physical Line(1001) = {12,15,16,17,9}; // bloch x left Line(17) = {13, 4};
Physical Line(1002) = {14,18,19,20,11}; // bloch x right Line(18) = {2, 16};
Physical Line(1003) = {10}; // top bound Line(19) = {16, 15};
Physical Line(1004) = {13}; // bot bound Line(20) = {15, 3};
Physical Line(1005) = {5,6,7,8}; // bound for lag Line(9) = {4, 12};
Physical Point(10000) = {1}; // Printpoint Line(10) = {12, 11};
Line(11) = {11, 3};
Line(12) = {1, 9};
Line(13) = {9, 10};
Line(14) = {10, 2};
Circle(21) = {30, 31, 29};
Circle(22) = {26, 28, 27};
Circle(23) = {24, 25, 23};
Circle(24) = {21, 22, 20};
Line(25) = {29, 26};
Line(26) = {27, 24};
Line(27) = {23, 20};
Line(28) = {21, 30};
Line Loop(1) = {12, 13, 14, -1};
Plane Surface(1) = {1};
Curve Loop(2) = {25, 22, 26, 23, 27, -24, 28, 21};
Plane Surface(2) = {2};
Curve Loop(3) = {20, 3, -17, -16, -15, 1, 18, 19};
Plane Surface(3) = {2, 3};
Line Loop(4) = {9, 10, 11, 3};
Plane Surface(4) = {-4};
// Rotate {{0, 0, 1}, {0, 0, 0}, 2.*Pi/180.} { Surface{ 2 } ; }
Periodic Line { 14,18,19,20,11 } = {12,15,16,17,9 } Translate {a_lat,0,0} ;
Physical Line("SCATBOUND",1005) = {21,22,23,24,25,26,27,28}; // bound for lag
EndIf
Physical Surface("SCAT",100) = {2}; // 1 dom in
Physical Surface("OUT",101) = {3}; // 2 dom out
Physical Surface("PMLBOT",102) = {1}; // PML bot
Physical Surface("PMLTOP",103) = {4}; // PML top
Physical Line("BLOCHXL",1001) = {12,15,16,17,9}; // bloch x left
Physical Line("BLOCHXR",1002) = {14,18,19,20,11}; // bloch x right
Physical Line("TOP",1003) = {10}; // top bound
Physical Line("BOT",1004) = {13}; // bot bound
Physical Point("PRINTPOINT",10000) = {1}; // Printpoint
EndIf EndIf
If (flag_Tmesh==1)
If (flag_Tmesh==1)
Point(1) = {-a_lat/2. ,-d_sq/2.-space2pml, 0. , lc_cell}; Point(1) = {-a_lat/2. ,-d_sq/2.-space2pml, 0. , lc_cell};
Point(2) = { a_lat/2. ,-d_sq/2.-space2pml, 0. , lc_cell}; Point(2) = { a_lat/2. ,-d_sq/2.-space2pml, 0. , lc_cell};
Point(3) = { a_lat/2. , d_sq/2.+space2pml, 0. , lc_cell}; Point(3) = { a_lat/2. , d_sq/2.+space2pml, 0. , lc_cell};
...@@ -369,3 +444,8 @@ If (flag_Tmesh==1) ...@@ -369,3 +444,8 @@ If (flag_Tmesh==1)
Physical Point(10000) = {1}; // Printpoint Physical Point(10000) = {1}; // Printpoint
EndIf EndIf
If (flag_o2g==1)
Mesh.ElementOrder = 2;
Else
Mesh.ElementOrder = 1;
EndIf
NonLinearEVP/NonLinearEVP.pro
+ 67
− 33
View file @ b43da713
...@@ -138,7 +138,11 @@ Integration { ...@@ -138,7 +138,11 @@ Integration {
Case { Case {
{ GeoElement Point ; NumberOfPoints 1 ; } { GeoElement Point ; NumberOfPoints 1 ; }
{ GeoElement Line ; NumberOfPoints 4 ; } { GeoElement Line ; NumberOfPoints 4 ; }
{ GeoElement Line2 ; NumberOfPoints 4 ; }
{ GeoElement Line3 ; NumberOfPoints 4 ; }
{ GeoElement Triangle ; NumberOfPoints 6 ; } { GeoElement Triangle ; NumberOfPoints 6 ; }
{ GeoElement Triangle2; NumberOfPoints 12 ; }
{ GeoElement Triangle3; NumberOfPoints 12 ; }
} }
} }
} }
...@@ -148,21 +152,50 @@ Integration { ...@@ -148,21 +152,50 @@ Integration {
FunctionSpace { FunctionSpace {
{ Name Hgrad_perp; Type Form1P; { Name Hgrad_perp; Type Form1P;
BasisFunction { BasisFunction {
{ Name un; NameOfCoef un; Function BF_PerpendicularEdge_1N; Support Region[Om]; Entity NodesOf[All]; } // // for second order mesh elements - Lagrange P3 :
{ Name un2; NameOfCoef un2; Function BF_PerpendicularEdge_2E; Support Region[Om]; Entity EdgesOf[All]; } // un : BF_Node // order 2 is automatic (since all nodes are in NodesOf, even mid-edge ones)
If(flag_o2==1) // un3e : BF_Node_3E // we add part of order 3 basis functions
{ Name un3; NameOfCoef un3; Function BF_PerpendicularEdge_2F; Support Region[Om]; Entity FacetsOf[All]; } // un3f : BF_Node_3F // we add the rest of order 3 basis functions
{ Name un4; NameOfCoef un4; Function BF_PerpendicularEdge_3E; Support Region[Om]; Entity EdgesOf[All]; } { Name sn; NameOfCoef un; Function BF_PerpendicularEdge_1N; Support Region[Om]; Entity NodesOf[All]; }
{ Name un5; NameOfCoef un5; Function BF_PerpendicularEdge_3F; Support Region[Om]; Entity FacetsOf[All]; } If(flag_o2g==0)
// for curved elements, mid-edge nodes are already included in BF_PerpendicularEdge_1N
// But for Mesh.ElementOrder=1, we need to add them explicitely
{ Name sn2e; NameOfCoef un2e; Function BF_PerpendicularEdge_2E; Support Region[Om]; Entity EdgesOf[All]; }
EndIf
If(flag_o2i==1)
{ Name sn3e; NameOfCoef un3e; Function BF_PerpendicularEdge_3E; Support Region[Om]; Entity EdgesOf[All]; }
{ Name sn3f; NameOfCoef un3f; Function BF_PerpendicularEdge_3F; Support Region[Om]; Entity FacetsOf[All]; }
EndIf EndIf
} }
Constraint { Constraint {
{ NameOfCoef un; EntityType NodesOf ; NameOfConstraint BlochX; } { NameOfCoef un; EntityType NodesOf ; NameOfConstraint BlochX; }
{ NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint BlochX; } If(flag_o2g==0)
If(flag_o2==1) { NameOfCoef un2e; EntityType EdgesOf ; NameOfConstraint BlochX; }
{ NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint BlochX; } EndIf
{ NameOfCoef un4; EntityType EdgesOf ; NameOfConstraint BlochX; } If(flag_o2i==1)
{ NameOfCoef un5; EntityType FacetsOf ; NameOfConstraint BlochX; } { NameOfCoef un3f; EntityType EdgesOf ; NameOfConstraint BlochX; }
EndIf
}
}
{ Name Hgrad; Type Form0;
BasisFunction {
{ Name sn; NameOfCoef un; Function BF_Node_1N; Support Region[Om]; Entity NodesOf[All]; }
If(flag_o2g==0) // for curved elements, mid-edge nodes are already included in BF_PerpendicularEdge_1N
{ Name sn2e; NameOfCoef un2e; Function BF_Node_2E; Support Region[Om]; Entity EdgesOf[All]; }
EndIf
If(flag_o2i==1)
{ Name sn3e; NameOfCoef un3e; Function BF_Node_3E; Support Region[Om]; Entity EdgesOf[All]; }
{ Name sn3f; NameOfCoef un3f; Function BF_Node_3F; Support Region[Om]; Entity FacetsOf[All]; }
EndIf
}
Constraint {
{ NameOfCoef un; EntityType NodesOf ; NameOfConstraint BlochX; }
If(flag_o2g==0)
{ NameOfCoef un2e; EntityType EdgesOf ; NameOfConstraint BlochX; }
EndIf
If(flag_o2i==1)
{ NameOfCoef un3f; EntityType EdgesOf ; NameOfConstraint BlochX; }
EndIf EndIf
} }
} }
...@@ -171,7 +204,7 @@ FunctionSpace { ...@@ -171,7 +204,7 @@ FunctionSpace {
BasisFunction { BasisFunction {
{ Name sn; NameOfCoef un; Function BF_Edge ; Support Region[Om]; Entity EdgesOf[All]; } { Name sn; NameOfCoef un; Function BF_Edge ; Support Region[Om]; Entity EdgesOf[All]; }
{ Name sn2; NameOfCoef un2; Function BF_Edge_2E; Support Region[Om]; Entity EdgesOf[All]; } { Name sn2; NameOfCoef un2; Function BF_Edge_2E; Support Region[Om]; Entity EdgesOf[All]; }
If(flag_o2==1) If(flag_o2i==1)
{ Name sn3; NameOfCoef un3; Function BF_Edge_3F_a; Support Region[Om]; Entity FacetsOf[All]; } { Name sn3; NameOfCoef un3; Function BF_Edge_3F_a; Support Region[Om]; Entity FacetsOf[All]; }
{ Name sn4; NameOfCoef un4; Function BF_Edge_3F_b; Support Region[Om]; Entity FacetsOf[All]; } { Name sn4; NameOfCoef un4; Function BF_Edge_3F_b; Support Region[Om]; Entity FacetsOf[All]; }
{ Name sn5; NameOfCoef un5; Function BF_Edge_4E ; Support Region[Om]; Entity EdgesOf[All]; } { Name sn5; NameOfCoef un5; Function BF_Edge_4E ; Support Region[Om]; Entity EdgesOf[All]; }
...@@ -183,7 +216,7 @@ FunctionSpace { ...@@ -183,7 +216,7 @@ FunctionSpace {
{ NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint BlochX; } { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint BlochX; }
{ NameOfCoef un; EntityType EdgesOf ; NameOfConstraint Dirichlet; } { NameOfCoef un; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
{ NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint Dirichlet; } { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
If(flag_o2==1) If(flag_o2i==1)
{ NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint BlochX; } { NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint BlochX; }
{ NameOfCoef un4; EntityType FacetsOf ; NameOfConstraint BlochX; } { NameOfCoef un4; EntityType FacetsOf ; NameOfConstraint BlochX; }
{ NameOfCoef un5; EntityType EdgesOf ; NameOfConstraint BlochX; } { NameOfCoef un5; EntityType EdgesOf ; NameOfConstraint BlochX; }
...@@ -198,7 +231,7 @@ FunctionSpace { ...@@ -198,7 +231,7 @@ FunctionSpace {
BasisFunction { BasisFunction {
{ Name sn; NameOfCoef un; Function BF_Edge ; Support Region[Om_1]; Entity EdgesOf[All]; } { Name sn; NameOfCoef un; Function BF_Edge ; Support Region[Om_1]; Entity EdgesOf[All]; }
{ Name sn2; NameOfCoef un2; Function BF_Edge_2E; Support Region[Om_1]; Entity EdgesOf[All]; } { Name sn2; NameOfCoef un2; Function BF_Edge_2E; Support Region[Om_1]; Entity EdgesOf[All]; }
If(flag_o2==1) If(flag_o2i==1)
{ Name sn3; NameOfCoef un3; Function BF_Edge_3F_a; Support Region[Om_1]; Entity FacetsOf[All]; } { Name sn3; NameOfCoef un3; Function BF_Edge_3F_a; Support Region[Om_1]; Entity FacetsOf[All]; }
{ Name sn4; NameOfCoef un4; Function BF_Edge_3F_b; Support Region[Om_1]; Entity FacetsOf[All]; } { Name sn4; NameOfCoef un4; Function BF_Edge_3F_b; Support Region[Om_1]; Entity FacetsOf[All]; }
{ Name sn5; NameOfCoef un5; Function BF_Edge_4E ; Support Region[Om_1]; Entity EdgesOf[All]; } { Name sn5; NameOfCoef un5; Function BF_Edge_4E ; Support Region[Om_1]; Entity EdgesOf[All]; }
...@@ -211,7 +244,7 @@ FunctionSpace { ...@@ -211,7 +244,7 @@ FunctionSpace {
BasisFunction { BasisFunction {
{ Name sn; NameOfCoef un; Function BF_Edge ; Support Region[{Om_1,Bound}]; Entity EdgesOf[All]; } { Name sn; NameOfCoef un; Function BF_Edge ; Support Region[{Om_1,Bound}]; Entity EdgesOf[All]; }
{ Name sn2; NameOfCoef un2; Function BF_Edge_2E; Support Region[{Om_1,Bound}]; Entity EdgesOf[All]; } { Name sn2; NameOfCoef un2; Function BF_Edge_2E; Support Region[{Om_1,Bound}]; Entity EdgesOf[All]; }
If(flag_o2==1) If(flag_o2i==1)
{ Name sn3; NameOfCoef un3; Function BF_Edge_3F_a; Support Region[{Om_1,Bound}]; Entity FacetsOf[All]; } { Name sn3; NameOfCoef un3; Function BF_Edge_3F_a; Support Region[{Om_1,Bound}]; Entity FacetsOf[All]; }
{ Name sn4; NameOfCoef un4; Function BF_Edge_3F_b; Support Region[{Om_1,Bound}]; Entity FacetsOf[All]; } { Name sn4; NameOfCoef un4; Function BF_Edge_3F_b; Support Region[{Om_1,Bound}]; Entity FacetsOf[All]; }
{ Name sn5; NameOfCoef un5; Function BF_Edge_4E ; Support Region[{Om_1,Bound}]; Entity EdgesOf[All]; } { Name sn5; NameOfCoef un5; Function BF_Edge_4E ; Support Region[{Om_1,Bound}]; Entity EdgesOf[All]; }
...@@ -222,7 +255,7 @@ FunctionSpace { ...@@ -222,7 +255,7 @@ FunctionSpace {
BasisFunction { BasisFunction {
{ Name sn; NameOfCoef un; Function BF_Edge ; Support Region[{Om_2,Bound}]; Entity EdgesOf[All]; } { Name sn; NameOfCoef un; Function BF_Edge ; Support Region[{Om_2,Bound}]; Entity EdgesOf[All]; }
{ Name sn2; NameOfCoef un2; Function BF_Edge_2E; Support Region[{Om_2,Bound}]; Entity EdgesOf[All]; } { Name sn2; NameOfCoef un2; Function BF_Edge_2E; Support Region[{Om_2,Bound}]; Entity EdgesOf[All]; }
If(flag_o2==1) If(flag_o2i==1)
{ Name sn3; NameOfCoef un3; Function BF_Edge_3F_a; Support Region[{Om_2,Bound}]; Entity FacetsOf[All]; } { Name sn3; NameOfCoef un3; Function BF_Edge_3F_a; Support Region[{Om_2,Bound}]; Entity FacetsOf[All]; }
{ Name sn4; NameOfCoef un4; Function BF_Edge_3F_b; Support Region[{Om_2,Bound}]; Entity FacetsOf[All]; } { Name sn4; NameOfCoef un4; Function BF_Edge_3F_b; Support Region[{Om_2,Bound}]; Entity FacetsOf[All]; }
{ Name sn5; NameOfCoef un5; Function BF_Edge_4E ; Support Region[{Om_2,Bound}]; Entity EdgesOf[All]; } { Name sn5; NameOfCoef un5; Function BF_Edge_4E ; Support Region[{Om_2,Bound}]; Entity EdgesOf[All]; }
...@@ -233,7 +266,7 @@ FunctionSpace { ...@@ -233,7 +266,7 @@ FunctionSpace {
{ NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint BlochX; } { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint BlochX; }
{ NameOfCoef un; EntityType EdgesOf ; NameOfConstraint Dirichlet; } { NameOfCoef un; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
{ NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint Dirichlet; } { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
If(flag_o2==1) If(flag_o2i==1)
{ NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint BlochX; } { NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint BlochX; }
{ NameOfCoef un4; EntityType FacetsOf ; NameOfConstraint BlochX; } { NameOfCoef un4; EntityType FacetsOf ; NameOfConstraint BlochX; }
{ NameOfCoef un5; EntityType EdgesOf ; NameOfConstraint BlochX; } { NameOfCoef un5; EntityType EdgesOf ; NameOfConstraint BlochX; }
...@@ -247,7 +280,7 @@ FunctionSpace { ...@@ -247,7 +280,7 @@ FunctionSpace {
BasisFunction { BasisFunction {
{ Name sn; NameOfCoef un; Function BF_Edge ; Support Region[Bound]; Entity EdgesOf[All]; } { Name sn; NameOfCoef un; Function BF_Edge ; Support Region[Bound]; Entity EdgesOf[All]; }
{ Name sn2; NameOfCoef un2; Function BF_Edge_2E; Support Region[Bound]; Entity EdgesOf[All]; } { Name sn2; NameOfCoef un2; Function BF_Edge_2E; Support Region[Bound]; Entity EdgesOf[All]; }
If(flag_o2==1) If(flag_o2i==1)
{ Name sn5; NameOfCoef un5; Function BF_Edge_4E ; Support Region[Bound]; Entity EdgesOf[All]; } { Name sn5; NameOfCoef un5; Function BF_Edge_4E ; Support Region[Bound]; Entity EdgesOf[All]; }
EndIf EndIf
} }
...@@ -319,13 +352,14 @@ Formulation { ...@@ -319,13 +352,14 @@ Formulation {
} }
} }
{Name modal_helmholtz_PEP_h; Type FemEquation; {Name modal_helmholtz_PEP_h; Type FemEquation;
Quantity {{ Name u ; Type Local; NameOfSpace Hgrad_perp ;}} // Quantity {{ Name u ; Type Local; NameOfSpace Hgrad_perp ;}}
Quantity {{ Name u ; Type Local; NameOfSpace Hgrad ;}}
Equation { Equation {
Galerkin { [-1/epsr_nod[]* om_d_1^2 * Dof{Curl u}, {Curl u} ]; In Om_2; Jacobian JSur; Integration Int_1;} Galerkin { [-1/epsr_nod[]* om_d_1^2 * Dof{d u}, {d u} ]; In Om_2; Jacobian JSur; Integration Int_1;}
Galerkin { Eig[ 1/epsr_nod[]*eps_oo_1*gam_1* Dof{Curl u}, {Curl u} ]; Order 1; In Om_2; Jacobian JSur; Integration Int_1;} Galerkin { Eig[ 1/epsr_nod[]*eps_oo_1*gam_1* Dof{d u}, {d u} ]; Order 1; In Om_2; Jacobian JSur; Integration Int_1;}
Galerkin { Eig[-1/epsr_nod[]*eps_oo_1 * Dof{Curl u}, {Curl u} ]; Order 2; In Om_2; Jacobian JSur; Integration Int_1;} Galerkin { Eig[-1/epsr_nod[]*eps_oo_1 * Dof{d u}, {d u} ]; Order 2; In Om_2; Jacobian JSur; Integration Int_1;}
Galerkin { Eig[ 1/epsr_nod[]*gam_1 * Dof{Curl u}, {Curl u} ]; Order 1; In Om_1; Jacobian JSur; Integration Int_1;} Galerkin { Eig[ 1/epsr_nod[]*gam_1 * Dof{d u}, {d u} ]; Order 1; In Om_1; Jacobian JSur; Integration Int_1;}
Galerkin { Eig[-1/epsr_nod[] * Dof{Curl u}, {Curl u} ]; Order 2; In Om_1; Jacobian JSur; Integration Int_1;} Galerkin { Eig[-1/epsr_nod[] * Dof{d u}, {d u} ]; Order 2; In Om_1; Jacobian JSur; Integration Int_1;}
Galerkin { Eig[-mur[]/cel^2 *om_d_1^2 * Dof{u}, {u} ]; Order 2; In Om ; Jacobian JSur; Integration Int_1;} Galerkin { Eig[-mur[]/cel^2 *om_d_1^2 * Dof{u}, {u} ]; Order 2; In Om ; Jacobian JSur; Integration Int_1;}
Galerkin { Eig[ mur[]/cel^2 *eps_oo_1*gam_1* Dof{u}, {u} ]; Order 3; In Om ; Jacobian JSur; Integration Int_1;} Galerkin { Eig[ mur[]/cel^2 *eps_oo_1*gam_1* Dof{u}, {u} ]; Order 3; In Om ; Jacobian JSur; Integration Int_1;}
Galerkin { Eig[-mur[]/cel^2 *eps_oo_1 * Dof{u}, {u} ]; Order 4; In Om ; Jacobian JSur; Integration Int_1;} Galerkin { Eig[-mur[]/cel^2 *eps_oo_1 * Dof{u}, {u} ]; Order 4; In Om ; Jacobian JSur; Integration Int_1;}
...@@ -371,7 +405,7 @@ Resolution { ...@@ -371,7 +405,7 @@ Resolution {
} }
Operation{ Operation{
GenerateSeparate[M1]; GenerateSeparate[M1];
EigenSolve[M1,neig,eig_target_re,eig_target_im]; // Print[M1]; EigenSolve[M1,neig,eig_target_re,eig_target_im]; Print[M1];
} }
} }
...@@ -482,7 +516,7 @@ PostOperation { ...@@ -482,7 +516,7 @@ PostOperation {
eig_tol = 1e-6; eig_tol = 1e-6;
slepc_options_nep_rat = Sprintf(" -nep_view "); slepc_options_nep_rat = Sprintf(" -nep_view ");
// // These slepc options are now by default // // These slepc options are now by default (see Kernel/EigenSolve_SLEPC.cpp for default options)
// slepc_options_st = " -st_type sinvert -st_ksp_type preonly -st_pc_type lu -st_pc_factor_mat_solver_type mumps "; // slepc_options_st = " -st_type sinvert -st_ksp_type preonly -st_pc_type lu -st_pc_factor_mat_solver_type mumps ";
// slepc_options_pep = Sprintf(" -pep_max_it 400 -pep_target_magnitude -pep_tol %.10f -pep_view ",eig_tol); // slepc_options_pep = Sprintf(" -pep_max_it 400 -pep_target_magnitude -pep_tol %.10f -pep_view ",eig_tol);
slepc_options_st = " "; slepc_options_st = " ";
......
NonLinearEVP/NonLinearEVP_data.geo
+ 8
− 4
View file @ b43da713
...@@ -34,7 +34,7 @@ DefineConstant[ a_lat = {50 , Name StrCat[pp0 , "1grating period d [nm]"] ...@@ -34,7 +34,7 @@ DefineConstant[ a_lat = {50 , Name StrCat[pp0 , "1grating period d [nm]"]
DefineConstant[ DefineConstant[
d_sq = {0.806 , Name StrCat[pp0 , "2sq [d]"] , Highlight Str[colorpp0] , Min 0.01, Max 0.99} , d_sq = {0.806 , Name StrCat[pp0 , "2square side size (fraction of period) [-]"] , Highlight Str[colorpp0] , Min 0.01, Max 0.99} ,
space2pml = {1 , Name StrCat[pp0 , "3PML distance to object [d]"] , Highlight Str[colorpp0], Min 0.1, Max 2} , space2pml = {1 , Name StrCat[pp0 , "3PML distance to object [d]"] , Highlight Str[colorpp0], Min 0.1, Max 2} ,
pmlsize = {3 , Name StrCat[pp0 , "4PML thickness [d]"] , Highlight Str[colorpp0], Min 0.5, Max 4.} , pmlsize = {3 , Name StrCat[pp0 , "4PML thickness [d]"] , Highlight Str[colorpp0], Min 0.5, Max 4.} ,
...@@ -55,12 +55,15 @@ DefineConstant[ ...@@ -55,12 +55,15 @@ DefineConstant[
paramaille = {4 , Name StrCat[pp4 , "0number of mesh elements per period []"] , Highlight Str[colorpp4], Min 2, Max 10} , paramaille = {4 , Name StrCat[pp4 , "0number of mesh elements per period []"] , Highlight Str[colorpp4], Min 2, Max 10} ,
flag_Tmesh = {0 , Name StrCat[pp4 , "2locally structured mesh?"] , Choices {0="unstruct",1="struct"} }, flag_Tmesh = {0 , Name StrCat[pp4 , "2locally structured mesh?"] , Choices {0="unstruct",1="struct"} },
flag_o2 = {1 , Name StrCat[pp4 , "3FEM order"] , Choices {0="o1",1="o2"} }, flag_o2g = {1 , Name StrCat[pp4 , "3Geometrical order"] , Choices {0="order 1 (linear)",1="order 2 (curved)"} },
flag_o2i = {1 , Name StrCat[pp4 , "4Interpolation order"] , Choices {0="order 1",1="full order 2"}, ServerAction "ResetDatabase"},
flag_rounding = {1 , Name StrCat[pp4 , "5Corner rounding/0Enable rounding of corners!"], Choices{0,1}, ServerAction "ResetDatabase"},
corner_rad_frac= {0.05 , Name StrCat[pp4 , "5Corner rounding/1corner radius (fraction of square side)"] , Highlight Str[colorpp2] , Min 0.005, Max 0.49},
flag_outEigvec = {1 , Name StrCat[pp4 , "7output eigenvector?"], Choices{0,1}},
flag_res = {2 , Name StrCat[pp5 , "0resolution type"], flag_res = {2 , Name StrCat[pp5 , "0resolution type"],
// Choices {0="Aux_E" ,1="PEP_E" ,2="NEP_E" ,3="Lag_E" ,4="PEP_h", 5="all"},ServerAction "ResetDatabase"}, // Choices {0="Aux_E" ,1="PEP_E" ,2="NEP_E" ,3="Lag_E" ,4="PEP_h", 5="all"},ServerAction "ResetDatabase"},
Choices {0="Aux_E" ,1="PEP_E" ,2="NEP_E" ,3="Lag_E" ,4="PEP_h"},ServerAction "ResetDatabase"}, Choices {0="Aux_E" ,1="PEP_E" ,2="NEP_E" ,3="Lag_E" ,4="PEP_h"},ServerAction "ResetDatabase"}
flag_outEigvec = {1 , Name StrCat[pp4, "output eigenvector?"], Choices{0,1}}
]; ];
// Normalized units so that 2*pi*c/a=1 // Normalized units so that 2*pi*c/a=1
...@@ -85,6 +88,7 @@ eig_min_im = eig_min_im / norm; ...@@ -85,6 +88,7 @@ eig_min_im = eig_min_im / norm;
eig_max_im = eig_max_im / norm; eig_max_im = eig_max_im / norm;
om_d_1 = om_d_1 / norm; om_d_1 = om_d_1 / norm;
gam_1 = gam_1 / norm; gam_1 = gam_1 / norm;
corner_rad = d_sq*corner_rad_frac;
eps_oo_2 = 1; eps_oo_2 = 1;
om_d_2 = 0; om_d_2 = 0;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment