diff --git a/libmeepgeom/meepgeom.cpp b/libmeepgeom/meepgeom.cpp index ee7f7acc8..68bb56aaa 100644 --- a/libmeepgeom/meepgeom.cpp +++ b/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 diff --git a/libmeepgeom/user-defined-material.cpp b/libmeepgeom/user-defined-material.cpp index b77306ced..412579474 100644 --- a/libmeepgeom/user-defined-material.cpp +++ b/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 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,6 +171,7 @@ 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