support susceptibilities in user-defined materials

libmeepgeom/meepgeom.cpp

@@ -76,6 +76,25 @@ static bool medium_struct_equal(const medium_struct *m1, const medium_struct *m2
           susceptibility_list_equal(m1->H_susceptibilities, m2->H_susceptibilities));
 }
 
+// garbage collection for susceptibility_list structures.
+// Assumes that the 'items' field, if non-empty, was allocated using new[];
+// this is automatically the case for python code but is not checked
+// for c++ code and will yield runtime errors if a user's user_material_func
+// uses e.g. malloc() instead.
+static void susceptibility_list_gc(susceptibility_list *sl)
+{ if ( !sl || !(sl->num_items) ) return;
+  delete[] sl->items;
+  sl->num_items=0;
+}
+
+// garbage collection for material structures: called to deallocate memory
+// allocated for susceptibilities in user-defined materials.
+static void material_gc(material_type m)
+ { if ( !m || m->which_subclass!=material_data::MATERIAL_USER ) return;
+   susceptibility_list_gc( &(m->medium.E_susceptibilities) );
+   susceptibility_list_gc( &(m->medium.H_susceptibilities) );
+ }
+
 static bool material_type_equal(const material_type m1, const material_type m2)
 {
     if (m1 == m2) return true;
@@ -636,15 +655,7 @@ void geom_epsilon::get_material_pt(material_type &material, const meep::vec &r)
      case material_data::MATERIAL_USER:
       md->medium = vacuum_medium;
       md->user_func(p, md->user_data, &(md->medium));
-      // TODO: update this to allow user's function to set
-      //       position-dependent susceptibilities. For now
-      //       it's an error if the user's function creates
-      //       any.
-      if (    (md->medium.E_susceptibilities.num_items>0)
-           || (md->medium.H_susceptibilities.num_items>0)
-         )
-       meep::abort("susceptibilities in user-defined-materials not yet supported");
-       return;
+      return;
 
      // position-independent material or metal: there is nothing to do
      case material_data::MEDIUM:
@@ -675,6 +686,7 @@ double geom_epsilon::chi1p1(meep::field_type ft, const meep::vec &r)
   material_type material;
   get_material_pt(material, r);
   material_epsmu(ft, material, &chi1p1, &chi1p1_inv);
+  material_gc(material);
 
   return (chi1p1.m00 + chi1p1.m11 + chi1p1.m22)/3;
 }
@@ -829,6 +841,7 @@ void geom_epsilon::eff_chi1inv_matrix(meep::component c, symmetric_matrix *chi1i
     get_material_pt(mat, v.center());
   trivial:
     material_epsmu(meep::type(c), mat, &meps, chi1inv_matrix);
+    material_gc(mat);
     return;
   }
 
@@ -994,6 +1007,7 @@ void geom_epsilon::fallback_chi1inv_row(meep::component c,
   material_type material;
   get_material_pt(material, v.center());
   material_epsmu(meep::type(c), material, &chi1p1, &chi1p1_inv);
+  material_gc(material);
   if (chi1p1.m01 != 0 || chi1p1.m02 != 0 || chi1p1.m12 != 0
       || chi1p1.m00 != chi1p1.m11 || chi1p1.m11 != chi1p1.m22 ||
       chi1p1.m00 != chi1p1.m22) {
@@ -1138,6 +1152,7 @@ double geom_epsilon::chi(meep::component c, const meep::vec &r, int p) {
        chi_val = 0;
    };
 
+  material_gc(material);
   return chi_val;
 }
 
@@ -1229,6 +1244,7 @@ double geom_epsilon::conductivity(meep::component c, const meep::vec &r) {
   default:
     cond_val = 0;
   }
+  material_gc(material);
 
   // if the user specified scalar absorbing layers, add their conductivities
   // to cond_val (isotropically, for both magnetic and electric conductivity).
@@ -1315,6 +1331,7 @@ void geom_epsilon::sigma_row(meep::component c, double sigrow[3],
 	}
 	break;
       }
+  material_gc(mat);
 }
 
 // add a polarization to the list if it is not already there

libmeepgeom/user-defined-material.cpp

@@ -20,6 +20,28 @@
   #define DATADIR "./"
 #endif
 
+/***************************************************************/
+/* geometric parameters ****************************************/
+/***************************************************************/
+#define R1X 0.5      // ellipsoid minor radius
+#define R1Y 1.0      // ellipsoid major radius
+#define R2  3.0      // cylinder radius
+
+/***************************************************************/
+/* parameters in (two-oscillator) Drude-Lorenz model of Ag     */
+/***************************************************************/
+#define EV_UM (1.0/1.23984193)
+#define AG_WP (9.01 * EV_UM)
+#define AG_F0 0.845
+
+#define AG_FRQ0 1.0e-10
+#define AG_GAM0 0.48*EV_UM
+#define AG_SIG0 (AG_F0*AG_WP*AG_WP)/(AG_FRQ0*AG_FRQ0)
+
+#define AG_FRQ1 0.065
+#define AG_GAM1 0.816*EV_UM
+#define AG_SIG1 (AG_F0*AG_WP*AG_WP)/(AG_FRQ1*AG_FRQ1)
+
 using namespace meep;
 
 typedef std::complex<double> cdouble;
