diff --git a/Antennas/dipole.pro b/Antennas/dipole.pro
index cf022c8075069ac03b6c35f819661abbd370d1f4..eb752cce39eb3b3feb2694e59153243862eb6221 100644
--- a/Antennas/dipole.pro
+++ b/Antennas/dipole.pro
@@ -130,8 +130,7 @@ Function {
eta0 = 120*Pi ; // eta0 = Sqrt(mu0/eps0)
- // FIXME: CG "Y[]" should be 0 in 3D !?
- dR[] = (Flag_3Dmodel) ? Unit[ Vector[Sin[Atan2[Z[],X[]]#1], Y[], -Cos[#1]] ] : Vector[0,0,-1] ;
+ dR[] = (Flag_3Dmodel) ? Unit[ Vector[Sin[Atan2[Z[],X[]]#1], 0, Cos[#1]] ] : Vector[0,0,1] ;
V0 = 1. ;
diff --git a/DiffractionGratings/grating2D.pro b/DiffractionGratings/grating2D.pro
index 84d0a5fdb2021ed47436d2aa83edecf85de27467..bbfaf79b97a2a34b0badbd41f0f3b3116c6811bb 100644
--- a/DiffractionGratings/grating2D.pro
+++ b/DiffractionGratings/grating2D.pro
@@ -486,11 +486,11 @@ PostOperation {
Print[ Q_layer_dep[layer_dep] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_layer_dep.txt"]];
Print[ Q_layer_cov[layer_cov] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_layer_cov.txt"]];
Print[ Hz_tot , OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Hz_tot_lambda%.2fnm_1.pos", lambda0/nm)] , Name Sprintf("Hz_tot_%.2fnm.pos", lambda0/nm)];
- // debug
- Print[ u1 , OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("u1%.2fnm_1.pos", lambda0/nm)] , Name Sprintf("u1%.2fnm.pos", lambda0/nm)];
- Print[ E1 , OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("E1%.2fnm_1.pos", lambda0/nm)] , Name Sprintf("E1%.2fnm.pos", lambda0/nm)];
+ // // debug
+ // Print[ u1 , OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("u1%.2fnm_1.pos", lambda0/nm)] , Name Sprintf("u1%.2fnm.pos", lambda0/nm)];
+ // Print[ E1 , OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("E1%.2fnm_1.pos", lambda0/nm)] , Name Sprintf("E1%.2fnm.pos", lambda0/nm)];
Echo[ Str["For i In {0:2}",
- " View[i].LineWidth = 4;View[i].ColormapNumber = 15;View[0].AxesFormatX = '%.1f';",
+ " View[i].LineWidth = 2;View[i].ColormapNumber = 15;View[0].AxesFormatX = '%.1f';",
"EndFor"], File StrCat[myDir,"tmp0.geo"]] ;
// View[i].RangeType = 2;View[i].CustomMin=0.;View[i].CustomMax=1.;
If(multiplot)
diff --git a/DiffractionGratings/grating2D_data_plasmonics.geo b/DiffractionGratings/grating2D_data_plasmonics.geo
index 61f52599e53260ce62bd3965a3ac33e672ac67aa..b6def6905998b05b0d66467a9b78ef0771260782 100644
--- a/DiffractionGratings/grating2D_data_plasmonics.geo
+++ b/DiffractionGratings/grating2D_data_plasmonics.geo
@@ -26,7 +26,7 @@ colorpp02 = "Ivory";
lambda_min = 600;
lambda_max = 800;
-nb_lambdas = 5;
+nb_lambdas = 15;
DefineConstant[
flag_Hparallel = {1 , Name StrCat[pp2, "2polarization case"], Choices {0="E //",1="H //"}, ServerAction "Reset GetDP/T0, GetDP/R0, GetDP/Lambda_step, GetDP/total absorption"},
@@ -89,7 +89,7 @@ DefineConstant[
paramaille_scale_layer_cov = {1 , Name StrCat[pp3, "6Custom mesh parameters/refine cover layer [-]" ] , ReadOnly 0, Highlight Str[colorpp], Closed 1, ServerAction "Reset GetDP/T0, GetDP/R0, GetDP/Lambda_step, GetDP/total absorption"},
paramaille_scale_sup = {1 , Name StrCat[pp3, "6Custom mesh parameters/refine superstrate [-]" ] , ReadOnly 0, Highlight Str[colorpp], Closed 1, ServerAction "Reset GetDP/T0, GetDP/R0, GetDP/Lambda_step, GetDP/total absorption"},
- multiplot = {0, Choices{0,1}, Name StrCat[pp4, "Plot solution on multiple periods"]},
+ multiplot = {1, Choices{0,1}, Name StrCat[pp4, "Plot solution on multiple periods"]},
plotRTgraphs = {1, Choices{0,1}, Name StrCat[pp4, "Plot R and T"]}
] ;
diff --git a/DiffractionGratings/grating3D.geo b/DiffractionGratings/grating3D.geo
index 947806f0281a04d99bc4ad2d4c1c2f475e6c7d75..73cf7ebbb85e3175efe6f3a375a1a9be696afc3d 100644
--- a/DiffractionGratings/grating3D.geo
+++ b/DiffractionGratings/grating3D.geo
@@ -1,21 +1,21 @@
// Copyright (C) 2019 Guillaume Demésy
//
// This file is part of the model grating3D.pro.
-//
+//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
-//
+//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
-//
+//
// You should have received a copy of the GNU General Public License
// along with This program. If not, see <https://www.gnu.org/licenses/>.
-Include StrCat["grating3D_data_",test_case,".geo"];
+Include "grating3D_data.geo";
SetFactory("OpenCASCADE");
e = 5*nm;
@@ -119,7 +119,7 @@ If (tag_geom==6)
Curve Loop(30) = {31, -23, 32, 4, 10};
Surface(18) = {30};
Surface Loop(3) = {5, 16, 18, 12, 11, 9, 8, 7, 10};
- Volume(1) = {3};
+ Volume(1) = {3};
Surface Loop(2) = {6, 15, 14, 10, 7, 8, 9, 13, 17};
Volume(2) = {2};
EndIf
@@ -288,4 +288,4 @@ If (tag_geom==3) // Split torus weird otherwise
EndIf
// Mesh.SurfaceEdges = 0;
Mesh.VolumeEdges = 0;
-Mesh.ElementOrder = og;
\ No newline at end of file
+Mesh.ElementOrder = og;
diff --git a/DiffractionGratings/grating3D.pro b/DiffractionGratings/grating3D.pro
index a7871b49b81a4ac1e4b637eb1b9a3ea8f3642258..1bc873bd94b2a0dc9929046ec2871f5b863e3d75 100644
--- a/DiffractionGratings/grating3D.pro
+++ b/DiffractionGratings/grating3D.pro
@@ -1,22 +1,22 @@
// Copyright (C) 2019 Guillaume Demésy
//
// This file is part of the model grating3D.pro.
-//
+//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
-//
+//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
-//
+//
// You should have received a copy of the GNU General Public License
// along with This program. If not, see <https://www.gnu.org/licenses/>.
-Include StrCat["grating3D_data_",test_case,".geo"];
-Include "grating3D_materials.pro"
+Include "grating3D_data.geo";
+Include "grating3D_materials.pro";
myDir = "res3D/";
@@ -40,7 +40,7 @@ Group {
SurfIntTop = Region[301];
SurfIntBot = Region[302];
-
+
SurfDirichlet = Region[{401,402}];
SurfBloch = Region[{SurfBlochXm,SurfBlochXp,SurfBlochYm,SurfBlochYp}];
@@ -49,7 +49,7 @@ Group {
Else
SurfExcludeFacets = Region[{SurfBloch}];
EndIf
-
+
L_1 = Region[{L_1_temp,SurfIntTop}];
L_6 = Region[{L_6_temp,SurfIntBot}];
@@ -59,7 +59,7 @@ Group {
Omega_super = Region[{Scat,L_1,L_2,L_3,L_4,L_5,PMLtop}];
Omega_subs = Region[{L_6,PMLbot}];
Omega_plot = Region[{L_6,L_5,L_4,L_3,L_2,L_1,Scat}];
-
+
// SurfNeumann = Region[{SurfBlochXm,SurfBlochXp,SurfBlochYm,SurfBlochYp}];
}
@@ -135,7 +135,7 @@ Function{
sz_bermutop[] = Complex[1,Damp_pml_top[]/k1norm[]];
sz_bermubot[] = Complex[1,Damp_pml_bot[]/k2norm[]];
sy = 1.;
-
+
epsr[PMLtop] = Re[epsr1[]]*TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
mur[PMLtop] = TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
epsr[PMLbot] = Re[epsr2[]]*TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
@@ -143,7 +143,7 @@ Function{
// epsr[PMLbot] = epsr2[];
// mur[PMLbot] = TensorDiag[1,1,1];
-
+
epsr_annex[PMLbot] = epsr[];
epsr_annex[PMLtop] = epsr[];
epsr_annex[Omega_source] = epsr1[] * TensorDiag[1,1,1];
@@ -201,7 +201,7 @@ Function{
H1d[Omega_subs] = Ht[];
source[] = (om0/cel)^2*(epsr[]-epsr_annex[])*E1[];
-
+
// Bloch phase shifts
skx1[] = k1x[];
// sky1[] = -k0*n1[]*Sin[theta0]*Sin[phi0+xsi];
@@ -250,7 +250,7 @@ Constraint {
Jacobian {
{ Name JVol ;
- Case {
+ Case {
{ Region All ; Jacobian Vol ; }
}
}
@@ -268,9 +268,9 @@ Jacobian {
Integration {
{ Name I1 ;
- Case {
+ Case {
{ Type Gauss ;
- Case {
+ Case {
{ GeoElement Point ; NumberOfPoints 1 ; }
{ GeoElement Line ; NumberOfPoints 4 ; }
{ GeoElement Triangle ; NumberOfPoints 12 ; }
@@ -334,18 +334,18 @@ Resolution {
System {
{ Name M; NameOfFormulation helmholtz_vector; Type ComplexValue; }
}
- Operation {
+ Operation {
CreateDir[Str[myDir]];
- Generate[M];
- Solve[M]; //SaveSolution[M];
+ Generate[M];
+ Solve[M]; //SaveSolution[M];
}
}
}
-PostProcessing {
+PostProcessing {
{ Name postpro_helmholtz_vector; NameOfFormulation helmholtz_vector; NameOfSystem M;
Quantity {
- { Name u ; Value { Local { [ {u} ]; In Omega; Jacobian JVol; } } }
+ { Name u ; Value { Local { [ {u} ]; In Omega; Jacobian JVol; } } }
{ Name Etot ; Value { Local { [ {u}+E1[] ]; In Omega; Jacobian JVol; } } }
{ Name Htot ; Value { Local { [ H1[]-I[]/(mur[]*mu0*om0)*{Curl u}]; In Omega; Jacobian JVol; } } }
{ Name Edif ; Value { Local { [ {u}+E1d[] ]; In Omega; Jacobian JVol; } } }
@@ -511,4 +511,4 @@ DefineConstant[
// related to? https://lists.mcs.anl.gov/mailman/htdig/petsc-users/2018-October/036377.html
// N.B. pastix works with petsc 3.12
-// -pc_factor_mat_solver_type pastix
\ No newline at end of file
+// -pc_factor_mat_solver_type pastix
diff --git a/DiffractionGratings/grating3D_data.geo b/DiffractionGratings/grating3D_data.geo
new file mode 100644
index 0000000000000000000000000000000000000000..d08505222fd57b476b19ca9d653d915070b2957b
--- /dev/null
+++ b/DiffractionGratings/grating3D_data.geo
@@ -0,0 +1,10 @@
+
+DefineConstant[
+ test_case = {"halfellipsoid",
+ Name "1Geometry",
+ Choices {"hole", "pyramid", "torus", "2Dlamellar", "solarcell", "conv", "skew"},
+ GmshOption "Reset", Autocheck 0
+ }
+];
+
+Include StrCat["grating3D_data_",test_case,".geo"];
diff --git a/ElectromagneticScattering/scattering.pro b/ElectromagneticScattering/scattering.pro
index 8768a0c1d6ba0932935c3c89b4f333b9a693ae0d..b4dde40106c78d868bf25f980451d950f0d1da2c 100644
--- a/ElectromagneticScattering/scattering.pro
+++ b/ElectromagneticScattering/scattering.pro
@@ -460,6 +460,7 @@ PostProcessing {
{ Name PW_postpro ; NameOfFormulation PW_helmholtz_vector; NameOfSystem P;
Quantity {
{ Name E_tot ; Value { Local { [{u}+Einc[]]; In All_domains; Jacobian JVol; } } }
+ { Name normE_tot ; Value { Local { [Norm[{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[],
@@ -566,13 +567,13 @@ PostOperation {
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"]];
+ Print [ E_tot , OnElementsOf Domain , File StrCat[myDir,"E_tot.pos"]];
+ Print [ E_tot , OnElementsOf Domain , File StrCat[myDir,"E_tot.pos"]];
+ Print [ normE_tot, OnElementsOf Scat_In , File StrCat[myDir,"normE_tot.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]}
diff --git a/ElectromagneticScattering/scattering_data.geo b/ElectromagneticScattering/scattering_data.geo
index 91e3c08c3ced34e22ad336208c257be16360a813..8cea03ed354067b62804b7b297c237710caf08fd 100644
--- a/ElectromagneticScattering/scattering_data.geo
+++ b/ElectromagneticScattering/scattering_data.geo
@@ -10,7 +10,7 @@ cel = 1.0/(Sqrt[epsilon0 * mu0]);
deg2rad = Pi/180.;
pp0 = "1Geometry/0";
-pp1 = "2Study Type/0";
+pp1 = "0Study Type/0";
pp2 = "3Electromagnetic parameters/0";
pp3 = "4Mesh size and PMLs parameters/0";
pp4 = "5Postpro options/0";
@@ -34,7 +34,7 @@ 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}];
+ Choices {ELL="ellispoid",PARALL="parallelepiped",CYL="cylinder",CONE="cone",TOR="split torus"},Closed 0}];
If (flag_shape==ELL)
DefineConstant[
@@ -103,7 +103,7 @@ DefineConstant[
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}];
+ 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},