new .geo + OpenCASCADE for photon trap geometry

tests/haroche/cavity.dat

+// Scaling
+mm = 1e-3;
+nm = 1e-9;
+
+// Physics
+C0   = 299792458;           // [m/s]
+Mu0  = 4 * Pi * 1e-7;       // [H/m]
+Eps0 = 1 / (C0 * C0 * Mu0); // [F/m]
+
+// First guess: resonance frequency, wavelength & wavenumber
+DefineConstant[
+  F0      = { 51.099e9,   Name "Input/00FirstGuess/00Frequency",  Units "Hz" },
+  Lambda0 = { C0/F0,      Name "Input/00FirstGuess/01Wavelength", Units "m",
+              ReadOnly 1 },
+  K0      = { 2*Pi*F0/C0, Name "Input/00FirstGuess/02Wavenumber", Units "rad/m",
+              ReadOnly 1}
+];
+
+// Cavity
+DefineConstant[
+  r             = { 39.400*mm, Name "Input/01Cavity/00Small curv.", Units "m" },
+  R             = { 40.600*mm, Name "Input/01Cavity/01Large curv.", Units "m" },
+  L_cav         = { 27.570*mm, Name "Input/01Cavity/02Length",      Units "m" },
+  radius_mirror = { 25.000*mm, Name "Input/01Cavity/03Radius",      Units "m" },
+  thick_mirror  = {  1.415*mm, Name "Input/01Cavity/04Thickness",   Units "m" }
+];
+
+// Air box and PML
+DefineConstant[
+  distXY = { 1, Name "Input/02PML/00Air box size (XY)", Units "wavelength" },
+  distZ  = { 1, Name "Input/02PML/01Air box size (Z)",  Units "wavelength" },
+  thick  = { 1, Name "Input/02PML/03Thickness",         Units "wavelength" }
+];
+
+dist2PML_xy =  distXY * Lambda0;
+dist2PML_z  =  distZ  * Lambda0;
+
+box_x =             radius_mirror + dist2PML_xy;
+box_y =             radius_mirror + dist2PML_xy;
+box_z = L_cav / 2. + thick_mirror + dist2PML_z;
+
+pml_x = thick * Lambda0;
+pml_y = thick * Lambda0;
+pml_z = thick * Lambda0;
+
+// Mesh
+DefineConstant[
+  mAir = { 8, Name "Input/03Mesh/00Air density",    Units "elem./wavelength" },
+  mPML = { 8, Name "Input/03Mesh/01PML density",    Units "elem./wavelength" },
+  mMir = { 8, Name "Input/03Mesh/02Mirror density", Units "elem./wavelength" },
+  mO   = { 2, Name "Input/03Mesh/03Order",          Units "-"                }
+];
+
+clPML = Lambda0 / mPML;
+clAir = Lambda0 / mAir;
+clMir = Lambda0 / mMir;

tests/haroche/cavity.geo