@@ -65,7 +87,6 @@ bool compare_hdf5_datasets(const char *file1, const char *name1,
   normDelta = sqrt(normDelta) / size;
 
   double norm12 = fmax( sqrt(norm1), sqrt(norm2) );
-printf("Norm{1,2,d} = {%e,%e,%e} (%e)\n",norm1,norm2,normDelta,normDelta/norm12);
   if (normDelta > rel_tol*norm12)
    return false;
   return true;
@@ -82,25 +103,62 @@ double dummy_eps(const vec &) { return 1.0; }
 /* material function that recreates the ellipsoid-in-cylinder  */
 /* configuration of the cyl-ellipsoid sample code              */
 /***************************************************************/
+typedef struct my_material_func_data
+  { 
+    double rxInner, ryInner, rOuter;
+    bool with_susceptibility; 
+  } my_material_func_data;
+
 void my_material_func(vector3 p, void *user_data, meep_geom::medium_struct *m)
 { 
-  (void) user_data;
-#define R1X 0.5
-#define R1Y 1.0
-#define R2  3.0
+  my_material_func_data *data=(my_material_func_data *)user_data;
+  double rxInner = data->rxInner, rxInner2=rxInner*rxInner;
+  double ryInner = data->ryInner, ryInner2=ryInner*ryInner;
+  double rOuter  = data->rOuter,  rOuter2 =rOuter*rOuter;
 
-  double x=p.x, y=p.y;
+  double x=p.x, x2=x*x, y=p.y, y2=y*y;
 
   // test for point inside inner ellipsoid
-  double nn;
-  if ( (x*x/(R1X*R1X) + y*y/(R1Y*R1Y)) < 1.0 )
-   nn=1.0;
-  else if ( ( x*x/(R2*R2) + y*y/(R2*R2) ) < 1.0 )
-   nn=3.5;
-  else 
-   nn=1.0;
+  bool innermost = (( x2/rxInner2 + y2/ryInner2 ) < 1.0 );
+  bool outermost = ( (x*x + y*y) > rOuter2 );
+  bool in_middle = (!innermost && !outermost);
 
+  // set permittivity
+  double nn = in_middle ? 3.5 : 1.0;
   m->epsilon_diag.x = m->epsilon_diag.y = m->epsilon_diag.z = nn*nn;
+
+  // add susceptibilities (two-oscillator model for Ag)
+  if (in_middle&& data->with_susceptibility)
+   { 
+      m->E_susceptibilities.num_items = 2;
+      m->E_susceptibilities.items = new meep_geom::susceptibility[2];
+
+      m->E_susceptibilities.items[0].sigma_offdiag.x = 0.0;
+      m->E_susceptibilities.items[0].sigma_offdiag.y = 0.0;
+      m->E_susceptibilities.items[0].sigma_offdiag.z = 0.0;
+      m->E_susceptibilities.items[0].sigma_diag.x    = AG_SIG0;
+      m->E_susceptibilities.items[0].sigma_diag.y    = AG_SIG0;
+      m->E_susceptibilities.items[0].sigma_diag.z    = AG_SIG0;
+      m->E_susceptibilities.items[0].frequency       = AG_FRQ0;
+      m->E_susceptibilities.items[0].gamma           = AG_GAM0;
+      m->E_susceptibilities.items[0].noise_amp       = 0.0;
+      m->E_susceptibilities.items[0].drude           = true;
+      m->E_susceptibilities.items[0].is_file         = false;
+
+      m->E_susceptibilities.items[1].sigma_offdiag.x = 0.0;
+      m->E_susceptibilities.items[1].sigma_offdiag.y = 0.0;
+      m->E_susceptibilities.items[1].sigma_offdiag.z = 0.0;
+      m->E_susceptibilities.items[1].sigma_diag.x    = AG_SIG1;
+      m->E_susceptibilities.items[1].sigma_diag.y    = AG_SIG1;
+      m->E_susceptibilities.items[1].sigma_diag.z    = AG_SIG1;
+      m->E_susceptibilities.items[1].frequency       = AG_FRQ1;
+      m->E_susceptibilities.items[1].gamma           = AG_GAM1;
+      m->E_susceptibilities.items[1].noise_amp       = 0.0;
+      m->E_susceptibilities.items[1].drude           = true;
+      m->E_susceptibilities.items[1].is_file         = false;
+
+   }
+
 }
 
 /***************************************************************/
@@ -113,13 +171,16 @@ int main(int argc, char *argv[])
   // simple argument parsing
   meep::component src_cmpt=Ez;
   std::string eps_ref_file = "cyl-ellipsoid-eps-ref.h5";
+  bool with_susceptibility=true;
   for(int narg=1; narg<argc; narg++)
    {
      if (argv[narg] && !strcmp(argv[narg],"--eps_ref_file"))
       { if (narg+1 == argc)
          abort("no option specified for --eps_ref_file");
         eps_ref_file=argv[++narg];
       }
+     else if (argv[narg] && !strcmp(argv[narg],"--without_susceptibility"))
+      with_susceptibility=false;
      else
       abort("unrecognized command-line option %s",argv[narg]);
    };
@@ -136,8 +197,14 @@ int main(int argc, char *argv[])
                                 : -mirror(X,gv) - mirror(Y,gv);
   structure the_structure(gv, dummy_eps, pml(1.0), sym);
 
+  my_material_func_data data;
+  data.with_susceptibility=with_susceptibility;
+  data.rxInner = R1X;
+  data.ryInner = R1Y;
+  data.rOuter  = R2;
+
   meep_geom::material_type my_material
-   = meep_geom::make_user_material(my_material_func, 0);
+   = meep_geom::make_user_material(my_material_func, (void *)&data);
 
   geometric_object_list g={0, 0};
   bool use_anisotropic_averaging=true;