diff --git a/NonLinearEVP/NonLinearEVP.geo b/NonLinearEVP/NonLinearEVP.geo
index 91eb52d22e230cb58d42399241096a0228010e13..6ad794673f3699239a0c489b65c4eed8143d0b18 100644
--- a/NonLinearEVP/NonLinearEVP.geo
+++ b/NonLinearEVP/NonLinearEVP.geo
@@ -16,135 +16,404 @@
// along with This program. If not, see <https://www.gnu.org/licenses/>.
Include "NonLinearEVP_data.geo";
-lc_cell = a_lat/paramaille;
-lc_sq = lc_cell/6;
-lc_pml = lc_cell*1.2;
-lc_sqa = lc_sq;
+lc_cell = a_lat/paramaille;
+lc_sq = lc_cell/2;
+lc_pml = lc_cell*1.2;
+lc_sqa = lc_sq;
lc_sq_inside = lc_sq*1.4;
epsc = lc_sq*0.9;
+refinec = 1;
+If (flag_geom==1)
+ If (flag_round==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};
+ Point(4) = {-a_lat/2. , d_sq/2.+space2pml, 0. , lc_cell};
+ Point(5) = {-d_sq/2.+corner_rad ,-d_sq/2. , 0. , lc_sq/refinec};
+ Point(6) = {-d_sq/2. ,-d_sq/2.+corner_rad, 0. , lc_sq/refinec};
+ Point(8) = { d_sq/2.-corner_rad ,-d_sq/2. , 0. , lc_sq/refinec};
+ Point(9) = { d_sq/2. ,-d_sq/2.+corner_rad, 0. , lc_sq/refinec};
+ Point(11) = { d_sq/2.-corner_rad , d_sq/2. , 0. , lc_sq/refinec};
+ Point(12) = { d_sq/2. , d_sq/2.-corner_rad, 0. , lc_sq/refinec};
+ Point(14) = {-d_sq/2.+corner_rad , d_sq/2. , 0. , lc_sq/refinec};
+ Point(15) = {-d_sq/2. , d_sq/2.-corner_rad, 0. , lc_sq/refinec};
+ // center
+ Point(7) = {-d_sq/2.+corner_rad ,-d_sq/2.+corner_rad, 0. , lc_sq/(refinec)};
+ Point(10) = { d_sq/2.-corner_rad ,-d_sq/2.+corner_rad, 0. , lc_sq/(refinec)};
+ Point(13) = { d_sq/2.-corner_rad , d_sq/2.-corner_rad, 0. , lc_sq/(refinec)};
+ Point(16) = {-d_sq/2.+corner_rad , d_sq/2.-corner_rad, 0. , lc_sq/(refinec)};
+
+ Point(17) = {-a_lat/2. ,-d_sq/2.-space2pml-pmlsize, 0. , lc_pml};
+ Point(18) = { a_lat/2. ,-d_sq/2.-space2pml-pmlsize, 0. , lc_pml};
+ Point(19) = { a_lat/2. , d_sq/2.+space2pml+pmlsize, 0. , lc_pml};
+ Point(20) = {-a_lat/2. , d_sq/2.+space2pml+pmlsize, 0. , lc_pml};
+
+ Point(21) = {-(d_sq/2.+(Sqrt[2]/2-1)*corner_rad) , -(d_sq/2.+(Sqrt[2]/2-1)*corner_rad), 0. , lc_sq/refinec};
+ Point(22) = { (d_sq/2.+(Sqrt[2]/2-1)*corner_rad) , -(d_sq/2.+(Sqrt[2]/2-1)*corner_rad), 0. , lc_sq/refinec};
+ Point(23) = { (d_sq/2.+(Sqrt[2]/2-1)*corner_rad) , (d_sq/2.+(Sqrt[2]/2-1)*corner_rad), 0. , lc_sq/refinec};
+ Point(24) = {-(d_sq/2.+(Sqrt[2]/2-1)*corner_rad) , (d_sq/2.+(Sqrt[2]/2-1)*corner_rad), 0. , lc_sq/refinec};
+
+ Point(25) = {-(d_sq/2) , 0, 0. , lc_sq};
+ Point(26) = { (d_sq/2) , 0, 0. , lc_sq};
+ Point(27) = { 0 , (d_sq/2.), 0. , lc_sq};
+ Point(28) = { 0 , (-d_sq/2.), 0. , lc_sq};
+ Point(29) = {0,0 , 0. , lc_sq};
+
+ Circle(15) = {15, 16, 24};
+ Circle(16) = {24, 16, 14};
+ Circle(17) = {11, 13, 23};
+ Circle(18) = {23, 13, 12};
+ Circle(19) = {9, 10, 22};
+ Circle(20) = {22, 10, 8};
+ Circle(21) = {5, 7, 21};
+ Circle(22) = {21, 7, 6};
+
+ Line(1) = {20, 19};
+ Line(2) = {4, 3};
+ Line(3) = {1, 2};
+ Line(4) = {17, 18};
+ Line(5) = {17, 1};
+ Line(6) = {1, 4};
+ Line(7) = {4, 20};
+ Line(8) = {18, 2};
+ Line(9) = {2, 3};
+ Line(10) = {3, 19};
+ Line(23) = {14, 27};
+ Line(24) = {27, 11};
+ Line(25) = {12, 26};
+ Line(26) = {26, 9};
+ Line(27) = {8, 28};
+ Line(28) = {28, 5};
+ Line(29) = {6, 25};
+ Line(30) = {25, 15};
+ Line Loop(1) = {7, 1, -10, -2};
+ Plane Surface(1) = {1};
+ Curve Loop(2) = {30, 15, 16, 23, 24, 17, 18, 25, 26, 19, 20, 27, 28, 21, 22, 29};
+ Plane Surface(2) = {2};
+ Line Loop(3) = {3, 9, -2, -6};
+ Plane Surface(3) = {2, -3};
+ Line Loop(4) = {5, 3, -8, -4};
+ Plane Surface(4) = {4};
+ Point{7,10,13,16,29} In Surface{2};
+
+ Periodic Line {8,9,10} = {5,6,7} Translate {a_lat,0,0} ;
+
+ Physical Surface(100) = {2}; // 1 dom in
+ Physical Surface(101) = {3}; // 2 dom out
+ Physical Surface(102) = {1}; // PML top
+ Physical Surface(103) = {4}; // PML bot
+ Physical Line(1001) = {5,6,7}; // bloch x left
+ Physical Line(1002) = {8,9,10}; // bloch x right
+ Physical Line(1003) = {1}; // top bound
+ Physical Line(1004) = {4}; // bot bound
+ Physical Line(1005) = {11,12,13,14,15,16,17,18}; // bound for lag
+ Physical Point(10000) = {1}; // Printpoint
+
+ Else
+ 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};
+ 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};
+
+ Line(1) = {1, 2};
+ Line(2) = {2, 3};
+ Line(3) = {3, 4};
+ Line(4) = {4, 1};
+ Line(5) = {5, 6};
+ Line(6) = {6, 7};
+ Line(7) = {7, 8};
+ Line(8) = {8, 5};
+
+ 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) = {4, 1, 2, 3};
+ Plane Surface(3) = {2, -3};
+ Line Loop(4) = {9, 10, 11, 3};
+ Plane Surface(4) = {-4};
+
+ // Periodic Line { 14,2,11} = { 12,4,9} ;
+ Periodic Line { 14,2,11 } = {12,4,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 top
+ Physical Surface(103) = {4}; // PML bot
+ Physical Line(1001) = {12,4,9}; // bloch x left
+ Physical Line(1002) = {14,2,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
+
+ EndIf
+EndIf
+
+If (flag_geom==2)
-If (flag_Tmesh==0)
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};
Point(4) = {-a_lat/2. , d_sq/2.+space2pml, 0. , lc_cell};
+
+ Point(5) = {-a_lat/2. ,-d_sq/2, 0. , lc_sq/refinec};
+ Point(6) = { a_lat/2. ,-d_sq/2, 0. , lc_sq/refinec};
+ Point(7) = { a_lat/2. , d_sq/2, 0. , lc_sq/refinec};
+ Point(8) = {-a_lat/2. , d_sq/2, 0. , lc_sq/refinec};
+
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};
- 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
- 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
+
+ Point(13) = {-a_lat/2. ,-d_sq/2-space2pml/3, 0. , lc_cell};
+ Point(14) = { a_lat/2. ,-d_sq/2-space2pml/3, 0. , lc_cell};
+ Point(15) = { a_lat/2. , d_sq/2+space2pml/3, 0. , lc_cell};
+ Point(16) = {-a_lat/2. , d_sq/2+space2pml/3, 0. , lc_cell};
+ Point(17) = {-a_lat/2. ,-d_sq/2+d_sq/7, 0. , lc_sq};
+ Point(18) = { a_lat/2. ,-d_sq/2+d_sq/7, 0. , lc_sq};
+ Point(19) = { a_lat/2. , d_sq/2-d_sq/7, 0. , lc_sq};
+ Point(20) = {-a_lat/2. , d_sq/2-d_sq/7, 0. , lc_sq};
+ Point(21) = { 0. , d_sq/2-d_sq/7, 0. , lc_sq};
+ Point(22) = { 0. ,-d_sq/2+d_sq/7, 0. , lc_sq};
+
+ Line(1) = {9, 1};
+ Line(2) = {1, 13};
+ Line(3) = {13, 5};
+ Line(4) = {5, 17};
+ Line(5) = {17, 20};
+ Line(6) = {20, 8};
+ Line(7) = {8, 16};
+ Line(8) = {16, 4};
+ Line(9) = {4, 12};
+ Line(10) = {10, 2};
+ Line(11) = {2, 14};
+ Line(12) = {14, 6};
+ Line(13) = {6, 18};
+ Line(14) = {18, 19};
+ Line(15) = {19, 7};
+ Line(16) = {7, 15};
+ Line(17) = {15, 3};
+ Line(18) = {3, 11};
+ Line(19) = {9, 10};
+ Line(20) = {1, 2};
+ Line(21) = {5, 6};
+ Line(22) = {8, 7};
+ Line(23) = {4, 3};
+ Line(24) = {12, 11};
+ Curve Loop(1) = {1, 20, -10, -19};
+ Plane Surface(1) = {1};
+ Curve Loop(2) = {2, 3, 21, -12, -11, -20};
+ Plane Surface(2) = {2};
+ Curve Loop(3) = {5, 6, 22, -15, -14, -13, -21, 4};
+ Plane Surface(3) = {3};
+ Curve Loop(4) = {8, 23, -17, -16, -22, 7};
+ Plane Surface(4) = {4};
+ Curve Loop(5) = {23, 18, -24, -9};
+ Plane Surface(5) = {5};
+
+ Periodic Line { 10,11,12,13,14,15,16,17,18 } = { 1,2,3,4,5,6,7,8,9 } Translate {a_lat,0,0} ;
+
+ Point{21,22} In Surface{3};
+
+ Physical Surface('IN',100) = {3}; // 1 dom in
+ Physical Surface('OUT',101) = {2,4}; // 2 dom out
+ Physical Surface('PMLTOP',102) = {5}; // PML top
+ Physical Surface('PMLBOT',103) = {1}; // PML bot
+ Physical Line('BXM',1001) = {1,2,3,4,5,6,7,8,9}; // bloch x left
+ Physical Line('BXP',1002) = {10,11,12,13,14,15,16,17,18}; // bloch x right
+ Physical Line(1003) = {24}; // top bound
+ Physical Line(1004) = {19}; // bot bound
+ Physical Line(1005) = {21,22}; // bound for lag
+ Physical Point(10000) = {1}; // Printpoint
+
+EndIf
+
+If (flag_geom==3)
+ 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};
+ Point(4) = {-a_lat/2. , d_sq/2.+space2pml, 0. , lc_cell};
+
+ Point(17) = {-a_lat/2. ,-d_sq/2.-space2pml-pmlsize, 0. , lc_pml};
+ Point(18) = { a_lat/2. ,-d_sq/2.-space2pml-pmlsize, 0. , lc_pml};
+ Point(19) = { a_lat/2. , d_sq/2.+space2pml+pmlsize, 0. , lc_pml};
+ Point(20) = {-a_lat/2. , d_sq/2.+space2pml+pmlsize, 0. , lc_pml};
+
+ Point(25) = {-(d_sq/2) , 0, 0. , lc_sq/refinec};
+ Point(26) = { (d_sq/2) , 0, 0. , lc_sq/refinec};
+ Point(27) = { 0 , (d_sq/2.), 0. , lc_sq/refinec};
+ Point(28) = { 0 , (-d_sq/2.), 0. , lc_sq/refinec};
+ Point(29) = {0,0 , 0. , lc_sq};
+
+
+ Line(1) = {20, 19};
+ Line(2) = {4, 3};
+ Line(3) = {1, 2};
+ Line(4) = {17, 18};
+ Line(5) = {17, 1};
+ Line(6) = {1, 4};
+ Line(7) = {4, 20};
+ Line(8) = {18, 2};
+ Line(9) = {2, 3};
+ Line(10) = {3, 19};
+ Circle(11) = {25, 29, 27};
+ Circle(12) = {27, 29, 26};
+ Circle(13) = {26, 29, 28};
+ Circle(14) = {28, 29, 25};
+
+ Line Loop(1) = {7, 1, -10, -2};
+ Plane Surface(1) = {1};
+ Curve Loop(2) = {11, 12, 13, 14};
+ Plane Surface(2) = {2};
+ Curve Loop(3) = {6, 2, -9, -3};
+ Plane Surface(3) = {2, 3};
+ Line Loop(4) = {5, 3, -8, -4};
+ Plane Surface(4) = {4};
+
+ Periodic Line {8,9,10} = {5,6,7} Translate {a_lat,0,0} ;
+
+ Physical Surface(100) = {2}; // 1 dom in
+ Physical Surface(101) = {3}; // 2 dom out
+ Physical Surface(102) = {1}; // PML top
+ Physical Surface(103) = {4}; // PML bot
+ Physical Line(1001) = {5,6,7}; // bloch x left
+ Physical Line(1002) = {8,9,10}; // bloch x right
+ Physical Line(1003) = {1}; // top bound
+ Physical Line(1004) = {4}; // bot bound
+ Physical Line(1005) = {11,12,13,14}; // bound for lag
+ Physical Point(10000) = {1}; // Printpoint
EndIf
-If (flag_Tmesh==1)
+If (flag_geom==4)
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};
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(5) = {-a_lat/2. ,-d_sq/2, 0. , lc_sq/refinec};
+ Point(6) = { a_lat/2. ,-d_sq/2, 0. , lc_sq/refinec};
+ Point(7) = { a_lat/2. , d_sq/2, 0. , lc_sq/refinec};
+ Point(8) = {-a_lat/2. , d_sq/2, 0. , lc_sq/refinec};
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. ,-d_sq/2-1.*lc_sq, 0. , lc_sq};
+ Point(14) = { a_lat/2. ,-d_sq/2-1.*lc_sq, 0. , lc_sq};
+ Point(15) = { a_lat/2. , d_sq/2+1.*lc_sq, 0. , lc_sq};
+ Point(16) = {-a_lat/2. , d_sq/2+1.*lc_sq, 0. , lc_sq};
+ Point(17) = {-a_lat/2. ,-d_sq/2+1.*lc_sq, 0. , lc_sq};
+ Point(18) = { a_lat/2. ,-d_sq/2+1.*lc_sq, 0. , lc_sq};
+ Point(19) = { a_lat/2. , d_sq/2-1.*lc_sq, 0. , lc_sq};
+ Point(20) = {-a_lat/2. , d_sq/2-1.*lc_sq, 0. , lc_sq};
+ // Point(21) = { 0. , d_sq/2-d_sq/7, 0. , lc_sq};
+ // Point(22) = { 0. ,-d_sq/2+d_sq/7, 0. , lc_sq};
+
+ Line(1) = {9, 1};
+ Line(2) = {1, 13};
+ Line(3) = {13, 5};
+ Line(4) = {5, 17};
+ Line(5) = {17, 20};
+ Line(6) = {20, 8};
+ Line(7) = {8, 16};
+ Line(8) = {16, 4};
+ Line(9) = {4, 12};
+ Line(10) = {10, 2};
+ Line(11) = {2, 14};
+ Line(12) = {14, 6};
+ Line(13) = {6, 18};
+ Line(14) = {18, 19};
+ Line(15) = {19, 7};
+ Line(16) = {7, 15};
+ Line(17) = {15, 3};
+ Line(18) = {3, 11};
+ Line(19) = {9, 10};
+ Line(20) = {1, 2};
+ Line(21) = {5, 6};
+ Line(22) = {8, 7};
+ Line(23) = {4, 3};
+ Line(24) = {12, 11};
+
+ Line(25) = {13, 14};
+ Line(26) = {17, 18};
+ Line(27) = {20, 19};
+ Line(28) = {16, 15};
+ Curve Loop(1) = {1, 20, -10, -19};
+ Plane Surface(1) = {1};
+ Curve Loop(2) = {2, 25, -11, -20};
+ Plane Surface(2) = {2};
+ Curve Loop(3) = {3, 21, -12, -25};
+ Plane Surface(3) = {3};
+ Curve Loop(4) = {4, 26, -13, -21};
+ Plane Surface(4) = {4};
+ Curve Loop(5) = {5, 27, -14, -26};
+ Plane Surface(5) = {5};
+ Curve Loop(6) = {6, 22, -15, -27};
+ Plane Surface(6) = {6};
+ Curve Loop(7) = {7, 28, -16, -22};
+ Plane Surface(7) = {7};
+ Curve Loop(8) = {8, 23, -17, -28};
+ Plane Surface(8) = {8};
+ Curve Loop(9) = {9, 24, -18, -23};
+ Plane Surface(9) = {9};
+
+ Transfinite Surface { 3 } AlternateLeft;
+ Transfinite Surface { 4 } AlternateRight;
+
+ Transfinite Surface { 6 } AlternateLeft;
+ Transfinite Surface { 7 } AlternateRight;
+
+ Periodic Line { 10,11,12,13,14,15,16,17,18 } = { 1,2,3,4,5,6,7,8,9 } Translate {a_lat,0,0} ;
+ // Point{21,22} In Surface{3};
+
+ Physical Surface('IN',100) = {4,5,6}; // 1 dom in
+ Physical Surface('OUT',101) = {2,3,7,8}; // 2 dom out
+ Physical Surface('PMLTOP',102) = {9}; // PML top
+ Physical Surface('PMLBOT',103) = {1}; // PML bot
+ Physical Line('BXM',1001) = {1,2,3,4,5,6,7,8,9}; // bloch x left
+ Physical Line('BXP',1002) = {10,11,12,13,14,15,16,17,18}; // bloch x right
+ Physical Line('LTOP',1003) = {24}; // top bound
+ Physical Line('LBOT',1004) = {19}; // bot bound
+ Physical Line('LAG',1005) = {21,22}; // bound for lag
+ Physical Point('PRINT',10000) = {1}; // Printpoint
+EndIf
+
+If (flag_geom==5)
+ epsc = lc_sq*0.5;
+ lc_sqa = lc_sq;
+ 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};
+ Point(4) = {-a_lat/2. , d_sq/2.+space2pml, 0. , lc_cell};
+ Point(5) = {-d_sq/2. ,-d_sq/2., 0. , lc_sq*15};
+ Point(6) = { d_sq/2. ,-d_sq/2., 0. , lc_sq*15};
+ Point(7) = { d_sq/2. , d_sq/2., 0. , lc_sq*15};
+ Point(8) = {-d_sq/2. , d_sq/2., 0. , lc_sq*15};
+ 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) = {-d_sq/2.-epsc , d_sq/2.+epsc, 0. , lc_sqa};
Point(14) = {-d_sq/2. , d_sq/2.+epsc, 0. , lc_sqa};
Point(15) = {-d_sq/2.+epsc , d_sq/2.+epsc, 0. , lc_sqa};
@@ -165,12 +434,10 @@ If (flag_Tmesh==1)
Point(30) = {-d_sq/2.-epsc ,-d_sq/2.+epsc, 0. , lc_sqa};
Point(31) = {-d_sq/2.-epsc , d_sq/2.-epsc, 0. , lc_sqa};
Point(32) = {-d_sq/2.-epsc , d_sq/2. , 0. , lc_sqa};
-
Point(33) = {-d_sq/2.+epsc ,-d_sq/2.+epsc, 0. , lc_sqa};
Point(34) = { d_sq/2.-epsc ,-d_sq/2.+epsc, 0. , lc_sqa};
Point(35) = { d_sq/2.-epsc , d_sq/2.-epsc, 0. , lc_sqa};
Point(36) = {-d_sq/2.+epsc , d_sq/2.-epsc, 0. , lc_sqa};
-
Point(37) = {-d_sq/2.+epsc ,-d_sq/2., 0. , lc_sqa};
Point(38) = { d_sq/2.-epsc ,-d_sq/2., 0. , lc_sqa};
Point(39) = { d_sq/2.-epsc , d_sq/2., 0. , lc_sqa};
@@ -179,12 +446,10 @@ If (flag_Tmesh==1)
Point(47) = { d_sq/2. ,-d_sq/2.+epsc, 0. , lc_sqa};
Point(48) = { d_sq/2. , d_sq/2.-epsc, 0. , lc_sqa};
Point(49) = {-d_sq/2. , d_sq/2.-epsc, 0. , lc_sqa};
-
Point(54) = {-a_lat/2. , a_lat/2., 0. , lc_sqa};
Point(55) = {-a_lat/2. ,-a_lat/2., 0. , lc_sqa};
Point(56) = { a_lat/2. , a_lat/2., 0. , lc_sqa};
Point(58) = { a_lat/2. ,-a_lat/2., 0. , lc_sqa};
-
Line(1) = {1, 2};
Line(3) = {3, 4};
Line(9) = {4, 12};
@@ -193,7 +458,6 @@ If (flag_Tmesh==1)
Line(12) = {1, 9};
Line(13) = {9, 10};
Line(14) = {10, 2};
-
Line(35) = {8, 45};
Line(36) = {45, 45};
Line(37) = {39, 39};
@@ -279,6 +543,7 @@ If (flag_Tmesh==1)
Line(117) = {58, 56};
Line(118) = {56, 3};
+
Line Loop(1) = {53, -55, -112};
Plane Surface(1) = {1};
Line Loop(2) = {55, 54, -93};
@@ -368,48 +633,215 @@ If (flag_Tmesh==1)
Plane Surface(45) = {-45,-46};
Line Loop(47) = {89, 90, 91, 92};
Plane Surface(46) = {-47};
-
-
Periodic Line { 14,116,117,118,11 } = {12,113, 114, 115,9 } Translate {a_lat,0,0} ;
-
Periodic Line {47} = {110} Translate {epsc,0,0} ;
Periodic Line {89} = {47} Translate {epsc,0,0} ;
Periodic Line {41} = {91} Translate {epsc,0,0} ;
Periodic Line {100} = {41} Translate {epsc,0,0} ;
-
Periodic Line {44} = {105} Translate {0,epsc,0} ;
Periodic Line {92} = {44} Translate {0,epsc,0} ;
Periodic Line {38} = {90} Translate {0,epsc,0} ;
Periodic Line {95} = {38} Translate {0,epsc,0} ;
+ Transfinite Surface { 9} Alternate;
+ Transfinite Surface {10} AlternateLeft;
+ Transfinite Surface {19} Alternate;
+ Transfinite Surface {20} AlternateLeft;
+ Transfinite Surface {29} Alternate;
+ Transfinite Surface {30} AlternateLeft;
+ Transfinite Surface {39} Alternate;
+ Transfinite Surface {40} AlternateLeft;
- // Transfinite Surface { expression-list } | "*" < = { expression-list } > < Left | Right | Alternate | AlternateRight | AlternateLeft > ;
- Transfinite Surface { 9 } AlternateLeft;
- Transfinite Surface { 10 } AlternateRight;
+ Physical Surface('IN',100) = {7,8,9,17,18,19,23,24,29,37,38,40,46}; // 1 dom in
+ Physical Surface('OUT',101) = {45, 1, 2, 3, 4, 10, 5, 6, 39, 36, 35, 34, 33, 32, 31, 30, 25, 26, 27, 28, 22, 21, 20, 16, 15, 14, 13, 12, 11}; // out
+ Physical Surface('PMLTOP',102) = {43}; // PML top
+ Physical Surface('PMLBOT',103) = {44}; // PML bot
+ Physical Line('BXM',1001) = {12,113,114,115,9}; // bloch x left
+ Physical Line('BXP',1002) = {14,116,117,118,11}; // bloch x right
+ Physical Line('TOP',1003) = {10}; // top bound
+ Physical Line('BOT',1004) = {13}; // bot bound
+ Physical Line('LAG',1005) = {38,39,40,41,42,43,44,45,46,47,48,35}; // bound for lag
+ Physical Point('PRINT',10000) = {1}; // Printpoint
+EndIf
- Transfinite Surface { 19 } AlternateLeft;
- Transfinite Surface { 20 } AlternateLeft;
+If (flag_geom==6)
+ epsc = lc_sq/2;
+ lc_sqa = lc_sq;
+ s2 = 1-Sqrt(2)/2;
+ c8 = epsc*Cos(Pi/8);
+ s8 = epsc*Sin(Pi/8);
+ 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};
+ Point(4) = {-a_lat/2. , d_sq/2.+space2pml, 0. , lc_cell};
+ Point(5) = {-d_sq/2.+s2*epsc,-d_sq/2.+s2*epsc, 0. , lc_sq*15};
+ Point(6) = { d_sq/2.-s2*epsc,-d_sq/2.+s2*epsc, 0. , lc_sq*15};
+ Point(7) = { d_sq/2.-s2*epsc, d_sq/2.-s2*epsc, 0. , lc_sq*15};
+ Point(8) = {-d_sq/2.+s2*epsc, d_sq/2.-s2*epsc, 0. , lc_sq*15};
+ 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(32) = {-d_sq/2.+epsc-c8*2 , d_sq/2.-epsc+s8*2, 0. , lc_sqa};
+ Point(14) = {-d_sq/2.+epsc-s8*2 , d_sq/2.-epsc+c8*2, 0. , lc_sqa};
+ Point(15) = {-d_sq/2.+epsc , d_sq/2.+epsc, 0. , lc_sqa};
+ Point(16) = { d_sq/2.-epsc , d_sq/2.+epsc, 0. , lc_sqa};
+ Point(17) = { d_sq/2.-epsc+s8*2 , d_sq/2.-epsc+c8*2, 0. , lc_sqa};
+ Point(19) = { d_sq/2.-epsc+c8*2 , d_sq/2.-epsc+s8*2, 0. , lc_sqa};
+ Point(20) = { d_sq/2.+epsc , d_sq/2.-epsc, 0. , lc_sqa};
+ Point(21) = { d_sq/2.+epsc ,-d_sq/2.+epsc, 0. , lc_sqa};
+ Point(13) = {-d_sq/2.-epsc+2*s2*epsc , d_sq/2.+epsc-2*s2*epsc, 0. , lc_sqa};
+ Point(23) = { d_sq/2.+epsc-2*s2*epsc ,-d_sq/2.-epsc+2*s2*epsc, 0. , lc_sqa};
+ Point(28) = {-d_sq/2.-epsc+2*s2*epsc ,-d_sq/2.-epsc+2*s2*epsc, 0. , lc_sqa};
+ Point(18) = { d_sq/2.+epsc-2*s2*epsc , d_sq/2.+epsc-2*s2*epsc, 0. , lc_sqa};
+ Point(24) = { d_sq/2.-epsc+s8*2 ,-d_sq/2.+epsc-c8*2, 0. , lc_sqa};
+ Point(22) = { d_sq/2.-epsc+c8*2 ,-d_sq/2.+epsc-s8*2, 0. , lc_sqa};
+ Point(25) = { d_sq/2.-epsc ,-d_sq/2.-epsc, 0. , lc_sqa};
+ Point(26) = {-d_sq/2.+epsc ,-d_sq/2.-epsc, 0. , lc_sqa};
+ Point(27) = {-d_sq/2.+epsc-s8*2 ,-d_sq/2.+epsc-c8*2, 0. , lc_sqa};
+ Point(29) = {-d_sq/2.+epsc-c8*2 ,-d_sq/2.+epsc-s8*2 , 0. , lc_sqa};
- Transfinite Surface { 29 } AlternateLeft;
- Transfinite Surface { 30 } AlternateLeft;
- Transfinite Surface { 39 } AlternateLeft;
- Transfinite Surface { 40 } AlternateRight;
+ Point(30) = {-d_sq/2.-epsc ,-d_sq/2.+epsc, 0. , lc_sqa};
+ Point(31) = {-d_sq/2.-epsc , d_sq/2.-epsc, 0. , lc_sqa};
+ Point(33) = {-d_sq/2.+epsc ,-d_sq/2.+epsc, 0. , lc_sqa};
+ Point(34) = { d_sq/2.-epsc ,-d_sq/2.+epsc, 0. , lc_sqa};
+ Point(35) = { d_sq/2.-epsc , d_sq/2.-epsc, 0. , lc_sqa};
+ Point(36) = {-d_sq/2.+epsc , d_sq/2.-epsc, 0. , lc_sqa};
+ Point(37) = {-d_sq/2.+epsc ,-d_sq/2., 0. , lc_sqa};
+ Point(38) = { d_sq/2.-epsc ,-d_sq/2., 0. , lc_sqa};
+ Point(39) = { d_sq/2.-epsc , d_sq/2., 0. , lc_sqa};
+ Point(45) = {-d_sq/2.+epsc , d_sq/2., 0. , lc_sqa};
+ Point(46) = {-d_sq/2. ,-d_sq/2.+epsc, 0. , lc_sqa};
+ Point(47) = { d_sq/2. ,-d_sq/2.+epsc, 0. , lc_sqa};
+ Point(48) = { d_sq/2. , d_sq/2.-epsc, 0. , lc_sqa};
+ Point(49) = {-d_sq/2. , d_sq/2.-epsc, 0. , lc_sqa};
+ Point(54) = {-a_lat/2. , a_lat/2., 0. , lc_sqa};
+ Point(55) = {-a_lat/2. ,-a_lat/2., 0. , lc_sqa};
+ Point(56) = { a_lat/2. , a_lat/2., 0. , lc_sqa};
+ Point(58) = { a_lat/2. ,-a_lat/2., 0. , lc_sqa};
+ Line(1) = {1, 2};
+ Line(3) = {3, 4};
+ 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(35) = {8, 45};
+ Circle(35) = {8, 36,45};
- Physical Surface(100) = {7,8,9,17,18,19,23,24,29,37,38,40,46}; // 1 dom in
- Physical Surface(101) = {45, 1, 2, 3, 4, 10, 5, 6, 39, 36, 35, 34, 33, 32, 31, 30, 25, 26, 27, 28, 22, 21, 20, 16, 15, 14, 13, 12, 11}; // out
- Physical Surface(102) = {43}; // PML top
- Physical Surface(103) = {44}; // PML bot
- Physical Line(1001) = {12,113,114,115,9}; // bloch x left
- Physical Line(1002) = {14,116,117,118,11}; // bloch x right
- Physical Line(1003) = {10}; // top bound
- Physical Line(1004) = {13}; // bot bound
- Physical Line(1005) = {38,39,40,41,42,43,44,45,46,47,48,35}; // bound for lag
- Physical Point(10000) = {1}; // Printpoint
-EndIf
+ Line(38) = {45, 39};
+ Circle(39) = {39, 35, 7};
+ Circle(40) = {7, 35, 48};
+
+ Line(41) = {48, 47};
+ Circle(42) = {47, 34, 6};
+ Circle(43) = {6, 34, 38};
+
+ Line(44) = {38, 37};
+ Circle(46) = {5, 33, 46};
+ Circle(45) = {37, 33, 5};
+
+ Line(47) = {46, 49};
+ Circle(48) = {49, 36, 8};
+
+ Line(49) = {31, 49};
+ Line(50) = {49, 36};
+ Line(51) = {36, 45};
+ Line(52) = {45, 15};
+ Line(59) = {16, 39};
+ Line(60) = {39, 35};
+ Line(61) = {35, 48};
+ Line(62) = {48, 20};
+ Line(69) = {34, 47};
+ Line(70) = {47, 21};
+ Line(71) = {34, 38};
+ Line(72) = {38, 25};
+ Line(79) = {30, 46};
+ Line(80) = {46, 33};
+ Line(81) = {33, 37};
+ Line(82) = {37, 26};
+ Line(89) = {33, 36};
+ Line(90) = {36, 35};
+ Line(91) = {35, 34};
+ Line(92) = {34, 33};
+ Line(95) = {15, 16};
+ Line(100) = {20, 21};
+ Line(105) = {25, 26};
+ Line(110) = {30, 31};
+ Line(113) = {1, 55};
+ Line(114) = {55, 54};
+ Line(115) = {54, 4};
+ Line(116) = {2, 58};
+ Line(117) = {58, 56};
+ Line(118) = {56, 3};
-If (flag_o2g==1)
- Mesh.ElementOrder = 2;
-Else
- Mesh.ElementOrder = 1;
+ Line Loop(9) = {51, 38, 60, -90};
+ Plane Surface(9) = {-9};
+ Line Loop(10) = {38, -59, -95, -52};
+ Plane Surface(10) = {10};
+ Line Loop(19) = {91, 69, -41, -61};
+ Plane Surface(19) = {19};
+ Line Loop(20) = {70, -100, -62, 41};
+ Plane Surface(20) = {20};
+ Line Loop(29) = {92, 81, -44, -71};
+ Plane Surface(29) = {29};
+ Line Loop(30) = {82, -105, -72, 44};
+ Plane Surface(30) = {30};
+ Line Loop(39) = {110, 49, -47, -79};
+ Plane Surface(39) = {-39};
+ Line Loop(40) = {50, -89, -80, 47};
+ Plane Surface(40) = {-40};
+ Line Loop(43) = {3, 9, 10, 11};
+ Plane Surface(43) = {-43};
+ Line Loop(44) = {1, -14, -13, -12};
+ Plane Surface(44) = {-44};
+ Line Loop(45) = {113, 114, 115, -3, -118, -117, -116, -1};
+ Curve Loop(48) = {48, 35, -51, -50};
+ Plane Surface(47) = {48};
+ Curve Loop(49) = {60, 61, -40, -39};
+ Plane Surface(48) = {49};
+ Curve Loop(50) = {71, -43, -42, -69};
+ Plane Surface(49) = {50};
+ Curve Loop(51) = {80, 81, 45, 46};
+ Plane Surface(50) = {51};
+
+ Curve Loop(52) = {95, 59, 39, 40, 62, 100, -70, 42, 43, 72, 105, -82, 45, 46, -79, 110, 49, 48, 35, 52};
+ Plane Surface(51) = {45, 52};
+ Line Loop(47) = {89, 90, 91, 92};
+ Plane Surface(46) = {-47};
+ Periodic Line { 14,116,117,118,11 } = {12,113, 114, 115,9 } Translate {a_lat,0,0} ;
+ Periodic Line {47} = {110} Translate {epsc,0,0} ;
+ Periodic Line {89} = {47} Translate {epsc,0,0} ;
+ Periodic Line {41} = {91} Translate {epsc,0,0} ;
+ Periodic Line {100} = {41} Translate {epsc,0,0} ;
+ Periodic Line {44} = {105} Translate {0,epsc,0} ;
+ Periodic Line {92} = {44} Translate {0,epsc,0} ;
+ Periodic Line {38} = {90} Translate {0,epsc,0} ;
+ Periodic Line {95} = {38} Translate {0,epsc,0} ;
+ Transfinite Surface { 9} AlternateLeft;
+ Transfinite Surface {10} Alternate;
+ Transfinite Surface {19} Alternate;
+ Transfinite Surface {20} Alternate;
+ Transfinite Surface {29} Alternate;
+ Transfinite Surface {30} Alternate;
+ Transfinite Surface {39} Alternate;
+ Transfinite Surface {40} AlternateLeft;
+
+ Point {13,14,32 , 17,18,19 , 22,23,24, 27,28,29 } In Surface{51};
+
+ Physical Surface('IN',100) = {7,8,9,17,18,19,23,24,29,37,38,40,46,47,48,49,50}; // 1 dom in
+ Physical Surface('OUT',101) = {51,45, 1, 2, 3, 4, 10, 5, 6, 39, 36, 35, 34, 33, 32, 31, 30, 25, 26, 27, 28, 22, 21, 20, 16, 15, 14, 13, 12, 11}; // out
+ Physical Surface('PMLTOP',102) = {43}; // PML top
+ Physical Surface('PMLBOT',103) = {44}; // PML bot
+ Physical Line('BXM',1001) = {12,113,114,115,9}; // bloch x left
+ Physical Line('BXP',1002) = {14,116,117,118,11}; // bloch x right
+ Physical Line('TOP',1003) = {10}; // top bound
+ Physical Line('BOT',1004) = {13}; // bot bound
+ Physical Line('LAG',1005) = {38,39,40,41,42,43,44,45,46,47,48,35}; // bound for lag
+ Physical Point('PRINT',10000) = {1}; // Printpoint
EndIf
+
+Mesh.ElementOrder = flag_og;
+
diff --git a/NonLinearEVP/NonLinearEVP.pro b/NonLinearEVP/NonLinearEVP.pro
index a27671f8dd7941f2ec6dffbe617b76a2bd0f801a..5dd74f50c0f54c504b897ab111a4a9b94762fc70 100644
--- a/NonLinearEVP/NonLinearEVP.pro
+++ b/NonLinearEVP/NonLinearEVP.pro
@@ -157,7 +157,7 @@ FunctionSpace {
// 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)
+ If(flag_og==1)
// 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]; }
@@ -169,7 +169,7 @@ FunctionSpace {
}
Constraint {
{ NameOfCoef un; EntityType NodesOf ; NameOfConstraint BlochX; }
- If(flag_o2g==0)
+ If(flag_og==1)
{ NameOfCoef un2e; EntityType EdgesOf ; NameOfConstraint BlochX; }
EndIf
If(flag_o2i==1)
@@ -181,7 +181,7 @@ FunctionSpace {
{ 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
+ If(flag_og==1) // 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)
@@ -191,7 +191,7 @@ FunctionSpace {
}
Constraint {
{ NameOfCoef un; EntityType NodesOf ; NameOfConstraint BlochX; }
- If(flag_o2g==0)
+ If(flag_og==1)
{ NameOfCoef un2e; EntityType EdgesOf ; NameOfConstraint BlochX; }
EndIf
If(flag_o2i==1)
diff --git a/NonLinearEVP/NonLinearEVP_data.geo b/NonLinearEVP/NonLinearEVP_data.geo
index abc94c3641f69cac413ab77091bfb799106296f9..a8a74de6c4f1eaba3d3d9b361f9615c26a59c228 100644
--- a/NonLinearEVP/NonLinearEVP_data.geo
+++ b/NonLinearEVP/NonLinearEVP_data.geo
@@ -34,6 +34,7 @@ DefineConstant[ a_lat = {50 , Name StrCat[pp0 , "1grating period d [nm]"]
DefineConstant[
+ flag_shape = {1 , Name StrCat[pp0 , "1Select the geometry"] , Choices {1="square",2="slab",3="circle"} },
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.} ,
@@ -54,11 +55,11 @@ DefineConstant[
eig_max_im = {4. , Name StrCat[pp3 , "6EV imag max [2Pic\a]"] , Highlight Str[colorpp3] , Closed close_menu} ,
paramaille = {5 , Name StrCat[pp4 , "0number of mesh elements per period []"] , Highlight Str[colorpp4], Min 2, Max 10} ,
- flag_Tmesh = {1 , Name StrCat[pp4 , "2locally structured mesh?"] , Choices {0="unstruct",1="struct"} },
- flag_o2g = {1 , Name StrCat[pp4 , "3Geometrical order"] , Choices {0="order 1 (linear)",1="order 2 (curved)"} },
+ flag_tm = {1 , Name StrCat[pp4 , "2structured symmetric mesh"] , Choices {0="no",1="yes"} },
+ flag_og = {1 , Name StrCat[pp4 , "3Geometrical order"] , Choices {1="order 1 (linear)",2="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_round = {0 , Name StrCat[pp4 , "5Corner rounding/0Enable rounding of corners!"], Choices{0,1}, ServerAction "ResetDatabase"},
+ corner_rad_frac= {0.04 , Name StrCat[pp4 , "5Corner rounding/1corner radius (fraction of period)"] , 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"],
@@ -66,6 +67,37 @@ DefineConstant[
Choices {0="Aux_E" ,1="PEP_E" ,2="NEP_E" ,3="Lag_E" ,4="PEP_h"},ServerAction "ResetDatabase"}
];
+
+If (flag_shape==1)
+ If (flag_round==0)
+ If (flag_tm==0)
+ flag_geom=1;
+ Else
+ flag_geom=5;
+ EndIf
+ Else
+ If (flag_tm==0)
+ flag_geom=1;
+ Else
+ flag_geom=6;
+ EndIf
+ EndIf
+EndIf
+
+If (flag_shape==2)
+ If (flag_tm==0)
+ flag_geom=2;
+ Else
+ flag_geom=4;
+ EndIf
+EndIf
+
+If (flag_shape==3)
+ flag_geom=3;
+EndIf
+
+
+
// Normalized units so that 2*pi*c/a=1
cel = a_lat/(2*Pi);
epsf = 8.854187817e-3;
@@ -88,7 +120,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;
+corner_rad = a_lat*corner_rad_frac;
eps_oo_2 = 1;
om_d_2 = 0;