+// OpenCASCADE //
+/////////////////
+SetFactory("OpenCASCADE");
+
+// Data //
+//////////
+Include "cavity.dat";
+
+// Mirror //
+////////////
+// Mirror is made by cutting a portion of a (self-intersecting) torus
+// --> V. Debierre et al., "Absorption in quantum electrodynamic cavities
+//     in terms of a quantum jump operator", Phys. Rev. A 90, 033806, 2014
+//     [DOI: 10.1103/PhysRevA.90.033806]
+Torus(1) = {0,0,-R+L_cav/2, R-r, r};
+Rotate{{1,0,0}, {0,0,-R+L_cav/2}, Pi/2}{Volume{1};}
+
+// This porition is obtained by a cylinder (which also determines
+// the thickness of the mirror)
+Cylinder(2) = {0,0,-R+L_cav/2, 0,0,R+thick_mirror, radius_mirror};
+
+// However, this leads to two connected components (the mirror itself
+// and its "negative"): a sphere is therefor used to keep the mirror only
+Sphere(3) = {0,0,-R+L_cav/2, r};
+
+// Mirror
+m[] = BooleanDifference{Volume{2}; Delete;}{Volume{1, 3}; Delete;};
+
+// Air box //
+/////////////
+v[] = {};
+v  += newv; Box(v[0]) = {0,0,0, box_x, box_y, box_z};
+a[] = BooleanDifference{Volume{v}; Delete;}{Volume{m[0]}; Delete;};
+
+// PML //
+/////////
+v  += newv; Box(v[1]) = {box_x,    0,    0,  +pml_x,+box_y,+box_z}; //   X(+)
+v  += newv; Box(v[2]) = {    0,box_y,    0,  +box_x,+pml_y,+box_z}; //   Y(+)
+v  += newv; Box(v[3]) = {    0,    0,box_z,  +box_x,+box_y,+pml_z}; //   Z(+)
+v  += newv; Box(v[4]) = {box_x,    0,box_z,  +pml_x,+box_y,+pml_z}; //  XZ(++)
+v  += newv; Box(v[5]) = {    0,box_y,box_z,  +box_x,+pml_y,+pml_z}; //  ZY(++)
+v  += newv; Box(v[6]) = {box_x,box_y,    0,  +pml_x,+pml_y,+box_z}; //  YX(++)
+v  += newv; Box(v[7]) = {box_x,box_y,box_z,  +pml_x,+pml_y,+pml_z}; // XYZ(+++)
+g[] = BooleanFragments{Volume{a[], v[]}; Delete;}{};                // Fragments
+
+// Entities //
+//////////////
+// Warning: might depend on OpenCASCADE and/or Gmsh versions...
+AIR[] = g[0];
+
+PML[]    = g[{1:7}];
+PMLx[]   = g[1];
+PMLy[]   = g[2];
+PMLz[]   = g[3];
+PMLxz[]  = g[4];
+PMLzy[]  = g[5];
+PMLyx[]  = g[6];
+PMLxyz[] = g[7];
+
+BND[] = CombinedBoundary{Volume{AIR[], PML[]};};
+OXZ[] = BND[{1,7,13,16}];
+OZY[] = BND[{0,9,12,18}];
+OYX[] = BND[{5,8,11,23}];
+MIR[] = BND[2];
+FRM[] = BND[{3,4}];
+OUT[] = BND[{6,10,14,15,17,19,20,21,22,24,25,26}];
+
+DMY[] = PointsOf{Volume{AIR[]};};
+
+// Mesh //
+//////////
+Characteristic Length{PointsOf{ Volume{PML[]};}} = clPML;
+Characteristic Length{PointsOf{ Volume{AIR[]};}} = clAir;
+Characteristic Length{PointsOf{Surface{MIR[]};}} = clMir;
+
+Mesh.Algorithm         = 6;  // Force Frontal-Delaunay for 2D mesh
+Mesh.Algorithm3D       = 10; // Force HXT for 3D mesh
+Mesh.HighOrderOptimize = 1;
+Mesh.ElementOrder      = mO;
+
+// Physical //
+//////////////
+Physical  Volume(138) = {AIR[]};
+
+Physical  Volume(139) = {PMLx[]};
+Physical  Volume(141) = {PMLy[]};
+Physical  Volume(142) = {PMLz[]};
+Physical  Volume(140) = {PMLyx[]};
+Physical  Volume(144) = {PMLxz[]};
+Physical  Volume(145) = {PMLzy[]};
+Physical  Volume(143) = {PMLxyz[]};
+
+Physical Surface(146) = {OXZ[]};
+Physical Surface(147) = {OZY[]};
+Physical Surface(149) = {OYX[]};
+Physical Surface(150) = {OUT[]};
+
+Physical Surface(148) = {MIR[]};
+Physical Surface(151) = {FRM[]};
+
+Physical Point(1000000) = {DMY[0]};
+
+// Write down PML data //
+/////////////////////////
+Printf("%.16e", pml_x)  > "pml.dat";
+Printf("%.16e", pml_y) >> "pml.dat";
+Printf("%.16e", pml_z) >> "pml.dat";
+Printf("%.16e", box_x) >> "pml.dat";
+Printf("%.16e", box_y) >> "pml.dat";
+Printf("%.16e", box_z) >> "pml.dat";
+Printf("%.16e",    K0) >> "pml.dat";