diff --git a/ElectromagneticScattering/README.md b/ElectromagneticScattering/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..248de22af1b92506c3adca3ae9e731d737b7ac24
--- /dev/null
+++ b/ElectromagneticScattering/README.md
@@ -0,0 +1,23 @@
+A Onelab model for 3D scattering problems in nanophotonics.
+
+## Synopsis
+
+This project contains a [Onelab](http://onelab.info/wiki/ONELAB) model for solving various 3D electromagnetic scattering problems on an isolated object :
+* T-matrix computation: See https://arxiv.org/abs/1802.00596 for details
+* Quasi-normal modes
+* Plane wave response
+* Green's function (LDOS)
+
+## Installation
+
+This model requires the following (open-source) programs:
+* [onelab](http://onelab.info/wiki/ONELAB), which bundles both [gmsh](http://www.gmsh.info/) and [getdp](http://www.getdp.info/)
+* python (v2.7.x or v3.6.x) with numpy, scipy and matplotlib
+
+## Running the model
+
+Open `scattering.pro` with Gmsh.
+
+## Authors
+
+Guillaume Demésy and Brian Stout
\ No newline at end of file
diff --git a/ElectromagneticScattering/scattering.geo b/ElectromagneticScattering/scattering.geo
new file mode 100644
index 0000000000000000000000000000000000000000..b51d28eb4bfe11e95a8a4c9d1f13cc6d1fbd812c
--- /dev/null
+++ b/ElectromagneticScattering/scattering.geo
@@ -0,0 +1,753 @@
+///////////////////////////////
+// Author : Guillaume Demesy //
+// scattering.geo //
+///////////////////////////////
+
+Include "scattering_data.geo";
+In_n = Sqrt[Fabs[epsr_In_re]];
+// Out_n = Sqrt[Fabs[epsr_Out_re]];
+
+paramaille_pml = paramaille/1.1;
+
+Out_lc = lambda_bg/paramaille;
+PML_lc = lambda_bg/paramaille_pml;
+In_lc = lambda0/(paramaille*In_n*refine_scat);
+CenterScat_lc = lambda0/(paramaille*In_n);
+
+If(flag_cartpml==1)
+ dom_x = r_pml_in*2;
+ dom_y = r_pml_in*2;
+ dom_z = r_pml_in*2;
+ PML_bot = pml_size;
+ PML_lat = pml_size;
+ PML_top = pml_size;
+ Point(1) = {-dom_x/2,-dom_y/2, -dom_z/2, Out_lc};
+ Point(2) = {-dom_x/2, dom_y/2, -dom_z/2, Out_lc};
+ Point(3) = { dom_x/2, dom_y/2, -dom_z/2, Out_lc};
+ Point(4) = { dom_x/2,-dom_y/2, -dom_z/2, Out_lc};
+ Point(5) = {-dom_x/2,-dom_y/2, dom_z/2, Out_lc};
+ Point(6) = {-dom_x/2, dom_y/2, dom_z/2, Out_lc};
+ Point(7) = { dom_x/2, dom_y/2, dom_z/2, Out_lc};
+ Point(8) = { dom_x/2,-dom_y/2, dom_z/2, Out_lc};
+ Point(9) = {-dom_x/2,-dom_y/2, -dom_z/2-PML_bot, PML_lc};
+ Point(10) = {-dom_x/2, dom_y/2, -dom_z/2-PML_bot, PML_lc};
+ Point(11) = { dom_x/2, dom_y/2, -dom_z/2-PML_bot, PML_lc};
+ Point(12) = { dom_x/2,-dom_y/2, -dom_z/2-PML_bot, PML_lc};
+ Point(13) = {-dom_x/2,-dom_y/2, dom_z/2+PML_top, PML_lc};
+ Point(14) = {-dom_x/2, dom_y/2, dom_z/2+PML_top, PML_lc};
+ Point(15) = { dom_x/2, dom_y/2, dom_z/2+PML_top, PML_lc};
+ Point(16) = { dom_x/2,-dom_y/2, dom_z/2+PML_top, PML_lc};
+ Point(17) = {-dom_x/2,-dom_y/2-PML_lat, -dom_z/2 , PML_lc};
+ Point(18) = {-dom_x/2, dom_y/2+PML_lat, -dom_z/2 , PML_lc};
+ Point(19) = { dom_x/2, dom_y/2+PML_lat, -dom_z/2 , PML_lc};
+ Point(20) = { dom_x/2,-dom_y/2-PML_lat, -dom_z/2 , PML_lc};
+ Point(21) = {-dom_x/2,-dom_y/2-PML_lat, dom_z/2 , PML_lc};
+ Point(22) = {-dom_x/2, dom_y/2+PML_lat, dom_z/2 , PML_lc};
+ Point(23) = { dom_x/2, dom_y/2+PML_lat, dom_z/2 , PML_lc};
+ Point(24) = { dom_x/2,-dom_y/2-PML_lat, dom_z/2 , PML_lc};
+ Point(25) = {-dom_x/2-PML_lat,-dom_y/2, -dom_z/2 , PML_lc};
+ Point(26) = {-dom_x/2-PML_lat, dom_y/2, -dom_z/2 , PML_lc};
+ Point(27) = { dom_x/2+PML_lat, dom_y/2, -dom_z/2 , PML_lc};
+ Point(28) = { dom_x/2+PML_lat,-dom_y/2, -dom_z/2 , PML_lc};
+ Point(29) = {-dom_x/2-PML_lat,-dom_y/2, dom_z/2 , PML_lc};
+ Point(30) = {-dom_x/2-PML_lat, dom_y/2, dom_z/2 , PML_lc};
+ Point(31) = { dom_x/2+PML_lat, dom_y/2, dom_z/2 , PML_lc};
+ Point(32) = { dom_x/2+PML_lat,-dom_y/2, dom_z/2 , PML_lc};
+ Point(33) = {-dom_x/2-PML_lat,-dom_y/2-PML_lat, -dom_z/2-PML_bot, PML_lc};
+ Point(34) = {-dom_x/2-PML_lat, dom_y/2+PML_lat, -dom_z/2-PML_bot, PML_lc};
+ Point(35) = { dom_x/2+PML_lat, dom_y/2+PML_lat, -dom_z/2-PML_bot, PML_lc};
+ Point(36) = { dom_x/2+PML_lat,-dom_y/2-PML_lat, -dom_z/2-PML_bot, PML_lc};
+ Point(37) = {-dom_x/2-PML_lat,-dom_y/2-PML_lat, dom_z/2+PML_top, PML_lc};
+ Point(38) = {-dom_x/2-PML_lat, dom_y/2+PML_lat, dom_z/2+PML_top, PML_lc};
+ Point(39) = { dom_x/2+PML_lat, dom_y/2+PML_lat, dom_z/2+PML_top, PML_lc};
+ Point(40) = { dom_x/2+PML_lat,-dom_y/2-PML_lat, dom_z/2+PML_top, PML_lc};
+ Point(41) = {-dom_x/2,-dom_y/2-PML_lat, -dom_z/2-PML_bot, PML_lc};
+ Point(42) = {-dom_x/2, dom_y/2+PML_lat, -dom_z/2-PML_bot, PML_lc};
+ Point(43) = { dom_x/2,-dom_y/2-PML_lat, -dom_z/2-PML_top, PML_lc};
+ Point(44) = { dom_x/2, dom_y/2+PML_lat, -dom_z/2-PML_top, PML_lc};
+ Point(45) = {-dom_x/2,-dom_y/2-PML_lat, dom_z/2+PML_bot, PML_lc};
+ Point(46) = {-dom_x/2, dom_y/2+PML_lat, dom_z/2+PML_bot, PML_lc};
+ Point(47) = { dom_x/2,-dom_y/2-PML_lat, dom_z/2+PML_top, PML_lc};
+ Point(48) = { dom_x/2, dom_y/2+PML_lat, dom_z/2+PML_top, PML_lc};
+ Point(49) = {-dom_x/2-PML_lat, -dom_y/2, -dom_z/2-PML_bot, PML_lc};
+ Point(50) = { dom_x/2+PML_lat, -dom_y/2, -dom_z/2-PML_bot, PML_lc};
+ Point(51) = {-dom_x/2-PML_lat, -dom_y/2, dom_z/2+PML_bot, PML_lc};
+ Point(52) = { dom_x/2+PML_lat, -dom_y/2, dom_z/2+PML_bot, PML_lc};
+ Point(53) = {-dom_x/2-PML_lat, dom_y/2, -dom_z/2-PML_bot, PML_lc};
+ Point(54) = { dom_x/2+PML_lat, dom_y/2, -dom_z/2-PML_bot, PML_lc};
+ Point(55) = {-dom_x/2-PML_lat, dom_y/2, dom_z/2+PML_bot, PML_lc};
+ Point(56) = { dom_x/2+PML_lat, dom_y/2, dom_z/2+PML_bot, PML_lc};
+ Point(57) = {-dom_x/2-PML_lat, -dom_y/2-PML_lat, -dom_z/2, PML_lc};
+ Point(58) = { dom_x/2+PML_lat, -dom_y/2-PML_lat, -dom_z/2, PML_lc};
+ Point(59) = {-dom_x/2-PML_lat, dom_y/2+PML_lat, -dom_z/2, PML_lc};
+ Point(60) = { dom_x/2+PML_lat, dom_y/2+PML_lat, -dom_z/2, PML_lc};
+ Point(61) = {-dom_x/2-PML_lat, -dom_y/2-PML_lat, dom_z/2, PML_lc};
+ Point(62) = { dom_x/2+PML_lat, -dom_y/2-PML_lat, dom_z/2, PML_lc};
+ Point(63) = {-dom_x/2-PML_lat, dom_y/2+PML_lat, dom_z/2, PML_lc};
+ Point(64) = { dom_x/2+PML_lat, dom_y/2+PML_lat, dom_z/2, PML_lc};
+ Point(73) = { 0,0,0, CenterScat_lc};
+
+ Line(1) = {33, 49};
+ Line(2) = {49, 53};
+ Line(3) = {53, 34};
+ Line(4) = {34, 42};
+ Line(5) = {42, 44};
+ Line(6) = {44, 35};
+ Line(7) = {35, 54};
+ Line(8) = {54, 11};
+ Line(9) = {11, 10};
+ Line(10) = {10, 53};
+ Line(11) = {49, 9};
+ Line(12) = {9, 12};
+ Line(13) = {12, 50};
+ Line(14) = {50, 54};
+ Line(15) = {50, 36};
+ Line(16) = {36, 43};
+ Line(17) = {43, 41};
+ Line(18) = {41, 33};
+ Line(19) = {41, 9};
+ Line(20) = {9, 10};
+ Line(21) = {10, 42};
+ Line(22) = {43, 12};
+ Line(23) = {12, 11};
+ Line(24) = {11, 44};
+ Line(25) = {57, 25};
+ Line(26) = {25, 26};
+ Line(27) = {26, 59};
+ Line(28) = {59, 18};
+ Line(29) = {18, 19};
+ Line(30) = {19, 60};
+ Line(31) = {60, 27};
+ Line(32) = {27, 28};
+ Line(33) = {28, 58};
+ Line(34) = {58, 20};
+ Line(35) = {20, 17};
+ Line(36) = {17, 57};
+ Line(37) = {25, 1};
+ Line(38) = {1, 4};
+ Line(39) = {4, 28};
+ Line(40) = {27, 3};
+ Line(41) = {3, 2};
+ Line(42) = {2, 26};
+ Line(43) = {17, 1};
+ Line(44) = {1, 2};
+ Line(45) = {2, 18};
+ Line(46) = {19, 3};
+ Line(47) = {3, 4};
+ Line(48) = {4, 20};
+ Line(49) = {61, 21};
+ Line(50) = {21, 24};
+ Line(51) = {24, 62};
+ Line(52) = {62, 32};
+ Line(53) = {32, 31};
+ Line(54) = {31, 64};
+ Line(55) = {64, 23};
+ Line(56) = {23, 22};
+ Line(57) = {22, 63};
+ Line(58) = {63, 30};
+ Line(59) = {30, 29};
+ Line(60) = {29, 61};
+ Line(61) = {21, 5};
+ Line(62) = {5, 6};
+ Line(63) = {6, 22};
+ Line(64) = {23, 7};
+ Line(65) = {7, 8};
+ Line(66) = {8, 24};
+ Line(67) = {32, 8};
+ Line(68) = {8, 5};
+ Line(69) = {5, 29};
+ Line(70) = {30, 6};
+ Line(71) = {6, 7};
+ Line(72) = {7, 31};
+ Line(73) = {37, 45};
+ Line(74) = {45, 47};
+ Line(75) = {47, 40};
+ Line(76) = {40, 52};
+ Line(77) = {52, 56};
+ Line(78) = {56, 39};
+ Line(79) = {39, 48};
+ Line(80) = {48, 46};
+ Line(81) = {46, 38};
+ Line(82) = {38, 55};
+ Line(83) = {55, 51};
+ Line(84) = {51, 37};
+ Line(85) = {45, 13};
+ Line(86) = {13, 14};
+ Line(87) = {14, 46};
+ Line(88) = {55, 14};
+ Line(89) = {14, 15};
+ Line(90) = {15, 56};
+ Line(91) = {52, 16};
+ Line(92) = {16, 13};
+ Line(93) = {13, 51};
+ Line(94) = {48, 15};
+ Line(95) = {15, 16};
+ Line(96) = {16, 47};
+ Line(97) = {37, 61};
+ Line(98) = {61, 57};
+ Line(99) = {57, 33};
+ Line(100) = {45, 21};
+ Line(101) = {21, 17};
+ Line(102) = {17, 41};
+ Line(103) = {47, 24};
+ Line(104) = {24, 20};
+ Line(105) = {20, 43};
+ Line(106) = {40, 62};
+ Line(107) = {62, 58};
+ Line(108) = {58, 36};
+ Line(109) = {52, 32};
+ Line(110) = {32, 28};
+ Line(111) = {28, 50};
+ Line(112) = {56, 31};
+ Line(113) = {31, 27};
+ Line(114) = {27, 54};
+ Line(115) = {39, 64};
+ Line(116) = {64, 60};
+ Line(117) = {60, 35};
+ Line(118) = {48, 23};
+ Line(119) = {23, 19};
+ Line(120) = {19, 44};
+ Line(121) = {46, 22};
+ Line(122) = {22, 18};
+ Line(123) = {18, 42};
+ Line(124) = {38, 63};
+ Line(125) = {63, 59};
+ Line(126) = {59, 34};
+ Line(127) = {55, 30};
+ Line(128) = {30, 26};
+ Line(129) = {26, 53};
+ Line(130) = {51, 29};
+ Line(131) = {29, 25};
+ Line(132) = {25, 49};
+ Line(133) = {14, 6};
+ Line(134) = {6, 2};
+ Line(135) = {2, 10};
+ Line(136) = {15, 7};
+ Line(137) = {7, 3};
+ Line(138) = {3, 11};
+ Line(139) = {16, 8};
+ Line(140) = {8, 4};
+ Line(141) = {4, 12};
+ Line(142) = {13, 5};
+ Line(143) = {5, 1};
+ Line(144) = {1, 9};
+
+ Line Loop(173) = {80, -87, 89, -94}; Plane Surface(174) = {173};
+ Line Loop(175) = {94, 90, 78, 79}; Plane Surface(176) = {175};
+ Line Loop(177) = {77, -90, 95, -91}; Plane Surface(178) = {177};
+ Line Loop(179) = {95, 92, 86, 89}; Plane Surface(180) = {179};
+ Line Loop(181) = {91, 96, 75, 76}; Plane Surface(182) = {181};
+ Line Loop(183) = {92, -85, 74, -96}; Plane Surface(184) = {183};
+ Line Loop(185) = {82, 88, 87, 81}; Plane Surface(186) = {185};
+ Line Loop(187) = {86, -88, 83, -93}; Plane Surface(188) = {187};
+ Line Loop(189) = {93, 84, 73, 85}; Plane Surface(190) = {189};
+ Line Loop(191) = {58, 70, 63, 57}; Plane Surface(192) = {191};
+ Line Loop(193) = {63, -56, 64, -71}; Plane Surface(194) = {193};
+ Line Loop(195) = {64, 72, 54, 55}; Plane Surface(196) = {195};
+ Line Loop(197) = {72, -53, 67, -65}; Plane Surface(198) = {197};
+ Line Loop(199) = {67, 66, 51, 52}; Plane Surface(200) = {199};
+ Line Loop(201) = {71, 65, 68, 62}; Plane Surface(202) = {201};
+ Line Loop(203) = {68, -61, 50, -66}; Plane Surface(204) = {203};
+ Line Loop(205) = {59, -69, 62, -70}; Plane Surface(206) = {205};
+ Line Loop(207) = {69, 60, 49, 61}; Plane Surface(208) = {207};
+ Line Loop(209) = {28, -45, 42, 27}; Plane Surface(210) = {209};
+ Line Loop(211) = {29, 46, 41, 45}; Plane Surface(212) = {211};
+ Line Loop(213) = {30, 31, 40, -46}; Plane Surface(214) = {213};
+ Line Loop(215) = {40, 47, 39, -32}; Plane Surface(216) = {215};
+ Line Loop(217) = {41, -44, 38, -47}; Plane Surface(218) = {217};
+ Line Loop(219) = {39, 33, 34, -48}; Plane Surface(220) = {219};
+ Line Loop(221) = {38, 48, 35, 43}; Plane Surface(222) = {221};
+ Line Loop(223) = {26, -42, -44, -37}; Plane Surface(224) = {223};
+ Line Loop(225) = {37, -43, 36, 25}; Plane Surface(226) = {225};
+ Line Loop(227) = {3, 4, -21, 10}; Plane Surface(228) = {227};
+ Line Loop(229) = {5, -24, 9, 21}; Plane Surface(230) = {229};
+ Line Loop(231) = {6, 7, 8, 24}; Plane Surface(232) = {231};
+ Line Loop(233) = {14, 8, -23, 13}; Plane Surface(234) = {233};
+ Line Loop(235) = {9, -20, 12, 23}; Plane Surface(236) = {235};
+ Line Loop(237) = {2, -10, -20, -11}; Plane Surface(238) = {237};
+ Line Loop(239) = {13, 15, 16, 22}; Plane Surface(240) = {239};
+ Line Loop(241) = {12, -22, 17, 19}; Plane Surface(242) = {241};
+ Line Loop(243) = {19, -11, -1, -18}; Plane Surface(244) = {243};
+ Line Loop(245) = {127, 70, -133, -88}; Plane Surface(246) = {245};
+ Line Loop(247) = {133, 71, -136, -89}; Plane Surface(248) = {247};
+ Line Loop(249) = {90, 112, -72, -136}; Plane Surface(250) = {249};
+ Line Loop(251) = {128, -42, -134, -70}; Plane Surface(252) = {251};
+ Line Loop(253) = {41, -134, 71, 137}; Plane Surface(254) = {253};
+ Line Loop(255) = {137, -40, -113, -72}; Plane Surface(256) = {255};
+ Line Loop(257) = {129, -10, -135, 42}; Plane Surface(258) = {257};
+ Line Loop(259) = {135, -9, -138, 41}; Plane Surface(260) = {259};
+ Line Loop(261) = {138, -8, -114, 40}; Plane Surface(262) = {261};
+ Line Loop(263) = {130, -69, -142, 93}; Plane Surface(264) = {263};
+ Line Loop(265) = {92, 142, -68, -139}; Plane Surface(266) = {265};
+ Line Loop(267) = {139, -67, -109, 91}; Plane Surface(268) = {267};
+ Line Loop(269) = {131, 37, -143, 69}; Plane Surface(270) = {269};
+ Line Loop(271) = {143, 38, -140, 68}; Plane Surface(272) = {271};
+ Line Loop(273) = {140, 39, -110, 67}; Plane Surface(274) = {273};
+ Line Loop(275) = {132, 11, -144, -37}; Plane Surface(276) = {275};
+ Line Loop(277) = {12, -141, -38, 144}; Plane Surface(278) = {277};
+ Line Loop(279) = {141, 13, -111, -39}; Plane Surface(280) = {279};
+ Line Loop(281) = {99, -18, -102, 36}; Plane Surface(282) = {281};
+ Line Loop(283) = {102, -17, -105, 35}; Plane Surface(284) = {283};
+ Line Loop(285) = {34, 105, -16, -108}; Plane Surface(286) = {285};
+ Line Loop(287) = {34, -104, 51, 107}; Plane Surface(288) = {287};
+ Line Loop(289) = {50, 104, 35, -101}; Plane Surface(290) = {289};
+ Line Loop(291) = {36, -98, 49, 101}; Plane Surface(292) = {291};
+ Line Loop(293) = {97, 49, -100, -73}; Plane Surface(294) = {293};
+ Line Loop(295) = {74, 103, -50, -100}; Plane Surface(296) = {295};
+ Line Loop(297) = {103, 51, -106, -75}; Plane Surface(298) = {297};
+ Line Loop(299) = {126, 4, -123, -28}; Plane Surface(300) = {299};
+ Line Loop(301) = {125, 28, -122, 57}; Plane Surface(302) = {301};
+ Line Loop(303) = {122, 29, -119, 56}; Plane Surface(304) = {303};
+ Line Loop(305) = {119, 30, -116, 55}; Plane Surface(306) = {305};
+ Line Loop(307) = {55, -118, -79, 115}; Plane Surface(308) = {307};
+ Line Loop(309) = {80, 121, -56, -118}; Plane Surface(310) = {309};
+ Line Loop(311) = {81, 124, -57, -121}; Plane Surface(312) = {311};
+ Line Loop(313) = {123, 5, -120, -29}; Plane Surface(314) = {313};
+ Line Loop(315) = {120, 6, -117, -30}; Plane Surface(316) = {315};
+ Line Loop(317) = {121, -63, -133, 87}; Plane Surface(318) = {317};
+ Line Loop(319) = {133, -62, -142, 86}; Plane Surface(320) = {319};
+ Line Loop(321) = {142, -61, -100, 85}; Plane Surface(322) = {321};
+ Line Loop(323) = {122, -45, -134, 63}; Plane Surface(324) = {323};
+ Line Loop(325) = {134, -44, -143, 62}; Plane Surface(326) = {325};
+ Line Loop(327) = {143, -43, -101, 61}; Plane Surface(328) = {327};
+ Line Loop(329) = {123, -21, -135, 45}; Plane Surface(330) = {329};
+ Line Loop(331) = {135, -20, -144, 44}; Plane Surface(332) = {331};
+ Line Loop(333) = {144, -19, -102, 43}; Plane Surface(334) = {333};
+ Line Loop(335) = {126, -3, -129, 27}; Plane Surface(336) = {335};
+ Line Loop(337) = {129, -2, -132, 26}; Plane Surface(338) = {337};
+ Line Loop(339) = {132, -1, -99, 25}; Plane Surface(340) = {339};
+ Line Loop(341) = {125, -27, -128, -58}; Plane Surface(342) = {341};
+ Line Loop(343) = {128, -26, -131, -59}; Plane Surface(344) = {343};
+ Line Loop(345) = {131, -25, -98, -60}; Plane Surface(346) = {345};
+ Line Loop(347) = {97, -60, -130, 84}; Plane Surface(348) = {347};
+ Line Loop(349) = {83, 130, -59, -127}; Plane Surface(350) = {349};
+ Line Loop(351) = {58, -127, -82, 124}; Plane Surface(352) = {351};
+ Line Loop(353) = {96, 103, -66, -139}; Plane Surface(354) = {353};
+ Line Loop(355) = {139, -65, -136, 95}; Plane Surface(356) = {355};
+ Line Loop(357) = {136, -64, -118, 94}; Plane Surface(358) = {357};
+ Line Loop(359) = {66, 104, -48, -140}; Plane Surface(360) = {359};
+ Line Loop(361) = {140, -47, -137, 65}; Plane Surface(362) = {361};
+ Line Loop(363) = {137, -46, -119, 64}; Plane Surface(364) = {363};
+ Line Loop(365) = {105, 22, -141, 48}; Plane Surface(366) = {365};
+ Line Loop(367) = {141, 23, -138, 47}; Plane Surface(368) = {367};
+ Line Loop(369) = {138, 24, -120, 46}; Plane Surface(370) = {369};
+ Line Loop(371) = {108, -15, -111, 33}; Plane Surface(372) = {371};
+ Line Loop(373) = {111, 14, -114, 32}; Plane Surface(374) = {373};
+ Line Loop(375) = {114, -7, -117, 31}; Plane Surface(376) = {375};
+ Line Loop(377) = {107, -33, -110, -52}; Plane Surface(378) = {377};
+ Line Loop(379) = {110, -32, -113, -53}; Plane Surface(380) = {379};
+ Line Loop(381) = {113, -31, -116, -54}; Plane Surface(382) = {381};
+ Line Loop(383) = {106, 52, -109, -76}; Plane Surface(384) = {383};
+ Line Loop(385) = {109, 53, -112, -77}; Plane Surface(386) = {385};
+ Line Loop(387) = {112, 54, -115, -78}; Plane Surface(388) = {387};
+
+ Surface Loop(393) = {182, 298, 384, 354, 200, 268}; Volume(394) = {393};
+ Surface Loop(395) = {184, 296, 354, 266, 204, 322}; Volume(396) = {395};
+ Surface Loop(397) = {190, 348, 294, 264, 208, 322}; Volume(398) = {397};
+ Surface Loop(399) = {188, 350, 264, 246, 320, 206}; Volume(400) = {399};
+ Surface Loop(401) = {186, 352, 312, 192, 318, 246}; Volume(402) = {401};
+ Surface Loop(403) = {180, 248, 266, 202, 356, 320}; Volume(404) = {403};
+ Surface Loop(405) = {310, 174, 318, 248, 194, 358}; Volume(406) = {405};
+ Surface Loop(407) = {176, 388, 308, 358, 250, 196}; Volume(408) = {407};
+ Surface Loop(409) = {178, 386, 198, 356, 268, 250}; Volume(410) = {409};
+ Surface Loop(411) = {378, 288, 200, 274, 360, 220}; Volume(412) = {411};
+ Surface Loop(413) = {216, 380, 198, 274, 256, 362}; Volume(414) = {413};
+ Surface Loop(415) = {306, 382, 214, 364, 196, 256}; Volume(416) = {415};
+ Surface Loop(417) = {304, 194, 364, 212, 254, 324}; Volume(418) = {417};
+ Surface Loop(419) = {302, 342, 192, 324, 252, 210}; Volume(420) = {419};
+ Surface Loop(421) = {344, 206, 252, 224, 270, 326}; Volume(422) = {421};
+ Surface Loop(423) = {292, 346, 328, 208, 270, 226}; Volume(424) = {423};
+ Surface Loop(425) = {290, 360, 204, 272, 328, 222}; Volume(426) = {425};
+ Surface Loop(427) = {240, 372, 286, 220, 366, 280}; Volume(428) = {427};
+ Surface Loop(429) = {234, 374, 216, 368, 262, 280}; Volume(430) = {429};
+ Surface Loop(431) = {376, 232, 316, 370, 214, 262}; Volume(432) = {431};
+ Surface Loop(433) = {314, 230, 260, 212, 330, 370}; Volume(434) = {433};
+ Surface Loop(435) = {228, 336, 300, 210, 258, 330}; Volume(436) = {435};
+ Surface Loop(437) = {236, 332, 368, 278, 260, 218}; Volume(438) = {437};
+ Surface Loop(439) = {284, 242, 366, 334, 222, 278}; Volume(440) = {439};
+ Surface Loop(441) = {282, 340, 244, 334, 226, 276}; Volume(442) = {441};
+ Surface Loop(443) = {238, 338, 258, 276, 332, 224}; Volume(444) = {443};
+
+ pt_sc=newp;
+ // Printf("newp scatterer=%g",pt_sc);
+ l_sc =newl;
+ ll_sc =newll;
+ s_sc =news;
+
+ If(flag_shape==ELL)
+ Point(pt_sc+20) = { 0, 0, ell_rz, In_lc}; Point(pt_sc+21) = { 0, 0,-ell_rz, In_lc};
+ Point(pt_sc+22) = { 0, ell_ry, 0, In_lc}; Point(pt_sc+23) = { 0,-ell_ry, 0, In_lc};
+ Point(pt_sc+25) = { ell_rx, 0, 0, In_lc}; Point(pt_sc+26) = {-ell_rx, 0, 0, In_lc};
+ Ellipse(l_sc+145) = {pt_sc+22, 73, pt_sc+26, pt_sc+26} Plane{0,0,1};
+ Ellipse(l_sc+146) = {pt_sc+26, 73, pt_sc+23, pt_sc+23} Plane{0,0,1};
+ Ellipse(l_sc+147) = {pt_sc+23, 73, pt_sc+25, pt_sc+25} Plane{0,0,1};
+ Ellipse(l_sc+148) = {pt_sc+25, 73, pt_sc+22, pt_sc+22} Plane{0,0,1};
+ Ellipse(l_sc+149) = {pt_sc+21, 73, pt_sc+26, pt_sc+26} Plane{0,0,1};
+ Ellipse(l_sc+150) = {pt_sc+26, 73, pt_sc+20, pt_sc+20} Plane{0,0,1};
+ Ellipse(l_sc+151) = {pt_sc+20, 73, pt_sc+25, pt_sc+25} Plane{0,0,1};
+ Ellipse(l_sc+152) = {pt_sc+25, 73, pt_sc+21, pt_sc+21} Plane{0,0,1};
+ Ellipse(l_sc+153) = {pt_sc+23, 73, pt_sc+20, pt_sc+20} Plane{0,0,1};
+ Ellipse(l_sc+154) = {pt_sc+20, 73, pt_sc+22, pt_sc+22} Plane{0,0,1};
+ Ellipse(l_sc+155) = {pt_sc+22, 73, pt_sc+21, pt_sc+21} Plane{0,0,1};
+ Ellipse(l_sc+156) = {pt_sc+21, 73, pt_sc+23, pt_sc+23} Plane{0,0,1};
+ Line Loop(ll_sc+157) = {l_sc+151, -147-l_sc, 153+l_sc}; Surface(158+s_sc) = {157+ll_sc};
+ Line Loop(ll_sc+159) = {l_sc+151, 148+l_sc, -154-l_sc}; Surface(160+s_sc) = {159+ll_sc};
+ Line Loop(ll_sc+161) = {l_sc+145, 150+l_sc, 154+l_sc}; Surface(162+s_sc) = {161+ll_sc};
+ Line Loop(ll_sc+163) = {l_sc+150, -153-l_sc, -146-l_sc}; Surface(164+s_sc) = {163+ll_sc};
+ Line Loop(ll_sc+165) = {l_sc+156, 147+l_sc, 152+l_sc}; Surface(166+s_sc) = {165+ll_sc};
+ Line Loop(ll_sc+167) = {l_sc+152, -155-l_sc, -148-l_sc}; Surface(168+s_sc) = {167+ll_sc};
+ Line Loop(ll_sc+169) = {l_sc+155, 149+l_sc, -145-l_sc}; Surface(170+s_sc) = {169+ll_sc};
+ Line Loop(ll_sc+171) = {l_sc+149, 146+l_sc, -156-l_sc}; Surface(172+s_sc) = {171+ll_sc};
+ Surface Loop(1072) = {160+s_sc, 162+s_sc, 164+s_sc, 158+s_sc, 166+s_sc, 170+s_sc, 172+s_sc, 168+s_sc};
+ Volume(1092) = {1072};
+ EndIf
+ If(flag_shape==PARALL)
+ Point(pt_sc+20) = {-par_ax/2,-par_ay/2,-par_az/2,In_lc};
+ Point(pt_sc+21) = { par_ax/2,-par_ay/2,-par_az/2,In_lc};
+ Point(pt_sc+22) = { par_ax/2, par_ay/2,-par_az/2,In_lc};
+ Point(pt_sc+23) = {-par_ax/2, par_ay/2,-par_az/2,In_lc};
+ Point(pt_sc+24) = {-par_ax/2,-par_ay/2, par_az/2,In_lc};
+ Point(pt_sc+25) = { par_ax/2,-par_ay/2, par_az/2,In_lc};
+ Point(pt_sc+26) = { par_ax/2, par_ay/2, par_az/2,In_lc};
+ Point(pt_sc+27) = {-par_ax/2, par_ay/2, par_az/2,In_lc};
+ Line(l_sc+1) = {pt_sc+20, pt_sc+24};
+ Line(l_sc+2) = {pt_sc+24, pt_sc+25};
+ Line(l_sc+3) = {pt_sc+25, pt_sc+21};
+ Line(l_sc+4) = {pt_sc+21, pt_sc+20};
+ Line(l_sc+5) = {pt_sc+23, pt_sc+27};
+ Line(l_sc+6) = {pt_sc+27, pt_sc+26};
+ Line(l_sc+7) = {pt_sc+26, pt_sc+22};
+ Line(l_sc+8) = {pt_sc+22, pt_sc+23};
+ Line(l_sc+9) = {pt_sc+23, pt_sc+20};
+ Line(l_sc+10) = {pt_sc+22, pt_sc+21};
+ Line(l_sc+11) = {pt_sc+26, pt_sc+25};
+ Line(l_sc+12) = {pt_sc+27, pt_sc+24};
+ Line Loop(ll_sc+13) = {l_sc+2, -11-l_sc, -6-l_sc, l_sc+12};Plane Surface(s_sc+14) = {ll_sc+13};
+ Line Loop(ll_sc+15) = {l_sc+1, -12-l_sc, -5-l_sc, l_sc+9};Plane Surface(s_sc+16) = {ll_sc+15};
+ Line Loop(ll_sc+17) = {l_sc+3, -10-l_sc, -7-l_sc, l_sc+11};Plane Surface(s_sc+18) = {ll_sc+17};
+ Line Loop(ll_sc+19) = {l_sc+8, 5+l_sc, 6+l_sc, l_sc+7};Plane Surface(s_sc+20) = {ll_sc+19};
+ Line Loop(ll_sc+21) = {l_sc+4, -9-l_sc, -8-l_sc, l_sc+10};Plane Surface(s_sc+22) = {ll_sc+21};
+ Line Loop(ll_sc+23) = {l_sc+1, 2+l_sc, 3+l_sc, l_sc+4};Plane Surface(s_sc+24) = {ll_sc+23};
+ Surface Loop(1072) = {s_sc+16, 24+s_sc, 14+s_sc, 18+s_sc, 22+s_sc, 20+s_sc};
+ Volume(1092) = {1072};
+ EndIf
+ If(flag_shape==CYL)
+ Point(pt_sc+20) = { 0, cyl_ry, -cyl_h/2, In_lc};
+ Point(pt_sc+21) = { 0,-cyl_ry, -cyl_h/2, In_lc};
+ Point(pt_sc+22) = { cyl_rx, 0, -cyl_h/2, In_lc};
+ Point(pt_sc+23) = {-cyl_rx, 0, -cyl_h/2, In_lc};
+ Point(pt_sc+24) = {0, 0, -cyl_h/2, In_lc};
+ Point(pt_sc+25) = { 0, cyl_ry, cyl_h/2, In_lc};
+ Point(pt_sc+26) = { 0,-cyl_ry, cyl_h/2, In_lc};
+ Point(pt_sc+27) = { cyl_rx, 0, cyl_h/2, In_lc};
+ Point(pt_sc+28) = {-cyl_rx, 0, cyl_h/2, In_lc};
+ Point(pt_sc+29) = {0, 0, cyl_h/2, In_lc};
+ Ellipse(l_sc+1) = {pt_sc+23, pt_sc+24, pt_sc+21, pt_sc+21};
+ Ellipse(l_sc+2) = {pt_sc+21, pt_sc+24, pt_sc+22, pt_sc+22};
+ Ellipse(l_sc+3) = {pt_sc+22, pt_sc+24, pt_sc+20, pt_sc+20};
+ Ellipse(l_sc+4) = {pt_sc+20, pt_sc+24, pt_sc+24, pt_sc+23};
+ Ellipse(l_sc+5) = {pt_sc+28, pt_sc+29, pt_sc+26, pt_sc+26};
+ Ellipse(l_sc+6) = {pt_sc+26, pt_sc+29, pt_sc+27, pt_sc+27};
+ Ellipse(l_sc+7) = {pt_sc+27, pt_sc+29, pt_sc+25, pt_sc+25};
+ Ellipse(l_sc+8) = {pt_sc+25, pt_sc+29, pt_sc+28, pt_sc+28};
+ Line(l_sc+9) = {pt_sc+26, pt_sc+21};
+ Line(l_sc+10) = {pt_sc+28, pt_sc+23};
+ Line(l_sc+11) = {pt_sc+25, pt_sc+20};
+ Line(l_sc+12) = {pt_sc+27, pt_sc+22};
+ Line Loop(ll_sc+13) = {l_sc+8, 5 +l_sc , 6+l_sc, 7+l_sc};Plane Surface(s_sc+14) = {ll_sc+13};
+ Line Loop(ll_sc+15) = {l_sc+4, 1 +l_sc , 2+l_sc, 3+l_sc};Plane Surface(s_sc+16) = {ll_sc+15};
+ Line Loop(ll_sc+17) = {l_sc+8, 10+l_sc , -4-l_sc, -11-l_sc};Surface(s_sc+18) = {ll_sc+17};
+ Line Loop(ll_sc+19) = {l_sc+5, 9 +l_sc , -1-l_sc, -10-l_sc};Surface(s_sc+20) = {ll_sc+19};
+ Line Loop(ll_sc+21) = {l_sc+9, 2 +l_sc , -12-l_sc, -6-l_sc};Surface(s_sc+22) = {ll_sc+21};
+ Line Loop(ll_sc+23) = {l_sc+3, -11-l_sc , -7-l_sc, 12+l_sc};Surface(s_sc+24) = {ll_sc+23};
+ Surface Loop(1072) = {s_sc+14, s_sc+18, s_sc+20, s_sc+22, s_sc+16, s_sc+24};
+ Volume(1092) = {1072};
+ EndIf
+ If(flag_shape==CONE)
+ Point(pt_sc+20) = { 0, cone_ry, -cone_h/3, In_lc};
+ Point(pt_sc+21) = { 0,-cone_ry, -cone_h/3, In_lc};
+ Point(pt_sc+22) = { cone_rx, 0, -cone_h/3, In_lc};
+ Point(pt_sc+23) = {-cone_rx, 0, -cone_h/3, In_lc};
+ Point(pt_sc+24) = {0, 0, -cone_h/3, In_lc};
+ Point(pt_sc+26) = {0, 0, 2*cone_h/3, In_lc};
+ Ellipse(l_sc+1) = {pt_sc+23, pt_sc+24, pt_sc+21, pt_sc+21};
+ Ellipse(l_sc+2) = {pt_sc+21, pt_sc+24, pt_sc+22, pt_sc+22};
+ Ellipse(l_sc+3) = {pt_sc+22, pt_sc+24, pt_sc+20, pt_sc+20};
+ Ellipse(l_sc+4) = {pt_sc+20, pt_sc+24, pt_sc+24, pt_sc+23};
+ Line(l_sc+5)={pt_sc+20,pt_sc+26};Line(l_sc+6)={pt_sc+21,pt_sc+26};
+ Line(l_sc+7)={pt_sc+22,pt_sc+26};Line(l_sc+8)={pt_sc+23,pt_sc+26};
+ Line Loop(ll_sc+9) = {l_sc+1, l_sc+2, l_sc+3, l_sc+4};Plane Surface(s_sc+10) = {ll_sc+9};
+ Line Loop(ll_sc+11) = {l_sc+6, -8-l_sc, 1+l_sc};Surface(s_sc+12) = {ll_sc+11};
+ Line Loop(ll_sc+13) = {l_sc+8, -5-l_sc, 4+l_sc};Surface(s_sc+14) = {ll_sc+13};
+ Line Loop(ll_sc+15) = {l_sc+5, -7-l_sc, 3+l_sc};Surface(s_sc+16) = {ll_sc+15};
+ Line Loop(ll_sc+17) = {l_sc+7, -6-l_sc, 2+l_sc};Surface(s_sc+18) = {ll_sc+17};
+ Surface Loop(1072) = {s_sc+12, s_sc+18, s_sc+16, s_sc+14, s_sc+10};
+ Volume(1092) = {1072};
+ EndIf
+ If (flag_shape==TOR)
+ Point(pt_sc+20) = { 0, tor_r1+tor_r2x+tor_r2x , 0, In_lc};
+ Point(pt_sc+21) = { 0, tor_r1+tor_r2x-tor_r2x , 0, In_lc};
+ Point(pt_sc+22) = { 0, tor_r1+tor_r2x , tor_r2z, In_lc};
+ Point(pt_sc+23) = { 0, tor_r1+tor_r2x , -tor_r2z, In_lc};
+ Point(pt_sc+24) = { 0, tor_r1+tor_r2x, 0, In_lc};
+ Ellipse(l_sc+1) = {pt_sc+21, pt_sc+24, pt_sc+22, pt_sc+22};
+ Ellipse(l_sc+2) = {pt_sc+22, pt_sc+24, pt_sc+20, pt_sc+20};
+ Ellipse(l_sc+3) = {pt_sc+20, pt_sc+24, pt_sc+23, pt_sc+23};
+ Ellipse(l_sc+4) = {pt_sc+23, pt_sc+24, pt_sc+21, pt_sc+21};
+ Line Loop(ll_sc+5) = {l_sc+2, l_sc+3, l_sc+4, l_sc+1};Plane Surface(s_sc+6) = {ll_sc+5};
+ Extrude {{0, 0, 1}, {0, 0, 0}, deg2rad*tor_angle/2} {Surface{s_sc+6};}
+ Extrude {{0, 0, 1}, {0, 0, 0},-deg2rad*tor_angle/2} {Surface{s_sc+6};}
+ Surface Loop(1072) = {472, 460, 464, 468, 473, 494, 482, 486, 490, 495};
+ EndIf
+ inerpml_ll=newll;
+ Surface Loop(inerpml_ll) = {326, 272, 218, 362, 202, 254};
+ Volume(1093) = {inerpml_ll,1072};
+
+ Physical Point(2000) = {1}; // PrintPoint
+ Physical Volume(1000) = {402, 398, 394, 408, 442, 428, 432, 436}; // PML XYZ
+ Physical Volume(1001) = {400, 410, 430, 444}; // PML XZ
+ Physical Volume(1002) = {406, 396, 434, 440}; // PML YZ
+ Physical Volume(1003) = {412, 416, 420, 424}; // PML XY
+ Physical Volume(1004) = {404, 438}; // PML Z
+ Physical Volume(1005) = {426, 418}; // PML Y
+ Physical Volume(1006) = {414, 422}; // PML X
+ Physical Volume(1007) = {1093}; // Out
+ If (flag_shape==TOR)
+ Physical Volume(1008) = {445,446}; // Scatterer
+ Else
+ Physical Volume(1008) = {1092}; // Scatterer
+ EndIf
+
+ // Physical Volume(1007) = {392}; // Out
+ // Physical Volume(1008) = {1039,1041,1043,1045,1047,1049,1051,1053}; // Scatterer
+EndIf
+If(flag_cartpml==0)
+ Point(1) = { r_pml_in, 0, 0 , Out_lc};
+ Point(2) = {-r_pml_in, 0, 0 , Out_lc};
+ Point(3) = {0, r_pml_in, 0 , Out_lc};
+ Point(4) = {0, -r_pml_in, 0 , Out_lc};
+ Point(5) = {0, 0, r_pml_in , Out_lc};
+ Point(6) = {0, 0, -r_pml_in , Out_lc};
+ Point(7) = { r_pml_out, 0, 0, PML_lc};
+ Point(8) = {-r_pml_out, 0, 0, PML_lc};
+ Point(9) = {0, r_pml_out, 0, PML_lc};
+ Point(10) = {0, -r_pml_out, 0, PML_lc};
+ Point(11) = {0, 0, r_pml_out, PML_lc};
+ Point(12) = {0, 0, -r_pml_out, PML_lc};
+ Point(13) = { 0, 0, 0, CenterScat_lc};
+ If(flag_shape==ELL)
+ Point(20) = { 0, 0, ell_rz, In_lc}; Point(21) = { 0, 0,-ell_rz, In_lc};
+ Point(22) = { 0, ell_ry, 0, In_lc}; Point(23) = { 0,-ell_ry, 0, In_lc};
+ Point(25) = { ell_rx, 0, 0, In_lc}; Point(26) = {-ell_rx, 0, 0, In_lc};
+ Ellipse(145) = {22, 13, 26, 26} Plane{0,0,1};
+ Ellipse(146) = {26, 13, 23, 23} Plane{0,0,1};
+ Ellipse(147) = {23, 13, 25, 25} Plane{0,0,1};
+ Ellipse(148) = {25, 13, 22, 22} Plane{0,0,1};
+ Ellipse(149) = {21, 13, 26, 26} Plane{0,0,1};
+ Ellipse(150) = {26, 13, 20, 20} Plane{0,0,1};
+ Ellipse(151) = {20, 13, 25, 25} Plane{0,0,1};
+ Ellipse(152) = {25, 13, 21, 21} Plane{0,0,1};
+ Ellipse(153) = {23, 13, 20, 20} Plane{0,0,1};
+ Ellipse(154) = {20, 13, 22, 22} Plane{0,0,1};
+ Ellipse(155) = {22, 13, 21, 21} Plane{0,0,1};
+ Ellipse(156) = {21, 13, 23, 23} Plane{0,0,1};
+ Line Loop(157) = {151, -147, 153}; Surface(158) = {157};
+ Line Loop(159) = {151, 148, -154}; Surface(160) = {159};
+ Line Loop(161) = {145, 150, 154}; Surface(162) = {161};
+ Line Loop(163) = {150, -153, -146}; Surface(164) = {163};
+ Line Loop(165) = {156, 147, 152}; Surface(166) = {165};
+ Line Loop(167) = {152, -155, -148}; Surface(168) = {167};
+ Line Loop(169) = {155, 149, -145}; Surface(170) = {169};
+ Line Loop(171) = {149, 146, -156}; Surface(172) = {171};
+ Surface Loop(1072) = {160, 162, 164, 158, 166, 170, 172, 168};
+ Volume(1092) = {1072};
+ EndIf
+ If(flag_shape==PARALL)
+ Point(20) = {-par_ax/2,-par_ay/2,-par_az/2,In_lc};
+ Point(21) = { par_ax/2,-par_ay/2,-par_az/2,In_lc};
+ Point(22) = { par_ax/2, par_ay/2,-par_az/2,In_lc};
+ Point(23) = {-par_ax/2, par_ay/2,-par_az/2,In_lc};
+ Point(24) = {-par_ax/2,-par_ay/2, par_az/2,In_lc};
+ Point(25) = { par_ax/2,-par_ay/2, par_az/2,In_lc};
+ Point(26) = { par_ax/2, par_ay/2, par_az/2,In_lc};
+ Point(27) = {-par_ax/2, par_ay/2, par_az/2,In_lc};
+ Line(1) = {20, 24};
+ Line(2) = {24, 25};
+ Line(3) = {25, 21};
+ Line(4) = {21, 20};
+ Line(5) = {23, 27};
+ Line(6) = {27, 26};
+ Line(7) = {26, 22};
+ Line(8) = {22, 23};
+ Line(9) = {23, 20};
+ Line(10) = {22, 21};
+ Line(11) = {26, 25};
+ Line(12) = {27, 24};
+ Line Loop(13) = {2, -11, -6, 12};Plane Surface(14) = {13};
+ Line Loop(15) = {1, -12, -5, 9};Plane Surface(16) = {15};
+ Line Loop(17) = {3, -10, -7, 11};Plane Surface(18) = {17};
+ Line Loop(19) = {8, 5, 6, 7};Plane Surface(20) = {19};
+ Line Loop(21) = {4, -9, -8, 10};Plane Surface(22) = {21};
+ Line Loop(23) = {1, 2, 3, 4};Plane Surface(24) = {23};
+ Surface Loop(1072) = {16, 24, 14, 18, 22, 20};
+ Volume(1092) = {1072};
+ EndIf
+ If(flag_shape==CYL)
+ Point(20) = { 0, cyl_ry, -cyl_h/2, In_lc};
+ Point(21) = { 0,-cyl_ry, -cyl_h/2, In_lc};
+ Point(22) = { cyl_rx, 0, -cyl_h/2, In_lc};
+ Point(23) = {-cyl_rx, 0, -cyl_h/2, In_lc};
+ Point(24) = {0, 0, -cyl_h/2, In_lc};
+ Point(25) = { 0, cyl_ry, cyl_h/2, In_lc};
+ Point(26) = { 0,-cyl_ry, cyl_h/2, In_lc};
+ Point(27) = { cyl_rx, 0, cyl_h/2, In_lc};
+ Point(28) = {-cyl_rx, 0, cyl_h/2, In_lc};
+ Point(29) = {0, 0, cyl_h/2, In_lc};
+ Ellipse(1) = {23, 24, 21, 21};
+ Ellipse(2) = {21, 24, 22, 22};
+ Ellipse(3) = {22, 24, 20, 20};
+ Ellipse(4) = {20, 24, 24, 23};
+ Ellipse(5) = {28, 29, 26, 26};
+ Ellipse(6) = {26, 29, 27, 27};
+ Ellipse(7) = {27, 29, 25, 25};
+ Ellipse(8) = {25, 29, 28, 28};
+ Line(9) = {26, 21};
+ Line(10) = {28, 23};
+ Line(11) = {25, 20};
+ Line(12) = {27, 22};
+ Line Loop(13) = {8, 5, 6, 7};Plane Surface(14) = {13};
+ Line Loop(15) = {4, 1, 2, 3};Plane Surface(16) = {15};
+ Line Loop(17) = {8, 10, -4, -11};Surface(18) = {17};
+ Line Loop(19) = {5, 9, -1, -10};Surface(20) = {19};
+ Line Loop(21) = {9, 2, -12, -6};Surface(22) = {21};
+ Line Loop(23) = {3, -11, -7, 12};Surface(24) = {23};
+ Surface Loop(1072) = {14, 18, 20, 22, 16, 24};
+ Volume(1092) = {1072};
+ EndIf
+ If(flag_shape==CONE)
+ Point(20) = { 0, cone_ry, -cone_h/3, In_lc};
+ Point(21) = { 0,-cone_ry, -cone_h/3, In_lc};
+ Point(22) = { cone_rx, 0, -cone_h/3, In_lc};
+ Point(23) = {-cone_rx, 0, -cone_h/3, In_lc};
+ Point(24) = {0, 0, -cone_h/3, In_lc};
+ Point(26) = {0, 0, 2*cone_h/3, In_lc};
+ Ellipse(1) = {23, 24, 21, 21};
+ Ellipse(2) = {21, 24, 22, 22};
+ Ellipse(3) = {22, 24, 20, 20};
+ Ellipse(4) = {20, 24, 24, 23};
+ Line(5)={20,26};Line(6)={21,26};
+ Line(7)={22,26};Line(8)={23,26};
+ Line Loop(9) = {1, 2, 3, 4};Plane Surface(10) = {9};
+ Line Loop(11) = {6, -8, 1};Surface(12) = {11};
+ Line Loop(13) = {8, -5, 4};Surface(14) = {13};
+ Line Loop(15) = {5, -7, 3};Surface(16) = {15};
+ Line Loop(17) = {7, -6, 2};Surface(18) = {17};
+ Surface Loop(1072) = {12, 18, 16, 14, 10};
+ Volume(1092) = {1072};
+ EndIf
+ If (flag_shape==TOR)
+ Point(20) = { 0, tor_r1+tor_r2x+tor_r2x , 0, In_lc};
+ Point(21) = { 0, tor_r1+tor_r2x-tor_r2x , 0, In_lc};
+ Point(22) = { 0, tor_r1+tor_r2x , tor_r2z, In_lc};
+ Point(23) = { 0, tor_r1+tor_r2x , -tor_r2z, In_lc};
+ Point(24) = { 0, tor_r1+tor_r2x, 0, In_lc};
+ Ellipse(1) = {21, 24, 22, 22};
+ Ellipse(2) = {22, 24, 20, 20};
+ Ellipse(3) = {20, 24, 23, 23};
+ Ellipse(4) = {23, 24, 21, 21};
+ Line Loop(5) = {2, 3, 4, 1};Plane Surface(6) = {5};
+ Extrude {{0, 0, 1}, {0, 0, 0}, deg2rad*tor_angle/2} {Surface{6};}
+ Extrude {{0, 0, 1}, {0, 0, 0},-deg2rad*tor_angle/2} {Surface{6};}
+ Surface Loop(1072) = {37, 41, 45, 49, 50, 15, 19, 23, 27, 28};
+ EndIf
+ // Rotate {{0, 0, 1}, {0, 0, 0}, rot_phi*deg2rad} {Volume{1092};}
+ // Rotate {{-Sin[rot_phi*deg2rad], Cos[rot_phi*deg2rad], 0}, {0, 0, 0}, rot_theta*deg2rad} {Volume{1092};}
+
+ Circle(114) = {2 , 13, 4};
+ Circle(115) = {4 , 13, 1};
+ Circle(116) = {1 , 13, 3};
+ Circle(117) = {3 , 13, 2};
+ Circle(118) = {5 , 13, 3};
+ Circle(119) = {3 , 13, 6};
+ Circle(120) = {6 , 13, 4};
+ Circle(121) = {4 , 13, 5};
+ Circle(122) = {6 , 13, 2};
+ Circle(123) = {2 , 13, 5};
+ Circle(124) = {5 , 13, 1};
+ Circle(125) = {1 , 13, 6};
+ Circle(126) = {8 , 13, 10};
+ Circle(127) = {10, 13, 7};
+ Circle(128) = {7 , 13, 9};
+ Circle(129) = {9 , 13, 8};
+ Circle(130) = {8 , 13, 11};
+ Circle(131) = {11, 13, 7};
+ Circle(132) = {7 , 13, 12};
+ Circle(133) = {12, 13, 8};
+ Circle(134) = {11, 13, 10};
+ Circle(135) = {10, 13, 12};
+ Circle(136) = {12, 13, 9};
+ Circle(137) = {9 , 13, 11};
+
+ Line Loop(1054) = {123, 118, 117};Surface(1055) = {1054};
+ Line Loop(1056) = {118, -116, -124};Surface(1057) = {1056};
+ Line Loop(1058) = {116, 119, -125};Surface(1059) = {1058};
+ Line Loop(1061) = {119, 122, -117};Surface(1062) = {1061};
+ Line Loop(1063) = {114, -120, 122};Surface(1064) = {1063};
+ Line Loop(1065) = {114, 121, -123};Surface(1066) = {1065};
+ Line Loop(1067) = {121, 124, -115};Surface(1068) = {1067};
+ Line Loop(1069) = {125, 120, 115};Surface(1070) = {1069};
+
+ Line Loop(1074) = {128, 137, 131}; Surface(1075) = {1074};
+ Line Loop(1076) = {129, 130, -137}; Surface(1077) = {1076};
+ Line Loop(1078) = {136, 129, -133}; Surface(1079) = {1078};
+ Line Loop(1080) = {128, -136, -132}; Surface(1081) = {1080};
+ Line Loop(1082) = {134, 127, -131}; Surface(1083) = {1082};
+ Line Loop(1084) = {127, 132, -135}; Surface(1085) = {1084};
+ Line Loop(1086) = {135, 133, 126}; Surface(1087) = {1086};
+ Line Loop(1088) = {126, -134, -130}; Surface(1089) = {1088};
+ // pml in
+ Surface Loop(1100) = {1057, 1055, 1066, 1064, 1070, 1059, 1062, 1068};
+ // pml out
+ Surface Loop(1101) = {1089, 1087, 1085, 1083, 1075, 1081, 1079, 1077};
+
+
+ Surface Loop(1071) = {1062, 1059, 1057, 1068, 1066, 1064, 1070, 1055};
+ Volume(1073) = {1071, 1072};
+ Surface Loop(1090) = {1075, 1081, 1079, 1077, 1089, 1087, 1085, 1083};
+ Volume(1091) = {1071, 1090};
+ If (flag_shape==TOR)
+ Physical Volume(1) = {1,2}; // Scatterer
+ Else
+ Physical Volume(1) = {1092}; // Scatterer
+ EndIf
+ Physical Volume(2) = {1073}; // Out - air
+ Physical Volume(3) = {1091}; // pml
+ Physical Point(2000) = {1}; // PrintPoint
+EndIf
+
+// Mesh.Algorithm = 1; // // 1=MeshAdapt, 5=Delaunay, 6=Frontal
+// Mesh.Algorithm3D = 1;//4; // // 1=Delaunay, 4=Frontal
+Mesh.Optimize = 1;
+Solver.AutoMesh=2;
+Geometry.Points = 1;
+// Mesh.SurfaceEdges = 0;
+Mesh.VolumeEdges = 0;
+
+Printf(Sprintf("nm = %.9e;",nm)) > "scattering_tmp.py";
+Printf(Sprintf("lambda0 = %.9e;",lambda0)) >> "scattering_tmp.py";
+Printf(Sprintf("rbb = %.9e;",rbb)) >> "scattering_tmp.py";
+Printf(Sprintf("epsr_In_re = %.9e;",epsr_In_re)) >> "scattering_tmp.py";
+Printf(Sprintf("epsr_In_im = %.9e;",epsr_In_im)) >> "scattering_tmp.py";
+Printf(Sprintf("epsr_Out_re = %.9e;",epsr_Out_re)) >> "scattering_tmp.py";
+Printf(Sprintf("epsr_Out_im = %.9e;",epsr_Out_im)) >> "scattering_tmp.py";
+Printf(Sprintf("sph_scan = %.9e;",sph_scan)) >> "scattering_tmp.py";
+Printf(Sprintf("r_sph_min = %.9e;",r_sph_min)) >> "scattering_tmp.py";
+Printf(Sprintf("r_sph_max = %.9e;",r_sph_max)) >> "scattering_tmp.py";
+Printf(Sprintf("nb_cuts = %g;" ,nb_cuts)) >> "scattering_tmp.py";
+Printf(Sprintf("npts_theta = %g;" ,npts_theta)) >> "scattering_tmp.py";
+Printf(Sprintf("npts_phi = %g;" ,npts_phi)) >> "scattering_tmp.py";
+Printf(Sprintf("flag_study = %g;" ,flag_study)) >> "scattering_tmp.py";
+Printf(Sprintf("n_max = %g;" ,n_max)) >> "scattering_tmp.py";
+Printf(Sprintf("p_max = %g;" ,p_max)) >> "scattering_tmp.py";
+Printf(Sprintf("siwt = %g;" ,siwt)) >> "scattering_tmp.py";
+
diff --git a/ElectromagneticScattering/scattering.pro b/ElectromagneticScattering/scattering.pro
new file mode 100644
index 0000000000000000000000000000000000000000..98159da821037281a6b4f75d9955f13db5d9f2ee
--- /dev/null
+++ b/ElectromagneticScattering/scattering.pro
@@ -0,0 +1,659 @@
+///////////////////////////////
+// Author : Guillaume Demesy //
+// scattering.pro //
+///////////////////////////////
+
+Include "scattering_data.geo";
+myDir = "run_results/";
+
+Group {
+ If(flag_cartpml==1)
+ PMLxyz = Region[1000];
+ PMLxz = Region[1001];
+ PMLyz = Region[1002];
+ PMLxy = Region[1003];
+ PMLz = Region[1004];
+ PMLy = Region[1005];
+ PMLx = Region[1006];
+ Scat_In = Region[1008];
+ Scat_Out = Region[1007];
+ // Domains
+ Domain = Region[{Scat_In,Scat_Out}];
+ PMLs = Region[{PMLxyz,PMLxy,PMLxz,PMLyz,PMLx,PMLy,PMLz}];
+ All_domains = Region[{Scat_In,Scat_Out,PMLxyz,PMLxy,PMLxz,PMLyz,PMLx,PMLy,PMLz}];
+ EndIf
+ If(flag_cartpml==0)
+ Scat_In = Region[1];
+ Scat_Out = Region[2];
+ Domain = Region[{Scat_In,Scat_Out}];
+ PMLs = Region[3];
+ All_domains = Region[{Scat_In,Scat_Out,PMLs}];
+ EndIf
+ PrintPoint = Region[2000];
+}
+
+Function{
+ // test={1.4,2.6,3.7};
+ // For i In {0:#test()-1}
+ // Printf("%g",test(i));
+ // EndFor
+ I[] = Complex[0,1];
+ avoid_sing = 1.e-14;
+ mu0 = 4*Pi*100.0*nm;
+ epsilon0 = 8.854187817e-3*nm;
+ cel = 1.0/Sqrt[epsilon0 * mu0];
+ Freq = cel/lambda0;
+ omega0 = 2.0*Pi*Freq;
+ k_Out = 2.0*Pi*Sqrt[epsr_Out_re]/lambda0;
+ r3D_sph[] = Sqrt[X[]^2+Y[]^2+Z[]^2];
+ r2D_sph[] = Sqrt[X[]^2+Y[]^2];
+ cos_theta[] = Z[]/(r3D_sph[]+avoid_sing);
+ theta[] = Acos[cos_theta[]];
+ phi[] = Atan2[-Y[],-X[]]+Pi;
+
+ For kr In {0:nb_cuts-1}
+ r_sph~{kr}=r_sph_min+kr*(r_sph_max-r_sph_min)/(nb_cuts-1);
+ EndFor
+
+ If (flag_study==0)
+ Ae = 1.0;
+ Ah = Ae*Sqrt[epsilon0*epsr_In_re/mu0];
+ alpha0 = k_Out*Sin[theta0]*Cos[phi0];
+ beta0 = k_Out*Sin[theta0]*Sin[phi0];
+ gamma0 = k_Out*Cos[theta0];
+ Ex0 = Ae * Cos[psi0]*Cos[theta0]*Cos[phi0] - Ae* Sin[psi0]*Sin[phi0];
+ Ey0 = Ae * Cos[psi0]*Cos[theta0]*Sin[phi0] + Ae* Sin[psi0]*Cos[phi0];
+ Ez0 = -Ae * Cos[psi0]*Sin[theta0];
+ Hx0 = -1/(omega0*mu0)*(beta0 * Ez0 - gamma0 * Ey0);
+ Hy0 = -1/(omega0*mu0)*(gamma0 * Ex0 - alpha0 * Ez0);
+ Hz0 = -1/(omega0*mu0)*(alpha0 * Ey0 - beta0 * Ex0);
+ Prop[] = Ae * Complex[ Cos[alpha0*X[]+beta0*Y[]+gamma0*Z[]] , siwt*Sin[alpha0*X[]+beta0*Y[]+gamma0*Z[]] ];
+ Einc[] = Vector[Ex0*Prop[],Ey0*Prop[],Ez0*Prop[]];
+ Hinc[] = Vector[Hx0*Prop[],Hy0*Prop[],Hz0*Prop[]];
+ Pinc = 0.5*Ae*Ae*Sqrt[epsilon0*epsr_In_re/mu0] * Cos[theta0];
+ EndIf
+ If (flag_study==2)
+ // Green Dyadic
+ normrr_p[] = Sqrt[(X[]-x_p)^2+(Y[]-y_p)^2+(Z[]-z_p)^2];
+ Gsca[] = Complex[Cos[-1.*siwt*k_Out*normrr_p[]],
+ Sin[-1.*siwt*k_Out*normrr_p[]]]
+ /(4.*Pi*normrr_p[]);
+ GDfact_diag[] = 1.+(I[]*k_Out*normrr_p[]-1.)/(k_Out*normrr_p[])^2;
+ GDfact_glob[] = (3.-3.*I[]*k_Out*normrr_p[]-(k_Out*normrr_p[])^2)/
+ (k_Out*normrr_p[]^2)^2;
+ Dyad_Green2[]=
+ Gsca[]*(
+ GDfact_diag[]*TensorDiag[1,1,1]
+ + GDfact_glob[]*
+ Tensor[(X[]-x_p)*(X[]-x_p),(X[]-x_p)*(Y[]-y_p),(X[]-x_p)*(Z[]-z_p),
+ (Y[]-y_p)*(X[]-x_p),(Y[]-y_p)*(Y[]-y_p),(Y[]-y_p)*(Z[]-z_p),
+ (Z[]-z_p)*(X[]-x_p),(Z[]-z_p)*(Y[]-y_p),(Z[]-z_p)*(Z[]-z_p)]);
+ // Getdp Func
+ Dyad_Green[] = DyadGreenHom[x_p,y_p,z_p,XYZ[],k_Out,siwt];
+ Einc_Green_0[] = Dyad_Green2[]*Vector[1,0,0];
+ Einc_Green_1[] = Dyad_Green2[]*Vector[0,1,0];
+ Einc_Green_2[] = Dyad_Green2[]*Vector[0,0,1];
+ CurlDyad_Green[] = CurlDyadGreenHom[x_p,y_p,z_p,XYZ[],k_Out,siwt];
+ Hinc_Green_0[] = siwt*I[]/(mu0*omega0)*CurlDyad_Green[]*Vector[1,0,0];
+ Hinc_Green_1[] = siwt*I[]/(mu0*omega0)*CurlDyad_Green[]*Vector[0,1,0];
+ Hinc_Green_2[] = siwt*I[]/(mu0*omega0)*CurlDyad_Green[]*Vector[0,0,1];
+ // test[] = DyadGreenHom[x_p,y_p,z_p,x_p+1e-1*nm,y_p+1e-1*nm,z_p+1e-1*nm,lambda0,Sqrt[epsr_Out_re]]*Vector[1,0,0];
+ EndIf
+
+ a_pml = 1.;
+ b_pml = -siwt;
+
+ eigtarget = (2.*Pi*cel/lambda0)^2;
+ eigfilter = 1e-4*Sqrt[eigtarget];
+ EigFilter[] = (Norm[$EigenvalueReal] > eigfilter);
+
+ If(flag_cartpml==1)
+ sx[Scat_In] = 1.;
+ sy[Scat_In] = 1.;
+ sz[Scat_In] = 1.;
+ sx[Scat_Out] = 1.;
+ sy[Scat_Out] = 1.;
+ sz[Scat_Out] = 1.;
+ sx[PMLxyz] = Complex[a_pml,b_pml];
+ sy[PMLxyz] = Complex[a_pml,b_pml];
+ sz[PMLxyz] = Complex[a_pml,b_pml];
+ sx[PMLxz] = Complex[a_pml,b_pml];
+ sy[PMLxz] = 1.0;
+ sz[PMLxz] = Complex[a_pml,b_pml];
+ sx[PMLyz] = 1.0;
+ sy[PMLyz] = Complex[a_pml,b_pml];
+ sz[PMLyz] = Complex[a_pml,b_pml];
+ sx[PMLxy] = Complex[a_pml,b_pml];
+ sy[PMLxy] = Complex[a_pml,b_pml];
+ sz[PMLxy] = 1.0;
+ sx[PMLx] = Complex[a_pml,b_pml];
+ sy[PMLx] = 1.0;
+ sz[PMLx] = 1.0;
+ sx[PMLy] = 1.0;
+ sy[PMLy] = Complex[a_pml,b_pml];
+ sz[PMLy] = 1.0;
+ sx[PMLz] = 1.0;
+ sy[PMLz] = 1.0;
+ sz[PMLz] = Complex[a_pml,b_pml];
+ Lxx[] = sy[]*sz[]/sx[];
+ Lyy[] = sz[]*sx[]/sy[];
+ Lzz[] = sx[]*sy[]/sz[];
+
+ epsilonr_In[] = Complex[epsr_In_re , epsr_In_im];
+ epsilonr_Out[] = Complex[epsr_Out_re , epsr_Out_im];
+
+ epsilonr[Scat_In] = epsilonr_In[] * TensorDiag[1.,1.,1.];
+ epsilonr[Scat_Out] = epsilonr_Out[] * TensorDiag[1.,1.,1.];
+ epsilonr[PMLs] = epsr_Out_re * TensorDiag[Lxx[],Lyy[],Lzz[]];
+
+ epsilonr1[Scat_In] = epsilonr_Out[] * TensorDiag[1.,1.,1.];
+ epsilonr1[Scat_Out] = epsilonr_Out[] * TensorDiag[1.,1.,1.];
+ epsilonr1[PMLs] = epsr_Out_re * TensorDiag[Lxx[],Lyy[],Lzz[]];
+
+ mur[Scat_In] = TensorDiag[1.,1.,1.];
+ mur[Scat_Out] = TensorDiag[1.,1.,1.];
+ mur[PMLs] = TensorDiag[Lxx[],Lyy[],Lzz[]];
+ EndIf
+
+ If(flag_cartpml==0)
+ sr[] = Complex[a_pml,b_pml];
+ sphi[] = Complex[1,0];
+ stheta[] = Complex[1,0];
+ r_tild[] = r_pml_in + (r3D_sph[] - r_pml_in) * sr[];
+ theta_tild[] = theta[];
+ pml_tens_temp1[] = TensorDiag[(r_tild[]/r3D_sph[])^2 * sphi[]*stheta[]/sr[],
+ (r_tild[]/r3D_sph[]) * sr[]*stheta[]/sphi[],
+ (r_tild[]/r3D_sph[]) * sphi[]*sr[]/stheta[]];
+ pml_tens_temp2[] = Rotate[pml_tens_temp1[],0,-theta[]-Pi/2,0];
+ pml_tens[] = Rotate[pml_tens_temp2[],0,0,-phi[]];
+
+ epsilonr_In[] = Complex[epsr_In_re , epsr_In_im];
+ epsilonr_Out[] = Complex[epsr_Out_re , epsr_Out_im];
+
+ epsilonr[Scat_In] = epsilonr_In[] * TensorDiag[1.,1.,1.];
+ epsilonr[Scat_Out] = epsilonr_Out[] * TensorDiag[1.,1.,1.];
+ epsilonr[PMLs] = pml_tens[];
+
+ epsilonr1[Scat_In] = epsilonr_Out[] * TensorDiag[1.,1.,1.];
+ epsilonr1[Scat_Out] = epsilonr_Out[] * TensorDiag[1.,1.,1.];
+ epsilonr1[PMLs] = pml_tens[];
+
+ mur[Scat_In] = TensorDiag[1.,1.,1.];
+ mur[Scat_Out] = TensorDiag[1.,1.,1.];
+ mur[PMLs] = pml_tens[];
+ EndIf
+
+ If(flag_study==RES_PW)
+ source[] = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Einc[];
+ EndIf
+ If(flag_study==RES_TMAT)
+ For pe In {1:p_max}
+ ne = Floor[Sqrt[pe]];
+ me = ne*(ne+1) - Floor[pe];
+ Mnm_source~{pe}[] = Mnm[1,ne,me,XYZ[],k_Out];
+ Nnm_source~{pe}[] = Nnm[1,ne,me,XYZ[],k_Out];
+ source_M~{pe}[] = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Mnm_source~{pe}[];
+ source_N~{pe}[] = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Nnm_source~{pe}[];
+ EndFor
+ EndIf
+
+ If (flag_study==RES_GREEN)
+ sourceG_0[] = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Einc_Green_0[];
+ sourceG_1[] = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Einc_Green_1[];
+ sourceG_2[] = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Einc_Green_2[];
+ EndIf
+}
+
+Constraint {
+ // {Name Dirichlet; Type Assign;
+ // Case {
+ // { Region SurfDirichlet; Value 0.; }
+ // }
+ // }
+}
+
+Jacobian {
+ { Name JVol ;
+ Case {
+ { Region All ; Jacobian Vol ; }
+ }
+ }
+ { Name JSur ;
+ Case {
+ { Region All ; Jacobian Sur ; }
+ }
+ }
+ { Name JLin ;
+ Case {
+ { Region All ; Jacobian Lin ; }
+ }
+ }
+}
+
+Integration {
+ { Name Int_1 ;
+ Case {
+ { Type Gauss ;
+ Case {
+ { GeoElement Point ; NumberOfPoints 1 ; }
+ { GeoElement Line ; NumberOfPoints 4 ; }
+ { GeoElement Triangle ; NumberOfPoints 6 ; }
+ { GeoElement Tetrahedron ; NumberOfPoints 15 ; }
+ }
+ }
+ }
+ }
+}
+
+FunctionSpace {
+ { Name Hcurl; Type Form1;
+ BasisFunction {
+ { Name sn; NameOfCoef un; Function BF_Edge;
+ Support Region[All_domains]; Entity EdgesOf[All]; }
+ { Name sn2; NameOfCoef un2; Function BF_Edge_2E;
+ Support Region[All_domains]; Entity EdgesOf[All]; }
+ If (is_FEM_o2==1)
+ { Name sn3; NameOfCoef un3; Function BF_Edge_3F_b;
+ Support Region[All_domains]; Entity FacetsOf[All]; }
+ { Name sn4; NameOfCoef un4; Function BF_Edge_3F_c;
+ Support Region[All_domains]; Entity FacetsOf[All]; }
+ { Name sn5; NameOfCoef un5; Function BF_Edge_4E;
+ Support Region[All_domains]; Entity EdgesOf[All]; }
+ EndIf
+ }
+ Constraint {
+ // { NameOfCoef un; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+ // { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+ }
+ }
+}
+
+Formulation {
+ If (flag_study==RES_QNM)
+ {Name QNM_helmholtz_vector; Type FemEquation;
+ Quantity {
+ { Name u; Type Local; NameOfSpace Hcurl;}
+ }
+ Equation {
+ Galerkin { [ -1/mur[] * Dof{Curl u} , {Curl u}]; In All_domains; Jacobian JVol; Integration Int_1; }
+ Galerkin { DtDtDof[-1. * 1/cel^2 * epsilonr[]*Dof{u} , {u} ]; In All_domains; Jacobian JVol; Integration Int_1; }
+ }
+ }
+ EndIf
+
+ If (flag_study==RES_PW)
+ {Name PW_helmholtz_vector; Type FemEquation;
+ Quantity {
+ { Name u; Type Local; NameOfSpace Hcurl;}
+ }
+ Equation {
+ Galerkin { [-1/mur[]*Dof{Curl u} , {Curl u}]; In All_domains; Jacobian JVol; Integration Int_1; }
+ Galerkin { [(omega0/cel)^2*epsilonr[]*Dof{u} , {u} ]; In All_domains; Jacobian JVol; Integration Int_1; }
+ Galerkin { [source[] , {u} ]; In All_domains; Jacobian JVol; Integration Int_1; }
+ }
+ }
+ EndIf
+ If (flag_study==RES_TMAT)
+ For pe In {1:p_max}
+ {Name VPWM_helmholtz_vector~{pe}; Type FemEquation;
+ Quantity {
+ { Name u; Type Local; NameOfSpace Hcurl;}
+ }
+ Equation {
+ Galerkin { [-1/mur[]*Dof{Curl u} , {Curl u}]; In All_domains; Jacobian JVol; Integration Int_1; }
+ Galerkin { [(omega0/cel)^2*epsilonr[]*Dof{u} , {u} ]; In All_domains; Jacobian JVol; Integration Int_1; }
+ Galerkin { [source_M~{pe}[] , {u} ]; In All_domains; Jacobian JVol; Integration Int_1; }
+ }
+ }
+ {Name VPWN_helmholtz_vector~{pe}; Type FemEquation;
+ Quantity {
+ { Name u; Type Local; NameOfSpace Hcurl;}
+ }
+ Equation {
+ Galerkin { [-1/mur[]*Dof{Curl u} , {Curl u}]; In All_domains; Jacobian JVol; Integration Int_1; }
+ Galerkin { [(omega0/cel)^2*epsilonr[]*Dof{u} , {u} ]; In All_domains; Jacobian JVol; Integration Int_1; }
+ Galerkin { [source_N~{pe}[] , {u} ]; In All_domains; Jacobian JVol; Integration Int_1; }
+ }
+ }
+ EndFor
+ EndIf
+ If (flag_study==RES_GREEN)
+ For ncomp In {0:2}
+ {Name GreenAll_helmholtz_vector~{ncomp}; Type FemEquation;
+ Quantity {
+ { Name u; Type Local; NameOfSpace Hcurl;}
+ }
+ Equation {
+ Galerkin { [-1/mur[]*Dof{Curl u} , {Curl u}]; In All_domains; Jacobian JVol; Integration Int_1; }
+ Galerkin { [(omega0/cel)^2*epsilonr[]*Dof{u} , {u} ]; In All_domains; Jacobian JVol; Integration Int_1; }
+ Galerkin { [sourceG~{ncomp}[] , {u} ]; In All_domains; Jacobian JVol; Integration Int_1; }
+ }
+ }
+ EndFor
+ EndIf
+}
+
+Resolution {
+ If (flag_study==RES_PW)
+ { Name res_PW_helmholtz_vector;
+ System {
+ { Name P; NameOfFormulation PW_helmholtz_vector; Type ComplexValue; Frequency Freq; }
+ }
+ Operation {
+ CreateDir[Str[myDir]];
+ Evaluate[Python[]{"scattering_init.py"}];
+ Generate[P];
+ Solve[P];
+ PostOperation[PW_postop];
+ Evaluate[Python[]{"scattering_post.py"}];
+ }
+ }
+ EndIf
+ If (flag_study==RES_TMAT)
+ { Name res_VPWall_helmholtz_vector;
+ System {
+ For pe In {1:p_max}
+ { Name M~{pe}; NameOfFormulation VPWM_helmholtz_vector~{pe}; Type ComplexValue; Frequency Freq; }
+ { Name N~{pe}; NameOfFormulation VPWN_helmholtz_vector~{pe}; Type ComplexValue; Frequency Freq; }
+ EndFor
+ }
+ Operation {
+ CreateDir[Str[myDir]];
+ Evaluate[Python[]{"scattering_init.py"}];
+ For pe In {1:p_max}
+ Generate[M~{pe}];
+ Solve[M~{pe}];
+ PostOperation[VPWM_postop~{pe}];
+ Generate[N~{pe}];
+ Solve[N~{pe}];
+ PostOperation[VPWN_postop~{pe}];
+ EndFor
+ Evaluate[Python[]{"scattering_post.py"}];
+ }
+ }
+ EndIf
+ If (flag_study==RES_GREEN)
+ { Name res_GreenAll_helmholtz_vector;
+ System {
+ For ncomp In {0:2}
+ { Name M~{ncomp}; NameOfFormulation GreenAll_helmholtz_vector~{ncomp}; Type ComplexValue; Frequency Freq; }
+ EndFor
+ }
+ Operation {
+ CreateDir[Str[myDir]];
+ Evaluate[Python[]{"scattering_init.py"}];
+ For ncomp In {0:2}
+ Generate[M~{ncomp}];
+ Solve[M~{ncomp}];
+ PostOperation[GreenAll_postop~{ncomp}];
+ EndFor
+ Evaluate[Python[]{"scattering_post.py"}];
+ }
+ }
+ EndIf
+ If (flag_study==RES_QNM)
+ { Name res_QNM_helmholtz_vector;
+ System{
+ { Name E; NameOfFormulation QNM_helmholtz_vector; Type ComplexValue; }
+ }
+ Operation{
+ CreateDir[Str[myDir]];
+ Evaluate[Python[]{"scattering_init.py"}];
+ GenerateSeparate[E];
+ EigenSolve[E,neig,eigtarget,0,EigFilter[]];
+ PostOperation[QNM_postop];
+ Evaluate[Python[]{"scattering_post.py"}];
+ }
+ }
+ EndIf
+
+}
+
+PostProcessing {
+ If (flag_study==RES_QNM)
+ { Name QNM_postpro; NameOfFormulation QNM_helmholtz_vector; NameOfSystem E;
+ Quantity {
+ { Name u ; Value { Local { [{u}]; In All_domains; Jacobian JVol; } } }
+ { Name EigenValuesReal; Value { Local{ [$EigenvalueReal]; In PrintPoint; Jacobian JVol; } } }
+ { Name EigenValuesImag; Value { Local{ [$EigenvalueImag]; In PrintPoint; Jacobian JVol; } } }
+ }
+ }
+ EndIf
+
+ If (flag_study==RES_PW)
+ { Name PW_postpro ; NameOfFormulation PW_helmholtz_vector;NameOfSystem P;
+ Quantity {
+ { Name E_tot ; Value { Local { [{u}+Einc[]]; In All_domains; Jacobian JVol; } } }
+ { Name E_scat ; Value { Local { [{u}]; In All_domains; Jacobian JVol; } } }
+ { Name E_scat_sph ; Value { Local { [Vector[
+ (X[]*CompX[{u}] + Y[]*CompY[{u}] + Z[]*CompZ[{u}])/r3D_sph[],
+ (X[]*Z[]*CompX[{u}] + Y[]*Z[]*CompY[{u}] - (X[]^2+Y[]^2)*CompZ[{u}])/(r3D_sph[]*r2D_sph[]),
+ (-Y[]*CompX[{u}] + X[]*CompY[{u}])/r2D_sph[]
+ ]];
+ In All_domains; Jacobian JVol; } } }
+ { Name H_scat ; Value { Local { [siwt*I[]/(mur[]*mu0*omega0)*{Curl u}] ; In All_domains; Jacobian JVol; } } }
+ { Name H_tot ; Value { Local { [ siwt*I[]/(mur[]*mu0*omega0)*{Curl u} +Hinc[]]; In All_domains; Jacobian JVol; } } }
+ { Name normalized_losses ; Value { Integral { [ 0.5*omega0*epsilon0*Fabs[Im[epsilonr_In[]]]*(SquNorm[{u}+Einc[]]) ] ; In Scat_In ; Integration Int_1 ; Jacobian JVol ; } } }
+ { Name Poy_dif ; Value { Local { [ 0.5*Re[Cross[{u} , Conj[ siwt*I[]/(mur[]*mu0*omega0)*{Curl u}]]] ]; In All_domains; Jacobian JVol; } } }
+ { Name Poy_tot ; Value { Local { [ 0.5*Re[Cross[{u}+Einc[] , Conj[Hinc[]+siwt*I[]/(mur[]*mu0*omega0)*{Curl u}]]] ]; In All_domains; Jacobian JVol; } } }
+ }
+ }
+ EndIf
+ If (flag_study==RES_TMAT)
+ For pe In {1:p_max}
+ { Name VPWM_postpro~{pe}; NameOfFormulation VPWM_helmholtz_vector~{pe};NameOfSystem M~{pe};
+ Quantity {
+ { Name E_scat ; Value { Local { [{u}]; In All_domains; Jacobian JVol; } } }
+ { Name E_scat_sph ; Value { Local { [Vector[
+ (X[]*CompX[{u}] + Y[]*CompY[{u}] + Z[]*CompZ[{u}])/r3D_sph[],
+ (X[]*Z[]*CompX[{u}] + Y[]*Z[]*CompY[{u}] - (X[]^2+Y[]^2)*CompZ[{u}])/(r3D_sph[]*r2D_sph[]),
+ (-Y[]*CompX[{u}] + X[]*CompY[{u}])/r2D_sph[]
+ ]];
+ In All_domains; Jacobian JVol; } } }
+ { Name H_scat ; Value { Local { [siwt*I[]/(mur[]*mu0*omega0)*{Curl u}]; In All_domains; Jacobian JVol; } } }
+ { Name Mnm_source~{pe} ; Value { Local { [ Mnm_source~{pe}[] ]; In All_domains; Jacobian JVol; } } }
+ }
+ }
+ { Name VPWN_postpro~{pe}; NameOfFormulation VPWN_helmholtz_vector~{pe};NameOfSystem N~{pe};
+ Quantity {
+ { Name E_scat ; Value { Local { [{u}]; In All_domains; Jacobian JVol; } } }
+ { Name E_scat_sph ; Value { Local { [Vector[
+ (X[]*CompX[{u}] + Y[]*CompY[{u}] + Z[]*CompZ[{u}])/r3D_sph[],
+ (X[]*Z[]*CompX[{u}] + Y[]*Z[]*CompY[{u}] - (X[]^2+Y[]^2)*CompZ[{u}])/(r3D_sph[]*r2D_sph[]),
+ (-Y[]*CompX[{u}] + X[]*CompY[{u}])/r2D_sph[]
+ ]];
+ In All_domains; Jacobian JVol; } } }
+ { Name H_scat ; Value { Local { [siwt*I[]/(mur[]*mu0*omega0)*{Curl u}]; In All_domains; Jacobian JVol; } } }
+ { Name Nnm_source~{pe} ; Value { Local { [ Nnm_source~{pe}[] ]; In All_domains; Jacobian JVol; } } }
+ }
+ }
+ EndFor
+ EndIf
+ If (flag_study==RES_GREEN)
+ For ncomp In {0:2}
+ { Name GreenAll_postpro~{ncomp}; NameOfFormulation GreenAll_helmholtz_vector~{ncomp};NameOfSystem M~{ncomp};
+ Quantity {
+ { Name E_tot~{ncomp} ; Value { Local {[{u}+Einc_Green~{ncomp}[]]; In All_domains; Jacobian JVol; } } }
+ { Name H_tot~{ncomp} ; Value { Local {[siwt*I[]/(mur[]*mu0*omega0)*{Curl u}+Hinc_Green~{ncomp}[]]; In All_domains; Jacobian JVol; } } }
+ { Name Einc_Green~{ncomp}; Value { Local {[ Einc_Green~{ncomp}[] ]; In All_domains; Jacobian JVol; } } }
+ { Name Hinc_Green~{ncomp}; Value { Local {[ Hinc_Green~{ncomp}[] ]; In All_domains; Jacobian JVol; } } }
+ { Name E_scat~{ncomp} ; Value { Local { [{u}]; In All_domains; Jacobian JVol; } } }
+ { Name E_tot_im~{ncomp} ; Value { Local {[Im[{u}+Einc_Green~{ncomp}[]]]; In All_domains; Jacobian JVol; } } }
+ { Name Einc_Green_im~{ncomp}; Value { Local {[Im[ Einc_Green~{ncomp}[] ]]; In All_domains; Jacobian JVol; } } }
+ // { Name Einc_Green_im~{ncomp}; Value { Local {[test[]]; In All_domains; Jacobian JVol; } } }
+ { Name E_scat_im~{ncomp} ; Value { Local {[Im[{u} ]]; In All_domains; Jacobian JVol; } } }
+ { Name E_tot_sph~{ncomp}; Value { Local { [Vector[
+ (X[]*CompX[{u}+Einc_Green~{ncomp}[]]
+ + Y[]*CompY[{u}+Einc_Green~{ncomp}[]]
+ + Z[]*CompZ[{u}+Einc_Green~{ncomp}[]])/r3D_sph[],
+ (X[]*Z[]*CompX[{u}+Einc_Green~{ncomp}[]]
+ + Y[]*Z[]*CompY[{u}+Einc_Green~{ncomp}[]]
+ - (X[]^2+Y[]^2)*CompZ[{u}+Einc_Green~{ncomp}[]])
+ /(r3D_sph[]*r2D_sph[]),
+ (-Y[]*CompX[{u}+Einc_Green~{ncomp}[]]
+ + X[]*CompY[{u}+Einc_Green~{ncomp}[]])
+ /r2D_sph[] ]];
+ In All_domains; Jacobian JVol; } } }
+ { Name E_scat_sph~{ncomp}; Value { Local { [Vector[
+ (X[]*CompX[{u}] + Y[]*CompY[{u}] + Z[]*CompZ[{u}])/r3D_sph[],
+ (X[]*Z[]*CompX[{u}] + Y[]*Z[]*CompY[{u}]
+ - (X[]^2+Y[]^2)*CompZ[{u}])/(r3D_sph[]*r2D_sph[]),
+ (-Y[]*CompX[{u}] + X[]*CompY[{u}])/r2D_sph[]
+ ]];
+ In All_domains; Jacobian JVol; } } }
+ { Name H_scat~{ncomp} ; Value { Local { [siwt*I[]/(mur[]*mu0*omega0)*{Curl u}]; In All_domains; Jacobian JVol; } } }
+ }
+ }
+ EndFor
+ EndIf
+}
+
+
+
+PostOperation {
+ If (flag_study==RES_QNM)
+ { Name QNM_postop; NameOfPostProcessing QNM_postpro ;
+ Operation {
+ Print [EigenValuesReal, OnElementsOf PrintPoint , Format TimeTable, File StrCat[myDir,"EigenValuesReal.txt"]];
+ Print [EigenValuesImag, OnElementsOf PrintPoint , Format TimeTable, File StrCat[myDir,"EigenValuesImag.txt"]];
+ Print [ u , OnElementsOf All_domains, File StrCat[myDir,"eigenVectors.pos"], EigenvalueLegend];
+ }
+ }
+ EndIf
+ If (flag_study==RES_PW)
+ {Name PW_postop; NameOfPostProcessing PW_postpro ;
+ Operation {
+ Print [ E_tot , OnElementsOf Domain , File StrCat[myDir,"E_tot.pos"]];
+ Print [ E_scat , OnElementsOf All_domains, File StrCat[myDir,"E_scat.pos"]];
+ // Print [ E_scat_sph , OnElementsOf All_domains, File StrCat[myDir,"E_scat_sph.pos"]];
+ Print [ H_scat , OnElementsOf Domain , File StrCat[myDir,"H_scat.pos"]];
+ Print [ H_tot , OnElementsOf All_domains, File StrCat[myDir,"H_tot.pos"]];
+ Print [ Poy_dif , OnElementsOf All_domains, File StrCat[myDir,"Poy_dif.pos"]];
+ Print [ Poy_tot , OnElementsOf Domain , File StrCat[myDir,"Poy_tot.pos"]];
+ For kr In {0:nb_cuts-1}
+ Print [ E_scat_sph , OnGrid
+ {r_sph~{kr}*Sin[$B]*Cos[$C], r_sph~{kr}*Sin[$B]*Sin[$C], r_sph~{kr}*Cos[$B]}
+ {r_sph~{kr}, {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_theta-1.0)},
+ {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_phi-1.0)} },
+ File StrCat[myDir,StrCat["E_scat_onsphere_sph_PW",Sprintf["_r%g.dat",kr]]], Format Table];
+ EndFor
+ }
+ }
+ EndIf
+ If (flag_study==RES_TMAT)
+ For pe In {1:p_max}
+ {Name VPWM_postop~{pe}; NameOfPostProcessing VPWM_postpro~{pe} ;
+ Operation {
+ For kr In {0:nb_cuts-1}
+ Print [ E_scat_sph , OnGrid
+ {r_sph~{kr}*Sin[$B]*Cos[$C], r_sph~{kr}*Sin[$B]*Sin[$C], r_sph~{kr}*Cos[$B]}
+ {r_sph~{kr}, {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_theta-1.0)},
+ {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_phi-1.0)} },
+ File StrCat[myDir,StrCat[StrCat["E_scat_onsphere_sph_M",Sprintf["%g",pe]],Sprintf["_r%g.dat",kr]]], Format Table];
+ EndFor
+ Print [ E_scat , OnGrid
+ {r_sph~{kr}*Sin[$B]*Cos[$C], r_sph~{kr}*Sin[$B]*Sin[$C], r_sph~{kr}*Cos[$B]}
+ {r_sph~{kr}, {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_plot_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_plot_theta-1.0)},
+ {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_plot_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_plot_phi-1.0)} },
+ File StrCat[myDir,StrCat[StrCat["E_scat_onsphere_cart_M",Sprintf["%g",pe]],Sprintf["_r%g.pos",kr]]],
+ Name StrCat["E_scat_onsphere_cart_M",Sprintf["%g",pe]]];
+ }
+ }
+ {Name VPWN_postop~{pe}; NameOfPostProcessing VPWN_postpro~{pe} ;
+ Operation {
+ For kr In {0:nb_cuts-1}
+ Print [ E_scat_sph , OnGrid
+ {r_sph~{kr}*Sin[$B]*Cos[$C], r_sph~{kr}*Sin[$B]*Sin[$C], r_sph~{kr}*Cos[$B]}
+ {r_sph~{kr}, {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_theta-1.0)},
+ {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_phi-1.0)} },
+ File StrCat[myDir,StrCat[StrCat["E_scat_onsphere_sph_N",Sprintf["%g",pe]],Sprintf["_r%g.dat",kr]]], Format Table];
+ EndFor
+ Print [ E_scat , OnGrid
+ {r_sph~{kr}*Sin[$B]*Cos[$C], r_sph~{kr}*Sin[$B]*Sin[$C], r_sph~{kr}*Cos[$B]}
+ {r_sph~{kr}, {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_plot_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_plot_theta-1.0)},
+ {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_plot_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_plot_phi-1.0)} },
+ File StrCat[myDir,StrCat[StrCat["E_scat_onsphere_cart_N",Sprintf["%g",pe]],Sprintf["_r%g.pos",kr]]],
+ Name StrCat["E_scat_onsphere_cart_N",Sprintf["%g",pe]]];
+ }
+ }
+ EndFor
+ EndIf
+ If (flag_study==RES_GREEN)
+ For ncomp In {0:2}
+ {Name GreenAll_postop~{ncomp}; NameOfPostProcessing GreenAll_postpro~{ncomp} ;
+ Operation {
+ Print [ E_tot_sph~{ncomp} , OnGrid
+ {r_pml_in *Sin[$B]*Cos[$C], r_pml_in*Sin[$B]*Sin[$C], r_pml_in*Cos[$B]}
+ {r_pml_in , {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_theta-1.0)},
+ {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_phi-1.0)} },
+ File StrCat[myDir,StrCat["E_tot_onsphere_sph_G",Sprintf["%g.dat",ncomp]]] , Format Table];
+ Print [ E_tot~{ncomp} , OnGrid
+ {r_pml_in *Sin[$B]*Cos[$C], r_pml_in*Sin[$B]*Sin[$C], r_pml_in*Cos[$B]}
+ {r_pml_in , {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_plot_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_plot_theta-1.0)},
+ {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_plot_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_plot_phi-1.0)} },
+ File StrCat[myDir,StrCat["E_tot_onsphere_cart_G",Sprintf["%g.pos",ncomp]]]];
+ Print [ E_scat_im~{ncomp} , OnPoint {x_p,y_p,z_p},
+ File StrCat[myDir,StrCat["E_scat_im_onpt_G",Sprintf["%g.dat",ncomp]]] , Format Table];
+ // Print [ Einc_Green_im~{ncomp} , OnPoint {x_p+0.1*nm,y_p+0.1*nm,z_p+0.1*nm},
+ // File StrCat[myDir,StrCat["E_inc_im_onpt_G",Sprintf["%g.dat",ncomp]]] , Format Table];
+ // If(flag_FF==1)
+ // Print [ E_tot~{ncomp}, OnElementsOf Scat_Out , File StrCat[myDir,StrCat["E_tot_",Sprintf["%g.pos",ncomp]]]];
+ // Print [ H_tot~{ncomp}, OnElementsOf Scat_Out , File StrCat[myDir,StrCat["H_tot_",Sprintf["%g.pos",ncomp]]]];
+ // Print [ E_tot~{ncomp} , OnGrid
+ // {r_pml_in *Sin[$B]*Cos[$C], r_pml_in*Sin[$B]*Sin[$C], r_pml_in*Cos[$B]}
+ // {r_pml_in , {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_theta-1.0)},
+ // {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_phi-1.0)} },
+ // File StrCat[myDir,StrCat["E_tot_onsphere_G",Sprintf["%g.pos",ncomp]]]];
+ // Print [ H_tot~{ncomp} , OnGrid
+ // {r_pml_in *Sin[$B]*Cos[$C], r_pml_in*Sin[$B]*Sin[$C], r_pml_in*Cos[$B]}
+ // {r_pml_in , {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_theta-1.0)},
+ // {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_phi-1.0)} },
+ // File StrCat[myDir,StrCat["H_tot_onsphere_G",Sprintf["%g.pos",ncomp]]]];
+ //
+ // Echo[
+ // Str["Plugin(NearToFarField).Wavenumber = 2*Pi/lambda0;",
+ // "Plugin(NearToFarField).PhiStart = 0;",
+ // "Plugin(NearToFarField).PhiEnd = 2*Pi;",
+ // "Plugin(NearToFarField).NumPointsPhi = npts_phi;",
+ // "Plugin(NearToFarField).ThetaStart = 0;",
+ // "Plugin(NearToFarField).ThetaEnd = Pi;",
+ // "Plugin(NearToFarField).NumPointsTheta = npts_theta;",
+ // "Plugin(NearToFarField).EView = PostProcessing.NbViews-2;",
+ // "Plugin(NearToFarField).HView = PostProcessing.NbViews-1;",
+ // "Plugin(NearToFarField).Normalize = 0;",
+ // "Plugin(NearToFarField).dB = 0;",
+ // "Plugin(NearToFarField).NegativeTime = 0;",
+ // "Plugin(NearToFarField).RFar = r_pml_in;",
+ // "Plugin(NearToFarField).Run;"], File StrCat[myDir,"tmp1.geo"]] ;
+ // EndIf
+ }
+ }
+ EndFor
+ EndIf
+}
+
+If (flag_study==RES_QNM)
+ DefineConstant[
+ R_ = {"res_QNM_helmholtz_vector", Name "GetDP/1ResolutionChoices", Visible 1},
+ C_ = {"-solve -pos -petsc_prealloc 500 -ksp_type preonly -pc_type lu -pc_factor_mat_solver_package mumps", Name "GetDP/9ComputeCommand", Visible 1},
+ P_ = {"", Name "GetDP/2PostOperationChoices", Visible 0}];
+EndIf
+
+If (flag_study==RES_TMAT)
+ DefineConstant[
+ R_ = {"res_VPWall_helmholtz_vector", Name "GetDP/1ResolutionChoices", Visible 1},
+ C_ = {"-solve -pos -petsc_prealloc 500 -ksp_type preonly -pc_type lu -pc_factor_mat_solver_package mumps", Name "GetDP/9ComputeCommand", Visible 1},
+ P_ = {"", Name "GetDP/2PostOperationChoices", Visible 0}];
+EndIf
+If (flag_study==RES_PW)
+ DefineConstant[
+ R_ = {"res_PW_helmholtz_vector", Name "GetDP/1ResolutionChoices", Visible 1},
+ C_ = {"-solve -pos -petsc_prealloc 500 -ksp_type preonly -pc_type lu -pc_factor_mat_solver_package mumps", Name "GetDP/9ComputeCommand", Visible 1},
+ P_ = {"", Name "GetDP/2PostOperationChoices", Visible 0}];
+EndIf
+If (flag_study==RES_GREEN)
+ DefineConstant[
+ R_ = {"res_GreenAll_helmholtz_vector", Name "GetDP/1ResolutionChoices", Visible 1},
+ C_ = {"-solve -pos -petsc_prealloc 500 -ksp_type preonly -pc_type lu -pc_factor_mat_solver_package mumps", Name "GetDP/9ComputeCommand", Visible 1},
+ P_ = {"", Name "GetDP/2PostOperationChoices", Visible 0}];
+EndIf
\ No newline at end of file
diff --git a/ElectromagneticScattering/scattering_data.geo b/ElectromagneticScattering/scattering_data.geo
new file mode 100644
index 0000000000000000000000000000000000000000..91e3c08c3ced34e22ad336208c257be16360a813
--- /dev/null
+++ b/ElectromagneticScattering/scattering_data.geo
@@ -0,0 +1,215 @@
+///////////////////////////////
+// Author : Guillaume Demesy //
+// scattering_data.geo //
+///////////////////////////////
+
+nm = 1.;
+epsilon0 = 8.854187817e-3*nm;
+mu0 = 400.*Pi*nm;
+cel = 1.0/(Sqrt[epsilon0 * mu0]);
+deg2rad = Pi/180.;
+
+pp0 = "1Geometry/0";
+pp1 = "2Study Type/0";
+pp2 = "3Electromagnetic parameters/0";
+pp3 = "4Mesh size and PMLs parameters/0";
+pp4 = "5Postpro options/0";
+pp5 = "6Post plot options/0";
+close_menu = 0;
+colorro = "LightGrey";
+colorppOK = "Ivory";
+colorppWA = "LightSalmon1";
+colorppNO = "LightSalmon4";
+
+ELL = 0;
+PARALL = 1;
+CYL = 2;
+CONE = 3;
+TOR = 4;
+
+RES_PW = 0;
+RES_TMAT = 1;
+RES_GREEN = 2;
+RES_QNM = 3;
+
+DefineConstant[
+ flag_shape = {ELL , Name StrCat[pp0, "0Scatterer shape"],
+ Choices {ELL="ellispoid",PARALL="parallelepiped",CYL="cylinder",CONE="cone",TOR="split torus"},Closed 0,GmshOption "Reset", Autocheck 0}];
+
+If (flag_shape==ELL)
+ DefineConstant[
+ ell_rx = {73.45 , Name StrCat[pp0 , "1ellipsoid X-radius [nm]"], Highlight Str[colorppOK] , Closed 0},
+ ell_ry = {73.45 , Name StrCat[pp0 , "2ellipsoid Y-radius [nm]"], Highlight Str[colorppOK] , Closed 0},
+ ell_rz = {73.45 , Name StrCat[pp0 , "3ellipsoid Z-radius [nm]"], Highlight Str[colorppOK] , Closed 0}
+ ];
+ // FIXME gmsh Max?
+ If (ell_rx>ell_ry)
+ rbb_tmp=ell_rx;
+ Else rbb_tmp=ell_ry;
+ EndIf
+ If (ell_rz>rbb_tmp)
+ rbb=ell_rz;
+ Else rbb=rbb_tmp;
+ EndIf
+EndIf
+If (flag_shape==PARALL)
+ DefineConstant[
+ par_ax = {300 , Name StrCat[pp0 , "1cube X-edge size [nm]"], Highlight Str[colorppOK] , Closed 0},
+ par_ay = {100 , Name StrCat[pp0 , "2cube Y-edge size [nm]"], Highlight Str[colorppOK] , Closed 0},
+ par_az = {600 , Name StrCat[pp0 , "3cube Z-edge size [nm]"], Highlight Str[colorppOK] , Closed 0}];
+ rbb = 0.5*Sqrt[par_ax^2+par_ay^2+par_az^2];
+EndIf
+If (flag_shape==CYL)
+ DefineConstant[
+ cyl_rx = {300 , Name StrCat[pp0 , "1cylinder X-radius [nm]"], Highlight Str[colorppOK] , Closed 0},
+ cyl_ry = {300 , Name StrCat[pp0 , "2cylinder Y-radius [nm]"], Highlight Str[colorppOK] , Closed 0},
+ cyl_h = {300 , Name StrCat[pp0 , "3cylinder height [nm]"] , Highlight Str[colorppOK] , Closed 0}];
+ // FIXME gmsh Max?
+ If (cyl_rx>cyl_ry)
+ rbb_tmp=cyl_rx;
+ Else rbb_tmp=cyl_ry;
+ EndIf
+ rbb = 0.5*Sqrt[rbb_tmp^2+cyl_h^2];
+EndIf
+If (flag_shape==CONE)
+ DefineConstant[
+ cone_rx = {300 , Name StrCat[pp0 , "1cone basis X-radius [nm]"], Highlight Str[colorppOK] , Closed 0},
+ cone_ry = {300 , Name StrCat[pp0 , "1cone basis Y-radius [nm]"], Highlight Str[colorppOK] , Closed 0},
+ cone_h = {300 , Name StrCat[pp0 , "2cone height [nm]"] , Highlight Str[colorppOK] , Closed 0}];
+ // FIXME gmsh Max?
+ If (cone_rx>cone_ry)
+ rbb_tmp=cone_rx;
+ Else rbb_tmp=cone_ry;
+ EndIf
+ If (2.0*cone_h/3.0>Sqrt[rbb_tmp^2+cone_h^2/9.0])
+ rbb=2.0*cone_h/3.0;
+ Else rbb=Sqrt[rbb_tmp^2+cone_h^2/9.0];
+ EndIf
+EndIf
+If (flag_shape==TOR)
+ DefineConstant[
+ tor_r1 = { 300 , Name StrCat[pp0 , "1torus radius 1 [nm]"], Highlight Str[colorppOK] , Closed 0},
+ tor_r2x = { 100 , Name StrCat[pp0 , "2torus radius 2x [nm]"], Highlight Str[colorppOK] , Closed 0},
+ tor_r2z = { 50 , Name StrCat[pp0 , "3torus radius 2z [nm]"], Highlight Str[colorppOK] , Closed 0},
+ tor_angle = { 340 , Name StrCat[pp0 , "4torus angle [deg]"] , Highlight Str[colorppOK] , Closed 0, Min 5, Max 355}];
+ rbb = tor_r1+tor_r2x;
+EndIf
+DefineConstant[
+ rot_theta = {0 , Name StrCat[pp0 , "5rotate scatterer (polar) [deg]"] , Highlight Str[colorppOK] , Closed 0, Min 0, Max 180},
+ rot_phi = {0 , Name StrCat[pp0 , "6rotate scatterer (azimut) [deg]"], Highlight Str[colorppOK] , Closed 0, Min 0, Max 360}];
+
+DefineConstant[
+ flag_study = { RES_TMAT , Name StrCat[pp1, "0Study type"],
+ Choices {RES_PW="Plane Wave Response",
+ RES_TMAT="T-Matrix",
+ RES_GREEN="Green's Tensor",
+ RES_QNM="Quasi-Normal Modes"},Closed 0,GmshOption "Reset", Autocheck 0}];
+
+DefineConstant[
+ epsr_In_re = { 9., Name StrCat[pp2 , "0scatterer permittivity (real) []"] , Highlight Str[colorppOK] , Closed 0},
+ epsr_In_im = { 0., Name StrCat[pp2 , "1scatterer permittivity (imag) []"] , Highlight Str[colorppOK] , Closed 0},
+ epsr_Out_re = { 1., Name StrCat[pp2 , "2background permittivity (real) []"], Highlight Str[colorppOK] , Closed 0},
+ epsr_Out_im = { 0., Name StrCat[pp2 , "3background permittivity (imag) []"], Highlight Str[colorppOK] , Closed 0}
+];
+
+If (flag_study!=RES_QNM)
+ DefineConstant[
+ lambda0 = {587.6 , Name StrCat[pp2 , "4wavelength [nm]"] , Highlight Str[colorppOK] , Closed 0,GmshOption "Reset", Autocheck 0}];
+Else
+ DefineConstant[
+ lambda0 = {1000 , Name StrCat[pp2 , "4wavelength target [nm]"] , Highlight Str[colorppOK] , Closed 0,GmshOption "Reset", Autocheck 0},
+ neig = {30 , Name StrCat[pp2 , "5number of eigenvalues []"] , Highlight Str[colorppOK] , Closed 0,GmshOption "Reset", Autocheck 0}
+ ];
+EndIf
+lambda_bg = lambda0/Sqrt[epsr_Out_re];
+
+If (flag_study==RES_PW)
+ DefineConstant[
+ theta0 = {0 , Name StrCat[pp2 , "5plane wave theta [deg]"], Highlight Str[colorppOK] , Closed 0, Min 0, Max 180},
+ phi0 = {0 , Name StrCat[pp2 , "6plane wave phi [deg]"] , Highlight Str[colorppOK] , Closed 0, Min 0, Max 360},
+ psi0 = {0 , Name StrCat[pp2 , "7plane wave psi [deg]"] , Highlight Str[colorppOK] , Closed 0, Min 0, Max 180}];
+ theta0 = theta0 *deg2rad;
+ phi0 = phi0 * deg2rad;
+ psi0 = psi0 * deg2rad;
+EndIf
+If (flag_study==RES_GREEN)
+ DefineConstant[
+ x_p = {0 , Name StrCat[pp2 , "5Green X-point [nm]"] , Highlight Str[colorppOK] , Closed 0},
+ y_p = {0 , Name StrCat[pp2 , "6Green Y-point [nm]"] , Highlight Str[colorppOK] , Closed 0},
+ z_p = {-rbb*1.1 , Name StrCat[pp2 , "7Green Z-point [nm]"] , Highlight Str[colorppOK] , Closed 0}];
+ x_p=x_p*nm;
+ y_p=y_p*nm;
+ z_p=z_p*nm;
+EndIf
+
+DefineConstant[
+ n_max = {1 , Name StrCat[pp2 , "8n_max integer"], Highlight Str[colorppOK] , Closed 0, Min 0, Max 5},
+ siwt = {0 , Name StrCat[pp2 , "9Time sign e^(+|-iwt)"], Choices {0="e^(-iwt)",2="e^(+iwt)"} , Closed 0}
+];
+DefineConstant[
+ flag_cartpml = {0 , Name StrCat[pp3 , "0PML type"] , Choices {0="spherical",1="cartesian"}, Closed 0},
+ pml_size = {lambda_bg/1.5, Name StrCat[pp3 , "1PML thickness [nm]"] , Highlight Str[colorppWA] , Closed 0, Min (lambda_bg/10), Max (3*lambda_bg)},
+ paramaille = {4. , Name StrCat[pp3 , "2mesh size"], Highlight Str[colorppWA] , Closed 0},
+ refine_scat = {2. , Name StrCat[pp3 , "3scatterer mesh refinement"], Highlight Str[colorppWA] , Closed 0}
+ is_FEM_o2 = {1 , Name StrCat[pp3 , "4Interpolation order "] , Choices {0="order 1",1="order 2"}, Closed 0}
+];
+DefineConstant[
+ space2pml = {lambda_bg/10.*4 , Name StrCat[pp4 , "0space around scatterer [nm]"] , Highlight Str[colorppNO] , Closed 1,Min (lambda_bg/10.), Max (3*lambda_bg)},
+ dist_rcut = {lambda_bg/10. , Name StrCat[pp4 , "1min dist from object & PML"] , Highlight Str[colorppNO] , Closed 1},
+ nb_cuts = {3 , Name StrCat[pp4 , "2number of cuts [integer]"] , Highlight Str[colorppNO] , Closed 1},
+ npts_theta = {50 , Name StrCat[pp4 , "3polar sampling [integer]"], Highlight Str[colorppNO] , Closed 0,Min 50, Max 300},
+ npts_phi = {100 , Name StrCat[pp4 , "4azimuthal sampling [integer]"], Highlight Str[colorppNO] , Closed 0,Min 100, Max 600}
+];
+DefineConstant[
+ flag_plotcuts = {0, Choices{0,1}, Name StrCat[pp5, "Plot radial cuts?"]},
+ flag_FF = {1, Choices{0,1}, Name StrCat[pp5, "Plot far field?"]}
+];
+
+
+
+// FIXME conditional for normalization
+If (flag_shape==ELL)
+ ell_rx = ell_rx*nm;
+ ell_ry = ell_ry*nm;
+ ell_rz = ell_rz*nm;
+EndIf
+If (flag_shape==PARALL)
+ par_ax = par_ax*nm;
+ par_ay = par_ay*nm;
+ par_az = par_az*nm;
+EndIf
+If (flag_shape==CYL)
+ cyl_rx = cyl_rx*nm;
+ cyl_ry = cyl_ry*nm;
+ cyl_h = cyl_h *nm;
+EndIf
+If (flag_shape==CONE)
+ cone_rx = cone_rx*nm;
+ cone_ry = cone_ry*nm;
+ cone_h = cone_h*nm;
+EndIf
+If (flag_shape==TOR)
+ tor_r1 = tor_r1*nm;
+ tor_r2x = tor_r2x*nm;
+ tor_r2z = tor_r2z*nm;
+EndIf
+
+lambda0 = lambda0*nm;
+lambda_bg = lambda_bg*nm;
+pml_size = pml_size*nm;
+dist_rcut = dist_rcut*nm;
+space2pml = space2pml*nm;
+rbb = rbb*nm;
+
+r_pml_in = space2pml+rbb;
+r_pml_out = space2pml+rbb+pml_size;
+r_sph_min = rbb+dist_rcut;
+r_sph_max = space2pml+rbb-dist_rcut;
+
+sph_scan = 0.00001;
+
+npts_plot_theta = 25;
+npts_plot_phi = 50;
+
+p_max = n_max*n_max+2*n_max;
+siwt=siwt-1;
diff --git a/ElectromagneticScattering/scattering_init.py b/ElectromagneticScattering/scattering_init.py
new file mode 100644
index 0000000000000000000000000000000000000000..2f7621c6a6b329b909510f8b92f14a5f33529dcf
--- /dev/null
+++ b/ElectromagneticScattering/scattering_init.py
@@ -0,0 +1,258 @@
+# -*- coding: utf-8 -*-
+
+"""
+///////////////////////////////
+// Author : Guillaume Demesy //
+// scattering_init.py //
+///////////////////////////////
+"""
+
+import sys,os
+import numpy as np
+from scipy.special import jv, yv, hankel1, hankel2
+sys.path.append(os.getcwd())
+from matplotlib import cm
+import pylab as pl
+from scattering_tmp import *
+np.set_printoptions(precision=2)
+pi = np.pi
+
+ps = np.linspace(1,p_max,p_max)
+r_sphs = np.linspace(r_sph_min,r_sph_max,nb_cuts)
+k_Out = 2*pi*np.sqrt(epsr_Out_re)/lambda0
+epsr_In = epsr_In_re+1j*epsr_In_im
+epsr_Out = epsr_Out_re+1j*epsr_Out_im
+
+phi_range = np.linspace(sph_scan,2*pi-sph_scan,npts_phi)
+theta_range = np.linspace(sph_scan, pi-sph_scan,npts_theta)
+[phi_sph,theta_sph]=np.meshgrid(phi_range,theta_range)
+cos_theta = np.cos(theta_range)
+sin_theta = np.sin(phi_range)
+
+def field_VSH_expansion(post_filename):
+ m_max = n_max
+ p_max = n_max**2 +2*n_max
+ FF_Xnm_t = np.zeros((npts_theta,npts_phi,p_max),dtype=complex)
+ FF_Xnm_p = np.zeros((npts_theta,npts_phi,p_max),dtype=complex)
+ FF_erCrossXnm_t = np.zeros((npts_theta,npts_phi,p_max),dtype=complex)
+ FF_erCrossXnm_p = np.zeros((npts_theta,npts_phi,p_max),dtype=complex)
+ #####################################
+ ##### sn, pn ,un
+ ##### Brian's recurrence relations
+ B_Pnpms = np.zeros((npts_theta,m_max+1,n_max+1))
+ B_unpms = np.zeros((npts_theta,m_max+1,n_max+1))
+ B_snpms = np.zeros((npts_theta,m_max+1,n_max+1))
+ ### Init
+ B_Pnpms[:,0,0] = np.sqrt(1./(4.*pi))
+ B_unpms[:,0,0] = 0.
+ B_unpms[:,1,1] = -0.25*np.sqrt(3./pi)
+ for k in range(npts_theta):
+ u=cos_theta[k]
+ for n in range(1,n_max+1):
+ B_Pnpms[k,n,n] = -np.sqrt((2.*float(n)+1.)/(2.*float(n))) * np.sqrt(1.-u**2) * B_Pnpms[k,n-1,n-1]
+ B_Pnpms[k,n-1,n] = u*np.sqrt(2.*float(n)+1.) * B_Pnpms[k,n-1,n-1]
+ for n in range(2,n_max+1):
+ B_unpms[k,n,n] = -np.sqrt( (float(n)*(2.*float(n)+1.)) / (2.*(float(n)+1.)*(float(n)-1.)) ) * np.sqrt(1.-u**2) * B_unpms[k,n-1,n-1]
+ B_unpms[k,n-1,n] = np.sqrt( (2.*float(n)+1.)*(float(n)-1.)/(float(n)+1.) ) * u * B_unpms[k,n-1,n-1]
+ for n in range(2,n_max+1):
+ for m in range(n-2+1):
+ B_Pnpms[k,m,n] = np.sqrt((4.*(float(n))**2-1.)/((float(n))**2-(float(m))**2)) * u * B_Pnpms[k,m,n-1] \
+ - np.sqrt( ( (2.*float(n)+1.)*((float(n)-1.)**2-(float(m))**2) ) \
+ / ( (2.*float(n)-3.)*((float(n))**2-(float(m))**2) ) )*B_Pnpms[k,m,n-2]
+
+ B_unpms[k,m,n] = np.sqrt( ((4.*(float(n))**2-1.)*(float(n)-1.))/((float(n)**2-float(m)**2)*(float(n)+1.)) ) * u * B_unpms[k,m,n-1] \
+ - np.sqrt( ((2.*float(n)+1.) * (float(n)-1.) * (float(n)-2.) * (float(n-m)-1.) *(float(n+m)-1.)) \
+ / ((2.*float(n)-3.) * (float(n)**2-float(m)**2)*float(n)*(float(n)+1.)))*B_unpms[k,m,n-2]
+ for n in range(0,n_max+1):
+ m=0
+ B_snpms[k,m,n] = 1./float(m+1)*np.sqrt((float(n+m)+1.)*(float(n-m))) * np.sqrt(1.-u**2) *\
+ B_unpms[k,m+1,n] + u*B_unpms[k,m,n]
+ for n in range(1,n_max+1):
+ for m in range(1,n+1):
+ B_snpms[k,m,n] = float(n)/float(m) * u * B_unpms[k,m,n] - float(m+n)/float(m) * \
+ np.sqrt( ( (2.*float(n)+1.)*(float(n)-float(m))*(float(n)-1.) ) / \
+ ( (2.*float(n)-1.)*(float(n)+float(m))*(float(n)+1.) ) )*B_unpms[k,m,n-1]
+ B_Pnmms = np.zeros_like(B_Pnpms)
+ B_unmms = np.zeros_like(B_unpms)
+ B_snmms = np.zeros_like(B_snpms)
+ for m in range(m_max+1):
+ B_Pnmms[:,m,:] = (-1.0)**m * B_Pnpms[:,m,:]
+ B_unmms[:,m,:] = (-1.0)**(m+1) * B_unpms[:,m,:]
+ B_snmms[:,m,:] = (-1.0)**(m) * B_snpms[:,m,:]
+ B_Pnmms=B_Pnmms[:,::-1,:]
+ B_unmms=B_unmms[:,::-1,:]
+ B_snmms=B_snmms[:,::-1,:]
+ B_Pnms = np.concatenate((B_Pnmms,B_Pnpms[:,1::,:]),axis=1)
+ B_unms = np.concatenate((B_unmms,B_unpms[:,1::,:]),axis=1)
+ B_snms = np.concatenate((B_snmms,B_snpms[:,1::,:]),axis=1)
+
+ #####################################
+ ##### sn, pn ,un
+ ##### Brian's recurrence relations
+ m_max = n_max
+ p_max = n_max*n_max+2*n_max
+ aM_nm = np.zeros(p_max,dtype=complex)
+ bN_nm = np.zeros(p_max,dtype=complex)
+ fenm_Y = np.zeros(p_max,dtype=complex)
+ fenm_Z = np.zeros(p_max,dtype=complex)
+ fhnm_X = np.zeros(p_max,dtype=complex)
+
+ B_Pnpms = np.zeros((npts_theta,m_max+1,n_max+1))
+ B_unpms = np.zeros((npts_theta,m_max+1,n_max+1))
+ B_snpms = np.zeros((npts_theta,m_max+1,n_max+1))
+ ### Init
+ B_Pnpms[:,0,0] = np.sqrt(1./(4.*pi))
+ B_unpms[:,0,0] = 0.
+ B_unpms[:,1,1] = -0.25*np.sqrt(3./pi)
+ for k in range(npts_theta):
+ u=cos_theta[k]
+ for n in range(1,n_max+1):
+ B_Pnpms[k,n,n] = -np.sqrt((2.*float(n)+1.)/(2.*float(n))) * np.sqrt(1.-u**2) * B_Pnpms[k,n-1,n-1]
+ B_Pnpms[k,n-1,n] = u*np.sqrt(2.*float(n)+1.) * B_Pnpms[k,n-1,n-1]
+ for n in range(2,n_max+1):
+ B_unpms[k,n,n] = -np.sqrt( (float(n)*(2.*float(n)+1.)) / (2.*(float(n)+1.)*(float(n)-1.)) ) * np.sqrt(1.-u**2) * B_unpms[k,n-1,n-1]
+ B_unpms[k,n-1,n] = np.sqrt( (2.*float(n)+1.)*(float(n)-1.)/(float(n)+1.) ) * u * B_unpms[k,n-1,n-1]
+ for n in range(2,n_max+1):
+ for m in range(n-2+1):
+ B_Pnpms[k,m,n] = np.sqrt((4.*(float(n))**2-1.)/((float(n))**2-(float(m))**2)) * u * B_Pnpms[k,m,n-1] \
+ - np.sqrt( ( (2.*float(n)+1.)*((float(n)-1.)**2-(float(m))**2) ) \
+ / ( (2.*float(n)-3.)*((float(n))**2-(float(m))**2) ) )*B_Pnpms[k,m,n-2]
+
+ B_unpms[k,m,n] = np.sqrt( ((4.*(float(n))**2-1.)*(float(n)-1.))/((float(n)**2-float(m)**2)*(float(n)+1.)) ) * u * B_unpms[k,m,n-1] \
+ - np.sqrt( ((2.*float(n)+1.) * (float(n)-1.) * (float(n)-2.) * (float(n-m)-1.) *(float(n+m)-1.)) \
+ / ((2.*float(n)-3.) * (float(n)**2-float(m)**2)*float(n)*(float(n)+1.)))*B_unpms[k,m,n-2]
+ for n in range(0,n_max+1):
+ m=0
+ B_snpms[k,m,n] = 1./float(m+1)*np.sqrt((float(n+m)+1.)*(float(n-m))) * np.sqrt(1.-u**2) * B_unpms[k,m+1,n] + u*B_unpms[k,m,n]
+ for n in range(1,n_max+1):
+ for m in range(1,n+1):
+ B_snpms[k,m,n] = float(n)/float(m) * u * B_unpms[k,m,n] - float(m+n)/float(m) * \
+ np.sqrt( ( (2.*float(n)+1.)*(float(n)-float(m))*(float(n)-1.) ) / ( (2.*float(n)-1.)*(float(n)+float(m))*(float(n)+1.) ) )*B_unpms[k,m,n-1]
+ B_Pnmms = np.zeros_like(B_Pnpms)
+ B_unmms = np.zeros_like(B_unpms)
+ B_snmms = np.zeros_like(B_snpms)
+ for m in range(m_max+1):
+ B_Pnmms[:,m,:] = (-1.0)**m * B_Pnpms[:,m,:]
+ B_unmms[:,m,:] = (-1.0)**(m+1) * B_unpms[:,m,:]
+ B_snmms[:,m,:] = (-1.0)**(m) * B_snpms[:,m,:]
+ B_Pnmms=B_Pnmms[:,::-1,:]
+ B_unmms=B_unmms[:,::-1,:]
+ B_snmms=B_snmms[:,::-1,:]
+ B_Pnms = np.concatenate((B_Pnmms,B_Pnpms[:,1::,:]),axis=1)
+ B_unms = np.concatenate((B_unmms,B_unpms[:,1::,:]),axis=1)
+ B_snms = np.concatenate((B_snmms,B_snpms[:,1::,:]),axis=1)
+ E_scat_onsphere_sph = np.array( [np.loadtxt(post_filename,usecols=[8])
+ + 1j*np.loadtxt(post_filename,usecols=[11]),
+ np.loadtxt(post_filename,usecols=[9])
+ + 1j*np.loadtxt(post_filename,usecols=[12]),
+ np.loadtxt(post_filename,usecols=[10])
+ + 1j*np.loadtxt(post_filename,usecols=[13])])
+ E_scat_onsphere_sph_r = E_scat_onsphere_sph[0,:].reshape(npts_phi,npts_theta,order='F').transpose()
+ E_scat_onsphere_sph_t = E_scat_onsphere_sph[1,:].reshape(npts_phi,npts_theta,order='F').transpose()
+ E_scat_onsphere_sph_p = E_scat_onsphere_sph[2,:].reshape(npts_phi,npts_theta,order='F').transpose()
+ for ko in range(p_max):
+ po = ps[ko]
+ n = int(np.sqrt(po))
+ m = n*(n+1) - int(po)
+ # print('=========>> po',po,'n',n,'m',m)
+ B_PnmN_costheta = np.tile( B_Pnms[:,m+m_max,n],(npts_phi,1)).transpose()
+ B_UnmN_costheta = np.tile( B_unms[:,m+m_max,n],(npts_phi,1)).transpose()
+ B_SnmN_costheta = np.tile( B_snms[:,m+m_max,n],(npts_phi,1)).transpose()
+ B_Ynm_r = B_PnmN_costheta * np.exp(1j*float(m)*phi_sph)
+ B_Ynm_t = np.zeros_like(phi_sph)
+ B_Ynm_p = np.zeros_like(phi_sph)
+ B_Xnm_r = np.zeros_like(phi_sph)
+ B_Xnm_t = 1j * B_UnmN_costheta * np.exp(1j*float(m)*phi_sph)
+ B_Xnm_p = -1. * B_SnmN_costheta * np.exp(1j*float(m)*phi_sph)
+ B_Znm_r = np.zeros_like(phi_sph)
+ B_Znm_t = 1. * B_SnmN_costheta * np.exp(1j*float(m)*phi_sph)
+ B_Znm_p = 1j * B_UnmN_costheta * np.exp(1j*float(m)*phi_sph)
+ B_erCrossXnm_r = np.zeros_like(phi_sph,dtype=complex)
+ B_erCrossXnm_t = -B_Xnm_p
+ B_erCrossXnm_p = B_Xnm_t
+ FF_Xnm_t[:,:,ko] = B_Xnm_t
+ FF_Xnm_p[:,:,ko] = B_Xnm_p
+ FF_erCrossXnm_t[:,:,ko] = B_erCrossXnm_t
+ FF_erCrossXnm_p[:,:,ko] = B_erCrossXnm_p
+
+
+ sph_bessel_n_ofkr = np.sqrt(pi/(2.*k_Out*r_sph))*outgoing_sph_hankel(float(n )+0.5,k_Out*r_sph)
+ sph_bessel_nminus1_ofkr = np.sqrt(pi/(2.*k_Out*r_sph))*outgoing_sph_hankel(float(n-1)+0.5,k_Out*r_sph)
+ dRicatti_dx_ofkr = (k_Out * r_sph * (sph_bessel_nminus1_ofkr-(float(n+1)/((k_Out*r_sph))) * sph_bessel_n_ofkr) + sph_bessel_n_ofkr)
+ B_Mnm_r = 0.
+ B_Mnm_t = sph_bessel_n_ofkr * B_Xnm_t
+ B_Mnm_p = sph_bessel_n_ofkr * B_Xnm_p
+ B_Nnm_r = 1./(k_Out*r_sph) * np.sqrt(float(n*(n+1))) * sph_bessel_n_ofkr * B_Ynm_r
+ B_Nnm_t = 1./(k_Out*r_sph) * dRicatti_dx_ofkr * B_Znm_t
+ B_Nnm_p = 1./(k_Out*r_sph) * dRicatti_dx_ofkr * B_Znm_p
+
+ B_EdotconjYnm = E_scat_onsphere_sph_r*B_Ynm_r.conjugate()
+ B_EdotconjZnm = E_scat_onsphere_sph_t*B_Znm_t.conjugate() + E_scat_onsphere_sph_p*B_Znm_p.conjugate()
+ B_EdotconjXnm = E_scat_onsphere_sph_t*B_Xnm_t.conjugate() + E_scat_onsphere_sph_p*B_Xnm_p.conjugate()
+ normalize_fhnm_X = 1./sph_bessel_n_ofkr
+ normalize_fenm_Y = k_Out*r_sph/(sph_bessel_n_ofkr*np.sqrt(float(n)*(float(n)+1.)) )
+ normalize_fenm_Z = k_Out*r_sph/dRicatti_dx_ofkr
+ fenm_Y[int(po)-1] = np.trapz(np.trapz((np.sin(theta_sph)*B_EdotconjYnm).transpose(),theta_sph[:,0]),phi_sph[0,:])*normalize_fenm_Y
+ fenm_Z[int(po)-1] = np.trapz(np.trapz((np.sin(theta_sph)*B_EdotconjZnm).transpose(),theta_sph[:,0]),phi_sph[0,:])*normalize_fenm_Z
+ fhnm_X[int(po)-1] = np.trapz(np.trapz((np.sin(theta_sph)*B_EdotconjXnm).transpose(),theta_sph[:,0]),phi_sph[0,:])*normalize_fhnm_X
+
+ EdotconjMnm = E_scat_onsphere_sph_r*B_Mnm_r.conjugate() + E_scat_onsphere_sph_t*B_Mnm_t.conjugate() + E_scat_onsphere_sph_p*B_Mnm_p.conjugate()
+ EdotconjNnm = E_scat_onsphere_sph_r*B_Nnm_r.conjugate() + E_scat_onsphere_sph_t*B_Nnm_t.conjugate() + E_scat_onsphere_sph_p*B_Nnm_p.conjugate()
+ MnmconjMnm = B_Mnm_r*B_Mnm_r.conjugate() + B_Mnm_t*B_Mnm_t.conjugate() + B_Mnm_p*B_Mnm_p.conjugate()
+ NnmconjNnm = B_Nnm_r*B_Nnm_r.conjugate() + B_Nnm_t*B_Nnm_t.conjugate() + B_Nnm_p*B_Nnm_p.conjugate()
+ MnmconjNnm = B_Mnm_r*B_Nnm_r.conjugate() + B_Mnm_t*B_Nnm_t.conjugate() + B_Mnm_p*B_Nnm_p.conjugate()
+ XnmconjXnm = B_Xnm_r*B_Xnm_r.conjugate() + B_Xnm_t*B_Xnm_t.conjugate() + B_Xnm_p*B_Xnm_p.conjugate()
+ YnmconjYnm = B_Ynm_r*B_Ynm_r.conjugate() + B_Ynm_t*B_Ynm_t.conjugate() + B_Ynm_p*B_Ynm_p.conjugate()
+ ZnmconjZnm = B_Znm_r*B_Znm_r.conjugate() + B_Znm_t*B_Znm_t.conjugate() + B_Znm_p*B_Znm_p.conjugate()
+ XnmconjYnm = B_Xnm_r*B_Ynm_r.conjugate() + B_Xnm_t*B_Ynm_t.conjugate() + B_Xnm_p*B_Ynm_p.conjugate()
+ YnmconjZnm = B_Ynm_r*B_Znm_r.conjugate() + B_Ynm_t*B_Znm_t.conjugate() + B_Ynm_p*B_Znm_p.conjugate()
+ ZnmconjXnm = B_Znm_r*B_Xnm_r.conjugate() + B_Znm_t*B_Xnm_t.conjugate() + B_Znm_p*B_Xnm_p.conjugate()
+ normalize_aM_nm2 = np.trapz(np.trapz((np.sin(theta_sph)*MnmconjMnm).transpose(),theta_sph[:,0]),phi_sph[0,:])
+ normalize_bN_nm2 = np.trapz(np.trapz((np.sin(theta_sph)*NnmconjNnm).transpose(),theta_sph[:,0]),phi_sph[0,:])
+ orth = np.trapz(np.trapz((np.sin(theta_sph)*MnmconjNnm).transpose(),theta_sph[:,0]),phi_sph[0,:])
+ orthX = np.trapz(np.trapz((np.sin(theta_sph)*XnmconjXnm).transpose(),theta_sph[:,0]),phi_sph[0,:])
+ orthY = np.trapz(np.trapz((np.sin(theta_sph)*YnmconjYnm).transpose(),theta_sph[:,0]),phi_sph[0,:])
+ orthZ = np.trapz(np.trapz((np.sin(theta_sph)*ZnmconjZnm).transpose(),theta_sph[:,0]),phi_sph[0,:])
+ orth1 = np.trapz(np.trapz((np.sin(theta_sph)*XnmconjYnm).transpose(),theta_sph[:,0]),phi_sph[0,:])
+ orth2 = np.trapz(np.trapz((np.sin(theta_sph)*YnmconjZnm).transpose(),theta_sph[:,0]),phi_sph[0,:])
+ orth3 = np.trapz(np.trapz((np.sin(theta_sph)*ZnmconjXnm).transpose(),theta_sph[:,0]),phi_sph[0,:])
+ aM_nm[int(po)-1] = np.trapz(np.trapz((np.sin(theta_sph)*EdotconjMnm).transpose(),theta_sph[:,0]),phi_sph[0,:]) / normalize_aM_nm2
+ bN_nm[int(po)-1] = np.trapz(np.trapz((np.sin(theta_sph)*EdotconjNnm).transpose(),theta_sph[:,0]),phi_sph[0,:]) / normalize_bN_nm2
+ return [fhnm_X,fenm_Y,fenm_Z,aM_nm,bN_nm,FF_Xnm_t,FF_Xnm_p,FF_erCrossXnm_t,FF_erCrossXnm_p]
+
+###########################
+def plot_farfield(far_field_sph,filename):
+ vmax_sph = np.max(far_field_sph)
+ vmin_sph = np.min(far_field_sph)
+ x_sph = far_field_sph/vmax_sph*np.sin(theta_sph) * np.cos(phi_sph)
+ y_sph = far_field_sph/vmax_sph*np.sin(theta_sph) * np.sin(phi_sph)
+ z_sph = far_field_sph/vmax_sph*np.cos(theta_sph)
+ fig = pl.figure()
+ ax = fig.add_subplot(111, projection='3d')
+ ax.set_aspect('equal')
+ surf=ax.plot_surface(x_sph,y_sph,z_sph,\
+ facecolors=cm.viridis(far_field_sph/vmax_sph),\
+ rstride=2,\
+ cstride=2,\
+ linewidth=0,\
+ vmin = vmin_sph,\
+ vmax = vmax_sph,\
+ shade=True,\
+ alpha=0.5,\
+ antialiased=False)
+ cset = ax.contourf(x_sph, y_sph, z_sph,1,fc='k', zdir='z', offset=-1)
+ cset = ax.contourf(x_sph, y_sph, z_sph,1,fc='k', zdir='x', offset=-1)
+ cset = ax.contourf(x_sph, y_sph, z_sph,1,fc='k', zdir='y', offset=1 )
+ surf.set_edgecolor('k')
+ max_range = 0.5*np.max(np.array([x_sph.max()-x_sph.min(), y_sph.max()-y_sph.min(), z_sph.max()-z_sph.min()]))
+ mid_x = (x_sph.max()+x_sph.min())*0.5
+ mid_y = (y_sph.max()+y_sph.min())*0.5
+ mid_z = (z_sph.max()+z_sph.min())*0.5
+ ax.set_xlim(mid_x-max_range, mid_x+max_range)
+ ax.set_ylim(mid_y-max_range, mid_y+max_range)
+ ax.set_zlim(mid_z-max_range, mid_z+max_range)
+ ax.xaxis.set_ticklabels([]);ax.xaxis.set_label('x')
+ ax.yaxis.set_ticklabels([]);ax.yaxis.set_label('y')
+ ax.zaxis.set_ticklabels([]);ax.zaxis.set_label('z')
+ pl.savefig(filename,bbox_inches='tight')
+ pl.close('all')
diff --git a/ElectromagneticScattering/scattering_post.py b/ElectromagneticScattering/scattering_post.py
new file mode 100644
index 0000000000000000000000000000000000000000..7e52c8538578872efa6ff9abdab109813e6ef708
--- /dev/null
+++ b/ElectromagneticScattering/scattering_post.py
@@ -0,0 +1,166 @@
+"""
+///////////////////////////////
+// Author : Guillaume Demesy //
+// scattering_post.py //
+///////////////////////////////
+"""
+
+Tmatrix_rfull_ab = np.zeros((2*p_max,2*p_max,nb_cuts),dtype=complex)
+Tmatrix_rfull_fy = np.zeros((2*p_max,2*p_max,nb_cuts),dtype=complex)
+Tmatrix_rfull_fz = np.zeros((2*p_max,2*p_max,nb_cuts),dtype=complex)
+tab_E_NTF_t_G = np.zeros((npts_theta,npts_phi,3),dtype=complex)
+tab_E_NTF_p_G = np.zeros((npts_theta,npts_phi,3),dtype=complex)
+tab_E_NTF_t_PW = np.zeros((npts_theta,npts_phi),dtype=complex)
+tab_E_NTF_p_PW = np.zeros((npts_theta,npts_phi),dtype=complex)
+
+my_dir='./run_results/'
+if flag_study==0:
+ nbNM = 1
+ p_max_in=1
+elif flag_study==1:
+ nbNM = 2
+ p_max_in=p_max
+elif flag_study==2:
+ nbNM = 1
+ p_max_in=p_max
+ nb_cuts = 3
+elif flag_study==3:
+ p_max_in=0
+
+if siwt==1:
+ outgoing_sph_hankel = hankel2
+else:
+ outgoing_sph_hankel = hankel1
+
+print('Postprocessing...')
+for ke in range(p_max_in):
+ for isN in range(nbNM):
+ pe = ps[ke]
+ ne =int(np.sqrt(pe))
+ me = ne*(ne+1) - int(pe)
+ aM_nm = np.zeros((nb_cuts,p_max),dtype=complex)
+ bN_nm = np.zeros((nb_cuts,p_max),dtype=complex)
+ fenm_Y = np.zeros((nb_cuts,p_max),dtype=complex)
+ fenm_Z = np.zeros((nb_cuts,p_max),dtype=complex)
+ fhnm_X = np.zeros((nb_cuts,p_max),dtype=complex)
+ for nr in range(nb_cuts):
+ r_sph = r_sphs[nr]
+ # print('======>> postprocessing r_sph',r_sph)
+ par_gmsh_getdp=open('parameters_gmsh_getdp.dat', 'a')
+ par_gmsh_getdp.write('r_sph = %3.15e;\n' %(r_sph) )
+ par_gmsh_getdp.close()
+ if flag_study==1:
+ if isN==0:post_filename = my_dir+'E_scat_onsphere_sph_M%g_r%g.dat'%((pe,nr))
+ else: post_filename = my_dir+'E_scat_onsphere_sph_N%g_r%g.dat'%((pe,nr))
+ elif flag_study==0:
+ post_filename = my_dir+'E_scat_onsphere_sph_PW_r%g.dat'%(nr)
+ elif flag_study==2:
+ post_filename = my_dir+'E_tot_onsphere_sph_G%g.dat'%(nr)
+ res=field_VSH_expansion(post_filename)
+ fhnm_X[nr,:] = res[0]
+ fenm_Y[nr,:] = res[1]
+ fenm_Z[nr,:] = res[2]
+ aM_nm[nr,:] = res[3]
+ bN_nm[nr,:] = res[4]
+ FF_Xnm_t = res[5]
+ FF_Xnm_p = res[6]
+ FF_erCrossXnm_t = res[7]
+ FF_erCrossXnm_p = res[8]
+ if flag_study==1:
+ if isN==0:
+ Tmatrix_rfull_ab[ke,0:p_max,:] = aM_nm.transpose()
+ Tmatrix_rfull_ab[ke,p_max:2*p_max,:] = bN_nm.transpose()
+ Tmatrix_rfull_fy[ke,0:p_max,:] = fhnm_X.transpose()
+ Tmatrix_rfull_fy[ke,p_max:2*p_max,:] = fenm_Y.transpose()
+ Tmatrix_rfull_fz[ke,0:p_max,:] = fhnm_X.transpose()
+ Tmatrix_rfull_fz[ke,p_max:2*p_max,:] = fenm_Z.transpose()
+ if isN==1:
+ Tmatrix_rfull_ab[p_max+ke,0:p_max,:] = aM_nm.transpose()
+ Tmatrix_rfull_ab[p_max+ke,p_max:2*p_max,:] = bN_nm.transpose()
+ Tmatrix_rfull_fy[p_max+ke,0:p_max,:] = fhnm_X.transpose()
+ Tmatrix_rfull_fy[p_max+ke,p_max:2*p_max,:] = fenm_Y.transpose()
+ Tmatrix_rfull_fz[p_max+ke,0:p_max,:] = fhnm_X.transpose()
+ Tmatrix_rfull_fz[p_max+ke,p_max:2*p_max,:] = fenm_Z.transpose()
+ Tmatrix_ab = np.mean(Tmatrix_rfull_ab,axis=2)
+ Tmatrix_fy = np.mean(Tmatrix_rfull_fy,axis=2)
+ Tmatrix_fz = np.mean(Tmatrix_rfull_fz,axis=2)
+ std_Tmatrix_ab = np.std(Tmatrix_rfull_ab,axis=2)
+ std_Tmatrix_fy = np.std(Tmatrix_rfull_fy,axis=2)
+ std_Tmatrix_fz = np.std(Tmatrix_rfull_fz,axis=2)
+ np.savez('T_matrix_ellips_sphPML_epsre_%.2lf_eps_im_%.2lf.npz'%(epsr_In.real,epsr_In.imag), \
+ Tmatrix_rfull_ab = Tmatrix_rfull_ab, \
+ Tmatrix_rfull_fy = Tmatrix_rfull_fy, \
+ Tmatrix_rfull_fz = Tmatrix_rfull_fz,\
+ r_sphs = r_sphs,\
+ n_max = n_max,\
+ ps = ps,\
+ p_max = p_max,\
+ lambda0 = lambda0, \
+ Tmatrix_ab = Tmatrix_ab, \
+ Tmatrix_fy = Tmatrix_fy, \
+ Tmatrix_fz = Tmatrix_fz, \
+ std_Tmatrix_ab = std_Tmatrix_ab,\
+ std_Tmatrix_fy = std_Tmatrix_fy,\
+ std_Tmatrix_fz = std_Tmatrix_fz)
+ if flag_study==0:
+ mean_aM_nm = np.mean(aM_nm ,axis=0)
+ mean_bN_nm = np.mean(bN_nm ,axis=0)
+ mean_fenm_Y = np.mean(fenm_Y,axis=0)
+ mean_fenm_Z = np.mean(fenm_Z,axis=0)
+ mean_fhnm_X = np.mean(fhnm_X,axis=0)
+ std_aM_nm = np.std( aM_nm ,axis=0)
+ std_bN_nm = np.std( bN_nm ,axis=0)
+ std_fenm_Y = np.std( fenm_Y,axis=0)
+ std_fenm_Z = np.std( fenm_Z,axis=0)
+ std_fhnm_X = np.std( fhnm_X,axis=0)
+ for ko in range(p_max):
+ tab_E_NTF_t_PW += (1j)**ko * mean_aM_nm[ko]*1j*FF_Xnm_t[:,:,ko] + (1j)**ko * mean_bN_nm[ko]*FF_erCrossXnm_t[:,:,ko]
+ tab_E_NTF_p_PW += (1j)**ko * mean_aM_nm[ko]*1j*FF_Xnm_p[:,:,ko] + (1j)**ko * mean_bN_nm[ko]*FF_erCrossXnm_p[:,:,ko]
+ I_NTF_PW = (np.abs(tab_E_NTF_t_PW))**2 + (np.abs(tab_E_NTF_p_PW))**2
+ plot_farfield(I_NTF_PW[:,:],my_dir+'farfield_PW.pdf')
+ if flag_study==2:
+ tens_Green_scat=np.zeros((3,3))
+ for nr in range(nb_cuts):
+ tens_Green_scat[:,nr] = np.loadtxt(my_dir+'E_scat_im_onpt_G%g.dat'%(nr))[8:11]
+ for ko in range(p_max):
+ tab_E_NTF_t_G[:,:,nr] += (1j)**ko * aM_nm[nr,ko]*1j*FF_Xnm_t[:,:,ko] + (1j)**ko * bN_nm[nr,ko]*FF_erCrossXnm_t[:,:,ko]
+ tab_E_NTF_p_G[:,:,nr] += (1j)**ko * aM_nm[nr,ko]*1j*FF_Xnm_p[:,:,ko] + (1j)**ko * bN_nm[nr,ko]*FF_erCrossXnm_p[:,:,ko]
+ I_NTF_G = (np.abs(tab_E_NTF_t_G))**2 + (np.abs(tab_E_NTF_p_G))**2
+ tens_Green_inc = -siwt*np.eye(3)*k_Out/(6.*pi)
+ tens_Green_tot = tens_Green_scat+tens_Green_inc
+ for k in range(3):plot_farfield(I_NTF_G[:,:,k],my_dir+'farfield_green%g.pdf'%(k))
+
+if flag_study==1:
+ print('Tmatrix_ab:')
+ print(Tmatrix_ab)
+ print('std_Tmatrix_ab:')
+ print(std_Tmatrix_ab)
+ print('Tmatrix_fy:')
+ print(Tmatrix_fy)
+ print('np.abs(2*Tmatrix_fy+1):')
+ print(np.abs(2*Tmatrix_fy+1))
+ print('Tmatrix_fz:')
+ print(Tmatrix_fz)
+ print('np.abs(2*Tmatrix_fz+1):')
+ print(np.abs(2*Tmatrix_fz+1))
+ print('np.abs(2*Tmatrix_ab+1):')
+ print(np.abs(2*Tmatrix_ab+1))
+elif flag_study==0:
+ print('anm:')
+ print(mean_aM_nm)
+ print('bnm')
+ print(mean_bN_nm)
+elif flag_study==2:
+ print('tens_Green_tot')
+ print(tens_Green_tot)
+ print('tens_Green_tot/G0')
+ print(tens_Green_tot/tens_Green_inc[0,0])
+
+elif flag_study==3:
+ eigs = np.loadtxt(my_dir+'EigenValuesReal.txt',usecols=[5]) + 1j*np.loadtxt(my_dir+'EigenValuesImag.txt',usecols=[5])
+ eigs *= (nm/1e-9)**2
+ pl.savez(my_dir+'eigs.npz',eigs=eigs)
+ pl.figure(figsize=(20,15))
+ pl.plot( eigs.real , eigs.imag,'+',mfc='none')
+ pl.grid()
+ pl.savefig(my_dir+'eigenvalues.pdf')
\ No newline at end of file
diff --git a/ElectromagneticScattering/screenshot1_512.png b/ElectromagneticScattering/screenshot1_512.png
new file mode 100644
index 0000000000000000000000000000000000000000..410e93a1a5a34c81866c3761356dd145698aeafe
Binary files /dev/null and b/ElectromagneticScattering/screenshot1_512.png differ
diff --git a/NonLinearEVP/NonLinearEVP.geo b/NonLinearEVP/NonLinearEVP.geo
new file mode 100644
index 0000000000000000000000000000000000000000..86769a45f1f6d8c124a88f549a9b3aada2c4eae5
--- /dev/null
+++ b/NonLinearEVP/NonLinearEVP.geo
@@ -0,0 +1,359 @@
+///////////////////////////////////
+//// Author : Guillaume Demesy ////
+//// .geo ////
+///////////////////////////////////
+
+Include "NonLinearEVP_data.geo";
+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;
+
+
+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) = {-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
+EndIf
+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};
+ 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};
+ Point(16) = { d_sq/2.-epsc , d_sq/2.+epsc, 0. , lc_sqa};
+ Point(17) = { d_sq/2. , d_sq/2.+epsc, 0. , lc_sqa};
+ Point(18) = { d_sq/2.+epsc , d_sq/2.+epsc, 0. , lc_sqa};
+ Point(19) = { d_sq/2.+epsc , d_sq/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(22) = { d_sq/2.+epsc ,-d_sq/2. , 0. , lc_sqa};
+ Point(23) = { d_sq/2.+epsc ,-d_sq/2.-epsc, 0. , lc_sqa};
+ Point(24) = { d_sq/2. ,-d_sq/2.-epsc, 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. ,-d_sq/2.-epsc, 0. , lc_sqa};
+ Point(28) = {-d_sq/2.-epsc ,-d_sq/2.-epsc, 0. , lc_sqa};
+ Point(29) = {-d_sq/2.-epsc ,-d_sq/2. , 0. , lc_sqa};
+ 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};
+ 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};
+ Line(36) = {45, 45};
+ Line(37) = {39, 39};
+ Line(38) = {45, 39};
+ Line(39) = {39, 7};
+ Line(40) = {7, 48};
+ Line(41) = {48, 47};
+ Line(42) = {47, 6};
+ Line(43) = {6, 38};
+ Line(44) = {38, 37};
+ Line(45) = {37, 5};
+ Line(46) = {5, 46};
+ Line(47) = {46, 49};
+ Line(48) = {49, 8};
+ Line(49) = {31, 49};
+ Line(50) = {49, 36};
+ Line(51) = {36, 45};
+ Line(52) = {45, 15};
+ Line(53) = {32, 8};
+ Line(54) = {8, 14};
+ Line(55) = {13, 8};
+ Line(56) = {8, 36};
+ Line(57) = {31, 8};
+ Line(58) = {8, 15};
+ Line(59) = {16, 39};
+ Line(60) = {39, 35};
+ Line(61) = {35, 48};
+ Line(62) = {48, 20};
+ Line(63) = {19, 7};
+ Line(64) = {7, 17};
+ Line(65) = {16, 7};
+ Line(66) = {7, 20};
+ Line(67) = {35, 7};
+ Line(68) = {7, 18};
+ Line(69) = {34, 47};
+ Line(70) = {47, 21};
+ Line(71) = {34, 38};
+ Line(72) = {38, 25};
+ Line(73) = {25, 6};
+ Line(74) = {6, 24};
+ Line(75) = {6, 22};
+ Line(76) = {6, 21};
+ Line(77) = {23, 6};
+ Line(78) = {6, 34};
+ Line(79) = {30, 46};
+ Line(80) = {46, 33};
+ Line(81) = {33, 37};
+ Line(82) = {37, 26};
+ Line(83) = {26, 5};
+ Line(84) = {5, 27};
+ Line(85) = {5, 28};
+ Line(86) = {5, 29};
+ Line(87) = {5, 30};
+ Line(88) = {5, 33};
+ Line(89) = {33, 36};
+ Line(90) = {36, 35};
+ Line(91) = {35, 34};
+ Line(92) = {34, 33};
+ Line(93) = {13, 14};
+ Line(94) = {14, 15};
+ Line(95) = {15, 16};
+ Line(96) = {16, 17};
+ Line(97) = {17, 18};
+ Line(98) = {18, 19};
+ Line(99) = {19, 20};
+ Line(100) = {20, 21};
+ Line(101) = {21, 22};
+ Line(102) = {22, 23};
+ Line(103) = {23, 24};
+ Line(104) = {24, 25};
+ Line(105) = {25, 26};
+ Line(106) = {26, 27};
+ Line(107) = {27, 28};
+ Line(108) = {28, 29};
+ Line(109) = {29, 30};
+ Line(110) = {30, 31};
+ Line(111) = {31, 32};
+ Line(112) = {32, 13};
+ Line(113) = {1, 55};
+ Line(114) = {55, 54};
+ Line(115) = {54, 4};
+ Line(116) = {2, 58};
+ Line(117) = {58, 56};
+ Line(118) = {56, 3};
+
+ Line Loop(1) = {53, -55, -112};
+ Plane Surface(1) = {1};
+ Line Loop(2) = {55, 54, -93};
+ Plane Surface(2) = {2};
+ Line Loop(3) = {54, 94, -58};
+ Plane Surface(3) = {-3};
+ Line Loop(4) = {35, 52, -58};
+ Plane Surface(4) = {4};
+ Line Loop(5) = {111, 53, -57};
+ Plane Surface(5) = {-5};
+ Line Loop(6) = {49, 48, -57};
+ Plane Surface(6) = {6};
+ Line Loop(7) = {48, 56, -50};
+ Plane Surface(7) = {-7};
+ Line Loop(8) = {56, 51, -35};
+ Plane Surface(8) = {8};
+ Line Loop(9) = {51, 38, 60, -90};
+ Plane Surface(9) = {-9};
+ Line Loop(10) = {38, -59, -95, -52};
+ Plane Surface(10) = {10};
+ Line Loop(11) = {59, 39, -65};
+ Plane Surface(11) = {11};
+ Line Loop(12) = {65, 64, -96};
+ Plane Surface(12) = {12};
+ Line Loop(13) = {64, 97, -68};
+ Plane Surface(13) = {-13};
+ Line Loop(14) = {68, 98, 63};
+ Plane Surface(14) = {-14};
+ Line Loop(15) = {63, 66, -99};
+ Plane Surface(15) = {15};
+ Line Loop(16) = {66, -62, -40};
+ Plane Surface(16) = {-16};
+ Line Loop(17) = {40, -61, 67};
+ Plane Surface(17) = {-17};
+ Line Loop(18) = {67, -39, 60};
+ Plane Surface(18) = {18};
+ Line Loop(19) = {91, 69, -41, -61};
+ Plane Surface(19) = {19};
+ Line Loop(20) = {70, -100, -62, 41};
+ Plane Surface(20) = {20};
+ Line Loop(21) = {42, 76, -70};
+ Plane Surface(21) = {21};
+ Line Loop(22) = {76, 101, -75};
+ Plane Surface(22) = {-22};
+ Line Loop(23) = {42, 78, 69};
+ Plane Surface(23) = {-23};
+ Line Loop(24) = {71, -43, 78};
+ Plane Surface(24) = {24};
+ Line Loop(25) = {72, 73, 43};
+ Plane Surface(25) = {25};
+ Line Loop(26) = {73, 74, 104};
+ Plane Surface(26) = {-26};
+ Line Loop(27) = {74, -103, 77};
+ Plane Surface(27) = {27};
+ Line Loop(28) = {77, 75, 102};
+ Plane Surface(28) = {-28};
+ Line Loop(29) = {92, 81, -44, -71};
+ Plane Surface(29) = {29};
+ Line Loop(30) = {82, -105, -72, 44};
+ Plane Surface(30) = {30};
+ Line Loop(31) = {83, -45, 82};
+ Plane Surface(31) = {-31};
+ Line Loop(32) = {83, 84, -106};
+ Plane Surface(32) = {32};
+ Line Loop(33) = {84, 107, -85};
+ Plane Surface(33) = {-33};
+ Line Loop(34) = {85, 108, -86};
+ Plane Surface(34) = {-34};
+ Line Loop(35) = {86, 109, -87};
+ Plane Surface(35) = {-35};
+ Line Loop(36) = {87, 79, -46};
+ Plane Surface(36) = {-36};
+ Line Loop(37) = {46, 80, -88};
+ Plane Surface(37) = {-37};
+ Line Loop(38) = {88, 81, 45};
+ Plane Surface(38) = {-38};
+ 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};
+ Line Loop(46) = {106, 107, 108, 109, 110, 111, 112, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105};
+ 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 { expression-list } | "*" < = { expression-list } > < Left | Right | Alternate | AlternateRight | AlternateLeft > ;
+ Transfinite Surface { 9 } Left;
+ Transfinite Surface { 10 } Right;
+
+ Transfinite Surface { 19 } Left;
+ Transfinite Surface { 20 } Left;
+
+ Transfinite Surface { 29 } Left;
+ Transfinite Surface { 30 } Left;
+
+ Transfinite Surface { 39 } Left;
+ Transfinite Surface { 40 } Right;
+
+
+ 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
+
diff --git a/NonLinearEVP/NonLinearEVP.pro b/NonLinearEVP/NonLinearEVP.pro
new file mode 100644
index 0000000000000000000000000000000000000000..fc8ac7652a477abf1909b0c706b7f07108c20205
--- /dev/null
+++ b/NonLinearEVP/NonLinearEVP.pro
@@ -0,0 +1,539 @@
+///////////////////////////////////
+//// Author : Guillaume Demesy ////
+//// .pro ////
+///////////////////////////////////
+
+Include "NonLinearEVP_data.geo";
+
+Group {
+ // Domains
+ Om_1 = Region[100] ; //in (dispersive)
+ Om_2inter = Region[101]; // out (non dispersive)
+ pmltop = Region[102]; // pml top
+ pmlbot = Region[103]; // pml top
+ Om = Region[{Om_1,Om_2inter,pmltop,pmlbot}];
+ Om_nopml = Region[{Om_1,Om_2inter}];
+ Om_pml = Region[{pmltop,pmlbot}];
+ Om_2 = Region[{Om_pml,Om_2inter}]; // non dispersive
+
+ // Boundaries
+ SurfBlochLeft = Region[1001];
+ SurfBlochRight = Region[1002];
+ Surfpml = Region[{1003,1004}];
+ Bound = Region[{1005}];
+ PrintPoint = Region[10000];
+}
+
+Function{
+ I[] = Complex[0,1];
+
+ siwt = -1.;
+ a_pml = 1.;
+ b_pml = -siwt;
+
+ sx[Om_pml] = 1.;
+ sy[Om_pml] = Complex[a_pml,b_pml];
+ sx[Om_nopml] = 1.;
+ sy[Om_nopml] = 1.;
+
+ sz = 1;
+
+ PML_Tensor[] = TensorDiag[sz*sy[]/sx[],sx[]*sz/sy[],sx[]*sy[]/sz];
+ mur[] = TensorDiag[1,1,1]*PML_Tensor[];
+ epsr_nod[] = TensorDiag[1,1,1]*PML_Tensor[];
+
+ term_o3_3cc_1[] = Complex[-eps_oo_1,0];
+ term_o3_2cc_1[] = Complex[0,-eps_oo_1*gam_1];
+ term_o3_1cc_1[] = Complex[om_d_1^2,0];
+ term_o3_1rr_1[] = cel^2*Complex[1,0];
+ term_o3_0rr_1[] = cel^2*Complex[0,gam_1];
+ term_o3_3cc_2[] = -1*epsr_nod[];
+ term_o3_2cc_2[] = -1*epsr_nod[]*Complex[0,gam_2];
+ term_o3_1cc_2[] = Complex[om_d_2^2,0];
+ term_o3_1rr_2[] = cel^2*Complex[1,0];
+ term_o3_0rr_2[] = cel^2*Complex[0,gam_2];
+ term_o3_3cc[Om_1] = term_o3_3cc_1[];
+ term_o3_2cc[Om_1] = term_o3_2cc_1[];
+ term_o3_1cc[Om_1] = term_o3_1cc_1[];
+ term_o3_1rr[Om_1] = term_o3_1rr_1[];
+ term_o3_0rr[Om_1] = term_o3_0rr_1[];
+ term_o3_3cc[Om_2] = term_o3_3cc_2[];
+ term_o3_2cc[Om_2] = term_o3_2cc_2[];
+ term_o3_1cc[Om_2] = term_o3_1cc_2[];
+ term_o3_1rr[Om_2] = term_o3_1rr_2[];
+ term_o3_0rr[Om_2] = term_o3_0rr_2[];
+
+
+ dephx[] = Complex[ Cos[kx*a_lat] , Sin[kx*a_lat]];
+ // dephy[] = Complex[ Cos[ky*a_lat] , Sin[ky*a_lat]];
+
+ // bidon
+ EigFilter[] = (Norm[$EigenvalueReal] > 1e-150);
+
+ //Auxialiary Field resolution coefficients
+ eps_inf[] = TensorDiag[1.,1.,1.]*eps_oo_1;
+ Om_D[] = TensorDiag[1,1,1]*Sqrt[om_d_1^2/( 2. )];
+ Gamma_d[] = TensorDiag[1,1,1]*gam_1;
+}
+
+Constraint {
+ {Name BlochX;
+ Case {
+ { Region SurfBlochRight; Type LinkCplx ; RegionRef SurfBlochLeft;
+ Coefficient dephx[]; Function Vector[$X-a_lat,$Y,$Z] ;
+ }
+ }
+ }
+
+ {Name Dirichlet; Type Assign;
+ Case {
+ { Region Surfpml ; Value 0.; }
+ }
+ }
+
+ // div condition
+ { Name Constraint_cInt; Type Assign;
+ Case {
+ }
+ }
+ { Name Constraint_Int; Type Assign;
+ Case {
+ { Region Bound; Value 0.; }
+ }
+ }
+}
+
+
+
+Jacobian {
+ {Name JVol ;
+ Case {
+ { Region All ; Jacobian Vol ; }
+ }
+ }
+ { Name JSur ;
+ Case{
+ { Region All ; Jacobian Sur ; }
+ }
+ }
+ { Name JLin ;
+ Case{
+ { Region All ; Jacobian Lin ; }
+ }
+ }
+}
+
+Integration {
+ { Name Int_1 ;
+ Case {
+ { Type Gauss ;
+ Case {
+ { GeoElement Point ; NumberOfPoints 1 ; }
+ { GeoElement Line ; NumberOfPoints 4 ; }
+ { GeoElement Triangle ; NumberOfPoints 6 ; }
+ }
+ }
+ }
+ }
+}
+
+FunctionSpace {
+ { 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]; }
+ 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; }
+ EndIf
+ }
+ }
+
+ { Name Hcurl; Type Form1;
+ 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)
+ { 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]; }
+ // BF_Edge_3F_a, BF_Edge_3F_b, BF_Edge_4E
+ EndIf
+ }
+ Constraint {
+ { NameOfCoef un; EntityType EdgesOf ; NameOfConstraint BlochX; }
+ { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint BlochX; }
+ { NameOfCoef un; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+ { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+ If(flag_o2==1)
+ { NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint BlochX; }
+ { NameOfCoef un4; EntityType FacetsOf ; NameOfConstraint BlochX; }
+ { NameOfCoef un5; EntityType EdgesOf ; NameOfConstraint BlochX; }
+ { NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint Dirichlet; }
+ { NameOfCoef un4; EntityType FacetsOf ; NameOfConstraint Dirichlet; }
+ { NameOfCoef un5; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+ EndIf
+ }
+ }
+
+ { Name Hcurl_aD; Type Form1;
+ 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)
+ { 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]; }
+ // BF_Edge_3F_a, BF_Edge_3F_b, BF_Edge_4E
+ EndIf
+ }
+ }
+
+ { Name Hcurl_1; Type Form1;
+ 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)
+ { 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]; }
+ // BF_Edge_3F_a, BF_Edge_3F_b, BF_Edge_4E
+ EndIf
+ }
+ }
+ { Name Hcurl_2; Type Form1;
+ 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)
+ { 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]; }
+ // BF_Edge_3F_a, BF_Edge_3F_b, BF_Edge_4E
+ EndIf
+ }
+ Constraint {
+ { NameOfCoef un; EntityType EdgesOf ; NameOfConstraint BlochX; }
+ { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint BlochX; }
+ { NameOfCoef un; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+ { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+ If(flag_o2==1)
+ { NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint BlochX; }
+ { NameOfCoef un4; EntityType FacetsOf ; NameOfConstraint BlochX; }
+ { NameOfCoef un5; EntityType EdgesOf ; NameOfConstraint BlochX; }
+ { NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint Dirichlet; }
+ { NameOfCoef un4; EntityType FacetsOf ; NameOfConstraint Dirichlet; }
+ { NameOfCoef un5; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+ EndIf
+ }
+ }
+ { Name Hcurl_l; Type Form1;
+ 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)
+ // { Name sn3; NameOfCoef un3; Function BF_Edge_3F_a; Support Region[Bound]; Entity FacetsOf[All]; }
+ // { Name sn4; NameOfCoef un4; Function BF_Edge_3F_b; Support Region[Bound]; Entity FacetsOf[All]; }
+ { Name sn5; NameOfCoef un5; Function BF_Edge_4E ; Support Region[Bound]; Entity EdgesOf[All]; }
+ EndIf
+ }
+ Constraint {
+ If(flag_o2==1)
+
+ EndIf
+ }
+ }
+
+
+}
+
+Formulation {
+ { Name modal_Aux_E; Type FemEquation;
+ Quantity {
+ { Name e ; Type Local; NameOfSpace Hcurl ;}
+ { Name eaD; Type Local; NameOfSpace Hcurl_aD;}
+ }
+ Equation {
+ Galerkin{ [ 1.* cel^2/mur[] * Dof{Curl e}, {Curl e} ] ; In Om ; Jacobian JSur; Integration Int_1;}
+ Galerkin{ Eig[ -1 *-epsr_nod[] * Dof{e} , {e} ] ; Order 2 ; In Om ; Jacobian JSur; Integration Int_1;}
+ Galerkin{ [ 1 * 2.*om_d_1/Sqrt[2] * Dof{e} , {eaD} ] ; In Om_1 ; Jacobian JSur; Integration Int_1;}
+ Galerkin{ [ -I[] * gam_1 * Dof{eaD} , {eaD} ] ; In Om_1 ; Jacobian JSur; Integration Int_1;}
+ Galerkin{ Eig[ -I[] * Om_D[] * Dof{eaD } , {e} ] ; Order 1 ; In Om_1 ; Jacobian JSur; Integration Int_1;}
+ Galerkin{ Eig[ -I[] *-Dof{eaD} , {eaD} ] ; Order 1 ; In Om_1 ; Jacobian JSur; Integration Int_1;}
+ }
+ }
+ { Name modal_helmholtz_Lag_E; Type FemEquation;
+ Quantity {
+ { Name u_1 ; Type Local ; NameOfSpace Hcurl_1 ;}
+ { Name u_2 ; Type Local ; NameOfSpace Hcurl_2 ;}
+ { Name lambda; Type Local ; NameOfSpace Hcurl_l ;}
+ }
+ Equation {
+ Galerkin { [ 1.* term_o3_0rr_1[]/mur[]* Dof{Curl u_1},{Curl u_1}]; In Om_1; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[-I[] * term_o3_1cc_1[] * Dof{u_1} ,{u_1} ]; Order 1; In Om_1; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[-I[] * term_o3_1rr_1[]/mur[]* Dof{Curl u_1},{Curl u_1}]; Order 1; In Om_1; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[ -1.* term_o3_2cc_1[] * Dof{u_1} ,{u_1} ]; Order 2; In Om_1; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[ I[] * term_o3_3cc_1[] * Dof{u_1} ,{u_1} ]; Order 3; In Om_1; Jacobian JSur; Integration Int_1;}
+ Galerkin { [ 1. * term_o3_0rr_2[]/mur[]* Dof{Curl u_2},{Curl u_2}]; In Om_2; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[-I[] * term_o3_1cc_2[] * Dof{u_2} ,{u_2} ]; Order 1;In Om_2; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[-I[] * term_o3_1rr_2[]/mur[]* Dof{Curl u_2},{Curl u_2}]; Order 1;In Om_2; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[ -1 * term_o3_2cc_2[] * Dof{u_2} ,{u_2} ]; Order 2;In Om_2; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[ I[] * term_o3_3cc_2[] * Dof{u_2} ,{u_2} ]; Order 3;In Om_2; Jacobian JSur; Integration Int_1;}
+
+ Galerkin { [ 1. * term_o3_0rr_1[] * 1e8 * Dof{lambda} , {u_1 } ] ; In Bound; Jacobian JLin ; Integration Int_1;}
+ Galerkin { [ 1. * -term_o3_0rr_2[] * 1e8 * Dof{lambda} , {u_2 } ] ; In Bound; Jacobian JLin ; Integration Int_1;}
+ Galerkin {Eig[ -I[] * term_o3_1rr_1[] * 1e8 * Dof{lambda} , {u_1 } ] ; Order 1; In Bound; Jacobian JLin ; Integration Int_1;}
+ Galerkin {Eig[ -I[] * -term_o3_1rr_2[] * 1e8 * Dof{lambda} , {u_2 } ] ; Order 1; In Bound; Jacobian JLin ; Integration Int_1;}
+ Galerkin { [ Dof{u_1} * 1e8 , {lambda} ] ; In Bound; Jacobian JLin ; Integration Int_1;}
+ Galerkin { [ -Dof{u_2} * 1e8 , {lambda} ] ; In Bound; Jacobian JLin ; Integration Int_1;}
+ }
+ }
+ { Name modal_helmholtz_PEP_E; Type FemEquation;
+ Quantity {{ Name u ; Type Local; NameOfSpace Hcurl ;}}
+ Equation {
+ Galerkin { [ 1.* cel^2/mur[]*I[]*gam_1 * Dof{Curl u}, {Curl u}]; In Om ; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[-I[]* cel^2/mur[] * Dof{Curl u}, {Curl u}]; Order 1; In Om ; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[-I[]* om_d_1^2 * Dof{u} , {u} ]; Order 1; In Om_1; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[-1. * I[]*-eps_oo_1*gam_1 * Dof{u} , {u} ]; Order 2; In Om_1; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[I[] * -eps_oo_1 * Dof{u} , {u} ]; Order 3; In Om_1; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[-1. * -epsr_nod[]*I[]*gam_1 * Dof{u} , {u} ]; Order 2; In Om_2; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[I[] * -epsr_nod[] * Dof{u} , {u} ]; Order 3; In Om_2; Jacobian JSur; Integration Int_1;}
+ }
+ }
+ {Name modal_helmholtz_PEP_E_nleig; Type FemEquation;
+ Quantity {{ Name u ; Type Local; NameOfSpace Hcurl ;}}
+ Equation {
+ Galerkin { Eig[ 1/mur[] * Dof{Curl u}, {Curl u} ]; Rational 1; In Om ; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[ -epsr_nod[]/cel^2 * Dof{u}, {u} ]; Rational 2; In Om_1; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[ -epsr_nod[]/cel^2 * Dof{u}, {u} ]; Rational 3; In Om_2; Jacobian JSur; Integration Int_1;}
+ }
+ }
+ {Name modal_helmholtz_PEP_h; Type FemEquation;
+ Quantity {{ Name u ; Type Local; NameOfSpace Hgrad_perp ;}}
+ Equation {
+ Galerkin { [ 1.* 1/epsr_nod[] * -om_d_1^2 * Dof{Curl u}, {Curl u} ]; In Om_2; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[-I[]* 1/epsr_nod[] * I[]*eps_oo_1*gam_1* Dof{Curl u}, {Curl u} ]; Order 1; In Om_2; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[ -1.* 1/epsr_nod[] * eps_oo_1 * Dof{Curl u}, {Curl u} ]; Order 2; In Om_2; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[-I[]* 1/epsr_nod[] *I[]*gam_1 * Dof{Curl u}, {Curl u} ]; Order 1; In Om_1; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[ -1.* 1/epsr_nod[] * Dof{Curl u}, {Curl u} ]; Order 2; In Om_1; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[ -1.* mur[]/cel^2 * om_d_1^2 * Dof{u}, {u} ]; Order 2; In Om; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[I[] * -mur[]/cel^2 * I[]*eps_oo_1*gam_1* Dof{u}, {u} ]; Order 3; In Om; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[ 1.* -mur[]/cel^2 * eps_oo_1 * Dof{u}, {u} ]; Order 4; In Om; Jacobian JSur; Integration Int_1;}
+ }
+ }
+}
+
+Resolution {
+ { Name res_PEP_E;
+ System{
+ { Name M1; NameOfFormulation modal_helmholtz_PEP_E ; Type ComplexValue; }
+ }
+ Operation{
+ GenerateSeparate[M1];EigenSolve[M1,neig,eig_target_re,eig_target_im];
+ // Print[M1];
+ }
+ }
+ { Name res_NEP_E;
+ System{{ Name M1; NameOfFormulation modal_helmholtz_PEP_E_nleig ; Type ComplexValue; }}
+ Operation{
+ GenerateSeparate[M1];
+ // // Fan rational (iw)
+ // (2*Pi*cel/lambda_vp)
+ EigenSolve[M1,neig,eig_target_re,eig_target_im,EigFilter[],
+ { {1}, {-eps_oo_1,gam_1*eps_oo_1, -om_d_1^2,0}, {-1,0,0} } ,
+ { {1}, {1,-gam_1}, {1} } ];
+ }
+ }
+
+ { Name res_Lag_E;
+ System{
+ { Name M1; NameOfFormulation modal_helmholtz_Lag_E ; Type ComplexValue; }
+ }
+ Operation{
+ GenerateSeparate[M1]; EigenSolve[M1,neig,eig_target_re,eig_target_im]; // Print[M1];
+ }
+ }
+ { Name res_PEP_h;
+ System{
+ { Name M1; NameOfFormulation modal_helmholtz_PEP_h ; Type ComplexValue; }
+ }
+ Operation{
+ GenerateSeparate[M1];
+ EigenSolve[M1,neig,eig_target_re,eig_target_im]; // Print[M1];
+ }
+ }
+
+ { Name res_Aux_E;
+ System{
+ { Name M1; NameOfFormulation modal_Aux_E ; Type ComplexValue; }
+ }
+ Operation{
+ GenerateSeparate[M1]; EigenSolve[M1,neig,eig_target_re,eig_target_im,EigFilter[]];
+ }
+ }
+}
+
+PostProcessing {
+ { Name postpro_eigenvectors_PEP_E; NameOfFormulation modal_helmholtz_PEP_E;
+ Quantity {
+ { Name intnormu_pml ; Value { Integral { [ Norm[{u}] ] ; In Om_pml ; Jacobian JSur; Integration Int_1 ;} } }
+ { Name intnormu_nopml ; Value { Integral { [ Norm[{u}] ] ; In Om_nopml; Jacobian JSur; Integration Int_1 ;} } }
+ { Name u ; Value { Local { [ {u} ] ; In Om; Jacobian JSur; } } }
+
+ { Name EigenValuesReal; Value { Local{ [$EigenvalueReal]; In PrintPoint; Jacobian JSur; } } }
+ { Name EigenValuesImag; Value { Local{ [$EigenvalueImag]; In PrintPoint; Jacobian JSur; } } }
+ }
+ }
+ { Name postpro_eigenvectors_NEP_E; NameOfFormulation modal_helmholtz_PEP_E_nleig;
+ Quantity {
+ { Name u ; Value { Local { [ {u} ] ; In Om; Jacobian JSur; } } }
+ { Name EigenValuesReal; Value { Local{ [$EigenvalueReal]; In PrintPoint; Jacobian JSur; } } }
+ { Name EigenValuesImag; Value { Local{ [$EigenvalueImag]; In PrintPoint; Jacobian JSur; } } }
+ }
+ }
+ { Name postpro_eigenvectors_Lag_E; NameOfFormulation modal_helmholtz_Lag_E;
+ Quantity {
+ { Name u_1 ; Value { Local { [ {u_1} ] ; In Om_1; Jacobian JSur; } } }
+ { Name u_2 ; Value { Local { [ {u_2} ] ; In Om_2; Jacobian JSur; } } }
+ { Name lambda ; Value { Local { [ {lambda} ] ; In Bound ; Jacobian JLin; } } }
+ { Name EigenValuesReal; Value { Local{ [$EigenvalueReal]; In PrintPoint; Jacobian JSur; } } }
+ { Name EigenValuesImag; Value { Local{ [$EigenvalueImag]; In PrintPoint; Jacobian JSur; } } }
+ }
+ }
+ { Name postpro_eigenvectors_PEP_h; NameOfFormulation modal_helmholtz_PEP_h;
+ Quantity {
+ { Name u ; Value { Local { [ {u} ] ; In Om; Jacobian JSur; } } }
+ { Name EigenValuesReal; Value { Local{ [$EigenvalueReal]; In PrintPoint; Jacobian JSur; } } }
+ { Name EigenValuesImag; Value { Local{ [$EigenvalueImag]; In PrintPoint; Jacobian JSur; } } }
+ }
+ }
+ { Name postpro_eigenvectors_Aux_E; NameOfFormulation modal_Aux_E;
+ Quantity {
+ { Name e ; Value { Local { [ {e} ] ; In Om; Jacobian JSur; } } }
+ { Name eaD ; Value { Local { [ {eaD} ] ; In Om_1; Jacobian JSur; } } }
+ { Name EigenValuesReal; Value { Local{ [$EigenvalueReal]; In PrintPoint; Jacobian JSur; } } }
+ { Name EigenValuesImag; Value { Local{ [$EigenvalueImag]; In PrintPoint; Jacobian JSur; } } }
+ }
+ }
+}
+
+PostOperation {
+ { Name postop_PEP_E; NameOfPostProcessing postpro_eigenvectors_PEP_E ;
+ Operation {
+ Print [EigenValuesReal, OnElementsOf PrintPoint, Format TimeTable, File "EigenValuesReal.dat"];
+ Print [EigenValuesImag, OnElementsOf PrintPoint, Format TimeTable, File "EigenValuesImag.dat"];
+ If (flag_outEigvec==1)
+ Print [ u , OnElementsOf Om,File "u_PEP_E.pos", Name "Electric eigen-field (PEP-E)"];
+ EndIf
+ }
+ }
+ { Name postop_NEP_E; NameOfPostProcessing postpro_eigenvectors_NEP_E ;
+ Operation {
+ Print [EigenValuesReal, OnElementsOf PrintPoint, Format TimeTable, File "EigenValuesReal.dat"];
+ Print [EigenValuesImag, OnElementsOf PrintPoint, Format TimeTable, File "EigenValuesImag.dat"];
+ If (flag_outEigvec==1)
+ Print [ u , OnElementsOf Om,File "u_NEP_E.pos", Name "Electric eigen-field NEP-E"];
+ EndIf
+ }
+ }
+ { Name postop_Lag_E; NameOfPostProcessing postpro_eigenvectors_Lag_E ;
+ Operation {
+ Print [EigenValuesReal, OnElementsOf PrintPoint, Format TimeTable, File "EigenValuesReal.dat"];
+ Print [EigenValuesImag, OnElementsOf PrintPoint, Format TimeTable, File "EigenValuesImag.dat"];
+ If (flag_outEigvec==1)
+ Print [ u_1 , OnElementsOf Om_1 , File "u1_Lag_E.pos", Name "Electric eigen-field 1 (Lag-E)"];
+ Print [ u_2 , OnElementsOf Om_2 , File "u2_Lag_E.pos", Name "Electric eigen-field 2 (Lag-E)"];
+ Print [ lambda , OnElementsOf Bound, File "lam_Lag_E.pos", Name "Lagrange multiplier (Lag-E)"];
+ EndIf
+ }
+ }
+ { Name postop_PEP_h; NameOfPostProcessing postpro_eigenvectors_PEP_h ;
+ Operation {
+ Print [EigenValuesReal, OnElementsOf PrintPoint, Format TimeTable, File "EigenValuesReal.dat"];
+ Print [EigenValuesImag, OnElementsOf PrintPoint, Format TimeTable, File "EigenValuesImag.dat"];
+ If (flag_outEigvec==1)
+ Print [ u , OnElementsOf Om, File "u_PEP_h.pos", Name "Magnetic eigen-field (PEP-h)"];
+ EndIf
+ }
+ }
+ { Name postop_Aux_E; NameOfPostProcessing postpro_eigenvectors_Aux_E ;
+ Operation {
+ Print [EigenValuesReal, OnElementsOf PrintPoint, Format TimeTable, File "EigenValuesReal.dat"];
+ Print [EigenValuesImag, OnElementsOf PrintPoint, Format TimeTable, File "EigenValuesImag.dat"];
+ If (flag_outEigvec==1)
+ Print [ eaD , OnElementsOf Om_1, File "u_Aux_EaD.pos", Name "Auxiliary Field Aux-E"];
+ Print [ e , OnElementsOf Om , File "u_Aux_E.pos", Name "E Aux-E"];
+ EndIf
+ }
+ }
+}
+
+// In the list of changes of PETSc 3.9
+// http://www.mcs.anl.gov/petsc/documentation/changes/39.html
+// MatSolverPackage is replaced with MatSolverType
+//
+// So you have to change in the command line option:
+// -st_pc_factor_mat_solver_package mumps
+// to:
+// -st_pc_factor_mat_solver_type mumps
+
+// fine tuning slecp options
+// -pep_ncv 600 (set size of the Krylov subspace)
+// -pep_mpd 600
+// -pep_toar_locking 0 ("do not take an EV for granted")
+
+
+eig_tol = 1e-6;
+slepc_options_nep_pol = " -nep_max_it 1000 -nep_type nleigs -nep_target_magnitude ";
+// // Warning slepc > 3.8 : -nep_nleigs_rational => -nep_rational
+slepc_options_nep_rat = Sprintf(" -nep_max_it 1000 -nep_type nleigs -nep_nleigs_rational -nep_target_magnitude -nep_tol %.10f -nep_view",eig_tol);
+// slepc_options_nep_rat = Sprintf(" -nep_max_it 1000 -nep_type nleigs -nep_rational -nep_target_magnitude -nep_tol %.10f ",eig_tol);
+
+// // Warning petsc > 3.9 st_pc_factor_mat_solver_package => st_pc_factor_mat_solver_type
+slepc_options_st = "-solve -st_type sinvert -st_ksp_type preonly -st_pc_type lu -st_pc_factor_mat_solver_package mumps -petsc_prealloc 200 -pos";
+// slepc_options_st = "-solve -st_type sinvert -st_ksp_type preonly -st_pc_type lu -st_pc_factor_mat_solver_type mumps -petsc_prealloc 200 -pos";
+
+slepc_options_pep = Sprintf(" -pep_max_it 400 -pep_target_magnitude -pep_tol %.10f ",eig_tol);
+
+If (flag_res==0)
+ which_res="Aux_E";
+ str_slepc_options=StrCat(slepc_options_st,slepc_options_pep);
+EndIf
+
+If (flag_res==1)
+ which_res="PEP_E";
+ str_slepc_options=StrCat(slepc_options_st,slepc_options_pep);
+EndIf
+
+If (flag_res==2)
+ which_res="NEP_E";
+ str_slepc_options=StrCat(slepc_options_nep_rat);
+EndIf
+
+If (flag_res==3)
+ which_res="Lag_E";
+ str_slepc_options=StrCat(slepc_options_st,slepc_options_pep);
+EndIf
+
+If (flag_res==4)
+ which_res="PEP_h";
+ str_slepc_options=StrCat(slepc_options_st,slepc_options_pep);
+EndIf
+
+// redirect converged eigenvalues at runtime to some textfile
+str_slepc_options = StrCat(str_slepc_options," -nep_monitor_all :temp_eigenvalues.txt ");
+
+DefineConstant[
+ R_ = {StrCat("res_",which_res), Name "GetDP/1ResolutionChoices", Visible 1},
+ C_ = {Sprintf(StrCat("-solve -pos -v2 -slepc ",str_slepc_options,slepc_options_rg)), Name "GetDP/9ComputeCommand", Visible 1},
+ P_ = {StrCat("postop_",which_res), Name "GetDP/2PostOperationChoices", Visible 1}];
diff --git a/NonLinearEVP/NonLinearEVP_data.geo b/NonLinearEVP/NonLinearEVP_data.geo
new file mode 100644
index 0000000000000000000000000000000000000000..035af222e3daef2819c881fc5a62b5c4410ac4bc
--- /dev/null
+++ b/NonLinearEVP/NonLinearEVP_data.geo
@@ -0,0 +1,83 @@
+///////////////////////////////////
+//// Author : Guillaume Demesy ////
+//// _data.geo ////
+///////////////////////////////////
+
+nm = 1.;
+epsilon0 = 8.854187817e-3*nm;
+mu0 = 400.*Pi*nm;
+cel = 1.0/(Sqrt[epsilon0 * mu0]);
+deg2rad = Pi/180;
+
+pp0 = "1Geometry/0";
+pp1 = "2Polarization-Bloch/0";
+pp3 = "3Eigenvalue problem parameters/0";
+pp2 = "4Drude Model parameters/0";
+pp4 = "5Discretization/0";
+pp5 = "6Formulations/0";
+pp6 = "7Output/0";
+close_menu = 0;
+colorpp0 = "MintCream";
+colorpp1 = "Ivory";
+colorpp2 = "PaleGoldenRod";
+colorpp3 = "Ivory";
+colorpp4 = "LightSalmon";
+colorpp5 = "Ivory";
+
+DefineConstant[ a_lat = {50 , Name StrCat[pp0 , "1grating period d [nm]"] , Highlight Str[colorpp0] , Closed close_menu} ];
+
+// normalization factor
+norm = a_lat/(2.*Pi*cel);
+
+DefineConstant[
+ d_sq = {0.806 , Name StrCat[pp0 , "2sq [d]"] , Highlight Str[colorpp0] , Closed close_menu} ,
+ space2pml = {1 , Name StrCat[pp0 , "3PML distance to object [d]"] , Highlight Str[colorpp0] , Closed close_menu} ,
+ pmlsize = {5 , Name StrCat[pp0 , "4PML thickness [d]"] , Highlight Str[colorpp0] , Closed close_menu} ,
+
+ flag_Hparallel = {1 , Name StrCat[pp1 , "1polarization case"] , Choices {0="E //",1="H //"} },
+ kx = {0.75 , Name StrCat[pp1 , "2kx [Pi\a]"] , Highlight Str[colorpp1] , Closed close_menu} ,
+
+ eps_oo_1 = {1 , Name StrCat[pp2 , "0eps_oo_1 [ - ]"] , Highlight Str[colorpp2] , Closed close_menu} ,
+ om_d_1 = {1.1 , Name StrCat[pp2 , "1om_d_1 [2Pic\a]"] , Highlight Str[colorpp2] , Closed close_menu} ,
+ gam_1 = {0.05 , Name StrCat[pp2 , "2gam_1 [2Pic\a]"] , Highlight Str[colorpp2] , Closed close_menu} ,
+
+ neig = {1 , Name StrCat[pp3 , "0Number of eigenvalues [int]"] , Highlight Str[colorpp3] , Closed close_menu} ,
+ eig_target_re = {0.0077, Name StrCat[pp3 , "1EV real part target [2Pic\a]"] , Highlight Str[colorpp3] , Closed close_menu} ,
+ eig_target_im = {0.2598, Name StrCat[pp3 , "2EV imag part target [2Pic\a]"] , Highlight Str[colorpp3] , Closed close_menu} ,
+ eig_min_re = {0.007 , Name StrCat[pp3 , "3EV real min [2Pic\a]"] , Highlight Str[colorpp3] , Closed close_menu} ,
+ eig_max_re = {0.009 , Name StrCat[pp3 , "4EV real max [2Pic\a]"] , Highlight Str[colorpp3] , Closed close_menu} ,
+ eig_min_im = {0.25 , Name StrCat[pp3 , "5EV imag min [2Pic\a]"] , Highlight Str[colorpp3] , Closed close_menu} ,
+ eig_max_im = {0.27 , Name StrCat[pp3 , "6EV imag max [2Pic\a]"] , Highlight Str[colorpp3] , Closed close_menu} ,
+
+ paramaille = {4 , Name StrCat[pp4 , "0number of mesh elements per period []"] , Highlight Str[colorpp4] , Closed close_menu} ,
+ 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_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"},
+ flag_outEigvec = {1 , Name StrCat[pp4, "output eigenvector?"], Choices{0,1}}
+];
+
+// Normalization
+d_sq = d_sq * a_lat;
+space2pml = space2pml * a_lat;
+pmlsize = pmlsize * a_lat;
+kx = kx * Pi/a_lat;
+eig_target_re = eig_target_re / norm;
+eig_target_im = eig_target_im / norm;
+eig_min_re = eig_min_re / norm;
+eig_max_re = eig_max_re / norm;
+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;
+
+eps_oo_2 = 1;
+om_d_2 = 0;
+gam_2 = 0;
+
+slepc_options_rg = StrCat(" -rg_interval_endpoints ",
+ Sprintf("%.8lf,",eig_min_re),
+ Sprintf("%.8lf,",eig_max_re),
+ Sprintf("%.8lf,",eig_min_im),
+ Sprintf("%.8lf",eig_max_im));
\ No newline at end of file
diff --git a/NonLinearEVP/README.TXT b/NonLinearEVP/README.TXT
new file mode 100644
index 0000000000000000000000000000000000000000..24cd01a307a68db1f8dae0fb4b35d9dcbca76ad0
--- /dev/null
+++ b/NonLinearEVP/README.TXT
@@ -0,0 +1,16 @@
+A demo of non-linear eigenvalue problems in GetDP.
+
+This model computes one selected eigenfrequency of a grating made of a
+frequency-dispersive material. Five different formulations are given,
+calling the polynomial and (rational) non-linear solvers of the SLEPc library
+thanks to the unified Eig operator.
+
+Quick start
+-----------
+
+Open `NonLinearEVP.pro' with Gmsh.
+
+Additional info
+---------------
+
+Some documentation is available here: https://arxiv.org/abs/1802.02363
\ No newline at end of file