diff --git a/DiffractionGratings/grating2D_data.geo b/DiffractionGratings/grating2D_data.geo
index 5f894f18119956fccf9a25ef48e20c6cc969cd40..5bf7f20e07571fb236b1da6fa009842c4b4b5c47 100644
--- a/DiffractionGratings/grating2D_data.geo
+++ b/DiffractionGratings/grating2D_data.geo
@@ -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_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_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_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"},
diff --git a/NonLinearEVP/NonLinearEVP.geo b/NonLinearEVP/NonLinearEVP.geo
index 0bc4185eb6a55eb815c5aafbcf4fa527b5d66e52..4a635fc031a7126f4ac7f5a7805aceab1cb6682a 100644
--- a/NonLinearEVP/NonLinearEVP.geo
+++ b/NonLinearEVP/NonLinearEVP.geo
@@ -20,7 +20,6 @@ lc_cell = a_lat/paramaille;
lc_sq = lc_cell/4;
lc_pml = lc_cell*1.2;
lc_sqa = lc_sq;
-
lc_sq_inside = lc_sq*1.4;
epsc = lc_sq*0.8;
@@ -30,67 +29,143 @@ If (flag_Tmesh==0)
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(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(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(12) = {-a_lat/2. , d_sq/2.+space2pml+pmlsize, 0. , lc_pml};
-
Point(13) = {-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(16) = { a_lat/2. ,-a_lat/2., 0. , lc_sqa};
-
- Line(1) = {1, 2};
- Line(3) = {3, 4};
- Line(5) = {5, 6};
- Line(6) = {6, 7};
- Line(7) = {7, 8};
- Line(8) = {8, 5};
-
- 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};
-
- Line Loop(1) = {12, 13, 14, -1};
- Plane Surface(1) = {1};
- Line Loop(2) = {5, 6, 7, 8};
- Plane Surface(2) = {2};
- Line Loop(3) = {17, -3, -20, -19, -18, -1, 15, 16};
- Plane Surface(3) = {2, -3};
- Line Loop(4) = {9, 10, 11, 3};
- Plane Surface(4) = {-4};
-
- Periodic Line { 14,18,19,20,11 } = {12,15,16,17,9 } Translate {a_lat,0,0} ;
-
- Physical Surface(100) = {2}; // 1 dom in
- Physical Surface(101) = {3}; // 2 dom out
- Physical Surface(102) = {1}; // PML bot
- Physical Surface(103) = {4}; // PML top
- Physical Line(1001) = {12,15,16,17,9}; // bloch x left
- Physical Line(1002) = {14,18,19,20,11}; // bloch x right
- Physical Line(1003) = {10}; // top bound
- Physical Line(1004) = {13}; // bot bound
- Physical Line(1005) = {5,6,7,8}; // bound for lag
- Physical Point(10000) = {1}; // Printpoint
+ 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(3) = {3, 4};
+ Line(5) = {5, 6};
+ Line(6) = {6, 7};
+ Line(7) = {7, 8};
+ Line(8) = {8, 5};
+ 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};
+ Line Loop(1) = {12, 13, 14, -1};
+ Plane Surface(1) = {1};
+ Line Loop(2) = {5, 6, 7, 8};
+ Plane Surface(2) = {2};
+ Line Loop(3) = {17, -3, -20, -19, -18, -1, 15, 16};
+ Plane Surface(3) = {2, -3};
+ Line Loop(4) = {9, 10, 11, 3};
+ Plane Surface(4) = {-4};
+ 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};
+
+ 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) = {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
-If (flag_Tmesh==1)
+If (flag_Tmesh==1)
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(3) = { a_lat/2. , d_sq/2.+space2pml, 0. , lc_cell};
@@ -369,3 +444,8 @@ If (flag_Tmesh==1)
Physical Point(10000) = {1}; // Printpoint
EndIf
+If (flag_o2g==1)
+ Mesh.ElementOrder = 2;
+Else
+ Mesh.ElementOrder = 1;
+EndIf
diff --git a/NonLinearEVP/NonLinearEVP.pro b/NonLinearEVP/NonLinearEVP.pro
index aeb360842e74c8af632218186ab489a3773c811d..a27671f8dd7941f2ec6dffbe617b76a2bd0f801a 100644
--- a/NonLinearEVP/NonLinearEVP.pro
+++ b/NonLinearEVP/NonLinearEVP.pro
@@ -136,9 +136,13 @@ Integration {
Case {
{ Type Gauss ;
Case {
- { GeoElement Point ; NumberOfPoints 1 ; }
- { GeoElement Line ; NumberOfPoints 4 ; }
- { GeoElement Triangle ; NumberOfPoints 6 ; }
+ { GeoElement Point ; NumberOfPoints 1 ; }
+ { GeoElement Line ; NumberOfPoints 4 ; }
+ { GeoElement Line2 ; NumberOfPoints 4 ; }
+ { GeoElement Line3 ; NumberOfPoints 4 ; }
+ { GeoElement Triangle ; NumberOfPoints 6 ; }
+ { GeoElement Triangle2; NumberOfPoints 12 ; }
+ { GeoElement Triangle3; NumberOfPoints 12 ; }
}
}
}
@@ -146,23 +150,52 @@ Integration {
}
FunctionSpace {
- { Name Hgrad_perp; Type Form1P;
+ { Name Hgrad_perp; Type Form1P;
BasisFunction {
- { Name un; NameOfCoef un; Function BF_PerpendicularEdge_1N; Support Region[Om]; Entity NodesOf[All]; }
- { Name un2; NameOfCoef un2; Function BF_PerpendicularEdge_2E; Support Region[Om]; Entity EdgesOf[All]; }
- If(flag_o2==1)
- { Name un3; NameOfCoef un3; Function BF_PerpendicularEdge_2F; Support Region[Om]; Entity FacetsOf[All]; }
- { Name un4; NameOfCoef un4; Function BF_PerpendicularEdge_3E; Support Region[Om]; Entity EdgesOf[All]; }
- { Name un5; NameOfCoef un5; Function BF_PerpendicularEdge_3F; Support Region[Om]; Entity FacetsOf[All]; }
+ // // for second order mesh elements - Lagrange P3 :
+ // un : BF_Node // order 2 is automatic (since all nodes are in NodesOf, even mid-edge ones)
+ // un3e : BF_Node_3E // we add part of order 3 basis functions
+ // un3f : BF_Node_3F // we add the rest of order 3 basis functions
+ { Name sn; NameOfCoef un; Function BF_PerpendicularEdge_1N; Support Region[Om]; Entity NodesOf[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
}
Constraint {
{ NameOfCoef un; EntityType NodesOf ; NameOfConstraint BlochX; }
- { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint BlochX; }
- If(flag_o2==1)
- { NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint BlochX; }
- { NameOfCoef un4; EntityType EdgesOf ; NameOfConstraint BlochX; }
- { NameOfCoef un5; EntityType FacetsOf ; NameOfConstraint BlochX; }
+ If(flag_o2g==0)
+ { NameOfCoef un2e; EntityType EdgesOf ; NameOfConstraint BlochX; }
+ EndIf
+ If(flag_o2i==1)
+ { 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
}
}
@@ -171,7 +204,7 @@ FunctionSpace {
BasisFunction {
{ 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]; }
- If(flag_o2==1)
+ If(flag_o2i==1)
{ 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 sn5; NameOfCoef un5; Function BF_Edge_4E ; Support Region[Om]; Entity EdgesOf[All]; }
@@ -183,7 +216,7 @@ FunctionSpace {
{ NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint BlochX; }
{ NameOfCoef un; 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 un4; EntityType FacetsOf ; NameOfConstraint BlochX; }
{ NameOfCoef un5; EntityType EdgesOf ; NameOfConstraint BlochX; }
@@ -198,7 +231,7 @@ FunctionSpace {
BasisFunction {
{ 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]; }
- 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 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]; }
@@ -211,7 +244,7 @@ FunctionSpace {
BasisFunction {
{ 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]; }
- 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 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]; }
@@ -222,7 +255,7 @@ FunctionSpace {
BasisFunction {
{ 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]; }
- 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 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]; }
@@ -233,7 +266,7 @@ FunctionSpace {
{ NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint BlochX; }
{ NameOfCoef un; 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 un4; EntityType FacetsOf ; NameOfConstraint BlochX; }
{ NameOfCoef un5; EntityType EdgesOf ; NameOfConstraint BlochX; }
@@ -247,7 +280,7 @@ FunctionSpace {
BasisFunction {
{ 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]; }
- If(flag_o2==1)
+ If(flag_o2i==1)
{ Name sn5; NameOfCoef un5; Function BF_Edge_4E ; Support Region[Bound]; Entity EdgesOf[All]; }
EndIf
}
@@ -319,16 +352,17 @@ Formulation {
}
}
{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 {
- Galerkin { [-1/epsr_nod[]* om_d_1^2 * Dof{Curl u}, {Curl 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 * Dof{Curl u}, {Curl 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[] * Dof{Curl u}, {Curl 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 *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 { [-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{d u}, {d u} ]; Order 1; 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{d u}, {d u} ]; Order 1; 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 *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;}
}
}
}
@@ -371,7 +405,7 @@ Resolution {
}
Operation{
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 {
eig_tol = 1e-6;
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_pep = Sprintf(" -pep_max_it 400 -pep_target_magnitude -pep_tol %.10f -pep_view ",eig_tol);
slepc_options_st = " ";
diff --git a/NonLinearEVP/NonLinearEVP_data.geo b/NonLinearEVP/NonLinearEVP_data.geo
index ea5e2cf0e38e665a31cea22d3c660763c28bf658..b66383ce830a3ba123fd456f4b62d5eee161b39d 100644
--- a/NonLinearEVP/NonLinearEVP_data.geo
+++ b/NonLinearEVP/NonLinearEVP_data.geo
@@ -34,7 +34,7 @@ DefineConstant[ a_lat = {50 , Name StrCat[pp0 , "1grating period d [nm]"]
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} ,
pmlsize = {3 , Name StrCat[pp0 , "4PML thickness [d]"] , Highlight Str[colorpp0], Min 0.5, Max 4.} ,
@@ -55,12 +55,15 @@ DefineConstant[
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_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"],
// 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"},
- flag_outEigvec = {1 , Name StrCat[pp4, "output eigenvector?"], Choices{0,1}}
+ Choices {0="Aux_E" ,1="PEP_E" ,2="NEP_E" ,3="Lag_E" ,4="PEP_h"},ServerAction "ResetDatabase"}
];
// Normalized units so that 2*pi*c/a=1
@@ -85,6 +88,7 @@ eig_min_im = eig_min_im / norm;
eig_max_im = eig_max_im / norm;
om_d_1 = om_d_1 / norm;
gam_1 = gam_1 / norm;
+corner_rad = d_sq*corner_rad_frac;
eps_oo_2 = 1;
om_d_2 = 0;