From 8013355f66f11050568474ec68691c26dd4ea754 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Mon, 30 Oct 2017 09:36:56 -0400 Subject: [PATCH 1/2] fixes for new material_type in libctlgeom 4 --- configure.ac | 2 +- libmeepgeom/array-slice-ll.cpp | 14 +- libmeepgeom/bend-flux-ll.cpp | 26 +-- libmeepgeom/cyl-ellipsoid-ll.cpp | 26 +-- libmeepgeom/material_data.hpp | 19 ++- libmeepgeom/meepgeom.cpp | 264 ++++++++++++++++++------------- libmeepgeom/ring-ll.cpp | 16 +- 7 files changed, 211 insertions(+), 156 deletions(-) diff --git a/configure.ac b/configure.ac index 17cc75a14..af1507c03 100644 --- a/configure.ac +++ b/configure.ac @@ -376,7 +376,7 @@ AC_SUBST(CTL_H_CPPFLAG) save_CPPFLAGS=$CPPFLAGS CPPFLAGS="$CPPFLAGS $CTL_H_CPPFLAG" # Check libctl version >= LIBCTL_MAJOR.LIBCTL_MINOR.LIBCTL_BUGFIX -LIBCTL_MAJOR=3; LIBCTL_MINOR=2; LIBCTL_BUGFIX=0 +LIBCTL_MAJOR=4; LIBCTL_MINOR=0; LIBCTL_BUGFIX=0 AC_MSG_CHECKING([whether libctl version is at least ${LIBCTL_MAJOR}.${LIBCTL_MINOR}.${LIBCTL_BUGFIX}]) AC_EGREP_CPP(yes, [[ #include diff --git a/libmeepgeom/array-slice-ll.cpp b/libmeepgeom/array-slice-ll.cpp index 7a14eedad..14df9a4b1 100644 --- a/libmeepgeom/array-slice-ll.cpp +++ b/libmeepgeom/array-slice-ll.cpp @@ -119,7 +119,7 @@ int main(int argc, char *argv[]) else { master_printf("unknown command-line option %s (aborting)",argv[narg]); usage(argv[0]); - }; + }; }; /***************************************************************/ @@ -142,8 +142,8 @@ int main(int argc, char *argv[]) gv.center_origin(); symmetry sym = use_symmetry ? -mirror(Y,gv) : identity(); structure the_structure(gv, dummy_eps, pml(dpml), sym); - material_type vacuum = meep_geom::vacuum; - material_type dielectric = meep_geom::make_dielectric(eps); + meep_geom::material_type vacuum = meep_geom::vacuum; + meep_geom::material_type dielectric = meep_geom::make_dielectric(eps); geometric_object objects[7]; vector3 origin = v3(0.0, 0.0, 0.0); vector3 xhat = v3(1.0, 0.0, 0.0); @@ -194,7 +194,7 @@ int main(int argc, char *argv[]) #define NX 126 #define NY 38 if (write_files) - { + { h5file *file = f.open_h5file(H5FILENAME); f.output_hdf5(Hz, v1d, file); f.output_hdf5(Sy, v2d, file); @@ -203,7 +203,7 @@ int main(int argc, char *argv[]) exit(0); } else - { + { // // read 1D and 2D array-slice data from HDF5 file // @@ -219,14 +219,14 @@ int main(int argc, char *argv[]) file_slice1d[n] = cdouble(rdata[n], idata[n]); delete[] rdata; delete[] idata; - + file_slice2d = file->read("sy", &rank, dims2D, 2); if (rank!=2 || dims2D[0]!=NX || dims2D[1]!=NY) abort("failed to read 2D reference data from file %s.h5",H5FILENAME); delete file; // - // generate 1D and 2D array slices and compare to + // generate 1D and 2D array slices and compare to // data read from file // rank=f.get_array_slice_dimensions(v1d, dims1D); diff --git a/libmeepgeom/bend-flux-ll.cpp b/libmeepgeom/bend-flux-ll.cpp index e7877eab3..5619c43bc 100644 --- a/libmeepgeom/bend-flux-ll.cpp +++ b/libmeepgeom/bend-flux-ll.cpp @@ -65,15 +65,15 @@ void bend_flux(bool no_bend) // (if no-bend? // (list // (list - // (make block + // (make block // (center wvg-xcen (* 0.5 pad)) // (size w (- sy pad) infinity) // (material (make dielectric (epsilon 12))))))) - // (make block + // (make block // (center 0 wvg-ycen) // (size infinity w infinity) // (material (make dielectric (epsilon 12))))) - // (make block + // (make block // (center (* -0.5 pad) wvg-ycen) // (size (- sx pad) w infinity) // (material (make dielectric (epsilon 12)))) @@ -81,7 +81,7 @@ void bend_flux(bool no_bend) vector3 e2=v3(0.0, 1.0, 0.0); vector3 e3=v3(0.0, 0.0, 1.0); - material_type dielectric = meep_geom::make_dielectric(12.0); + meep_geom::material_type dielectric = meep_geom::make_dielectric(12.0); if (no_bend) { geometric_object objects[1]; @@ -101,13 +101,13 @@ void bend_flux(bool no_bend) e1, e2, e3, v3(sx-pad, w, ENORMOUS) ); - + objects[1] = make_block(dielectric, v3(wvg_xcen, 0.5*pad), e1, e2, e3, v3(w, sy-pad, ENORMOUS) ); - + geometric_object_list g={ 2, objects }; meep_geom::set_materials_from_geometry(&the_structure, g); }; @@ -115,7 +115,7 @@ void bend_flux(bool no_bend) fields f(&the_structure); // (set! sources (list - // (make source + // (make source // (src (make gaussian-src (frequency fcen) (fwidth df))) // (component Ez) // (center (+ 1 (* -0.5 sx)) wvg-ycen) @@ -137,16 +137,16 @@ void bend_flux(bool no_bend) double f_start = fcen-0.5*df; double f_end = fcen+0.5*df; int nfreq = 100; // number of frequencies at which to compute flux - + volume *trans_volume= no_bend ? new volume(vec(0.5*sx-1.5, wvg_ycen), vec(0.0, 2.0*w)) : new volume(vec(wvg_xcen, 0.5*sy-1.5), vec(2.0*w, 0.0)); volume_list trans_vl = volume_list(*trans_volume, Sz); dft_flux trans=f.add_dft_flux(&trans_vl, f_start, f_end, nfreq); - + //(define refl ; reflected flux // (add-flux fcen df nfreq - // (make flux-region + // (make flux-region // (center (+ (* -0.5 sx) 1.5) wvg-ycen) (size 0 (* w 2))))) // volume refl_volume( vec(-0.5*sx+1.5, wvg_ycen), vec(0.0,2.0*w)); @@ -161,7 +161,7 @@ void bend_flux(bool no_bend) refl.scale_dfts(-1.0); }; - //(run-sources+ + //(run-sources+ // (stop-when-fields-decayed 50 Ez // (if no-bend? // (vector3 (- (/ sx 2) 1.5) wvg-ycen) @@ -177,7 +177,7 @@ void bend_flux(bool no_bend) double max_abs=0.0, cur_max=0.0; bool Done=false; do - { + { f.step(); // manually check fields-decayed condition @@ -197,7 +197,7 @@ void bend_flux(bool no_bend) //(if no-bend? (save-flux "refl-flux" refl)) if (no_bend) refl.save_hdf5(f, dataname); - + //(display-fluxes trans refl) printf("%11s | %12s | %12s\n", " Time ", " trans flux", " refl flux"); double f0=fcen-0.5*df, fstep=df/(nfreq-1); diff --git a/libmeepgeom/cyl-ellipsoid-ll.cpp b/libmeepgeom/cyl-ellipsoid-ll.cpp index 53edfed04..6bb0b335a 100644 --- a/libmeepgeom/cyl-ellipsoid-ll.cpp +++ b/libmeepgeom/cyl-ellipsoid-ll.cpp @@ -30,7 +30,7 @@ bool compare_hdf5_datasets(const char *file1, const char *name1, double rel_tol=1.0e-4, double abs_tol=1.0e-8) { h5file f1(file1, h5file::READONLY, false); - int rank1; + int rank1; int *dims1=new int[expected_rank]; double *data1 = f1.read(name1, &rank1, dims1, expected_rank); if (!data1) return false; @@ -75,13 +75,13 @@ int main(int argc, char *argv[]) { initialize mpi(argc, argv); - // simple argument parsing + // simple argument parsing meep::component src_cmpt=Ez; std::string eps_ref_file = "cyl-ellipsoid-eps-ref.h5"; for(int narg=1; narg (*user_material_func)(vector3 x, void *user_data); typedef struct medium_struct { @@ -70,11 +74,6 @@ typedef struct medium_struct { vector3 B_conductivity_diag; } medium_struct; -typedef struct material_type_list - { material_type *items; - int num_items; - } material_type_list; - typedef struct material_data_struct { enum { MATERIAL_TYPE_SELF, MATERIAL_FUNCTION, PERFECT_METAL, MEDIUM } which_subclass; @@ -85,9 +84,15 @@ typedef struct material_data_struct // used only if which_subclass==MEDIUM medium_struct *medium; - + } material_data; + typedef material_data *material_type; + typedef struct material_type_list + { material_type *items; + int num_items; + } material_type_list; + // global variables and exported functions extern material_type vacuum; material_type make_dielectric(double epsilon); diff --git a/libmeepgeom/meepgeom.cpp b/libmeepgeom/meepgeom.cpp index f2d4e3e97..e9c10297e 100644 --- a/libmeepgeom/meepgeom.cpp +++ b/libmeepgeom/meepgeom.cpp @@ -38,10 +38,10 @@ medium_struct vacuum_medium = {0.0, 0.0, 0.0}, /* D_conductivity_diag */ {0.0, 0.0, 0.0} /* B_conductivity_diag */ }; -material_data vacuum_material_data = +material_data vacuum_material_data = { material_data::MEDIUM, 0, 0, &vacuum_medium }; -material_type vacuum = { (void *)&vacuum_material_data }; +material_type vacuum = &vacuum_material_data; /***************************************************************/ /***************************************************************/ @@ -57,14 +57,63 @@ material_type make_dielectric(double epsilon) md->medium->epsilon_diag.x=epsilon; md->medium->epsilon_diag.y=epsilon; md->medium->epsilon_diag.z=epsilon; + return md; +} - material_type mt = { (void *)md }; - return mt; +static void material_type_destroy(material_type m) +{ + (void) m; // unused + // fixme: deallocate m->medium at least +} + +static bool susceptibility_equal(const susceptibility &s1, const susceptibility &s2) +{ + return (vector3_equal(s1.sigma_diag, s2.sigma_diag) && + vector3_equal(s1.sigma_offdiag, s2.sigma_offdiag) && + s1.frequency == s2.frequency && s1.gamma == s2.gamma && + s1.noise_amp == s2.noise_amp && s1.drude == s2.drude && s1.is_self == s2.is_self); +} + +static bool susceptibility_list_equal(const susceptibility_list &s1, const susceptibility_list &s2) +{ + if (s1.num_items != s2.num_items) return false; + for (int i = 0; i < s1.num_items; ++i) + if (!susceptibility_equal(s1.items[i], s2.items[i])) + return false; + return true; +} + +static bool medium_struct_equal(const medium_struct *m1, const medium_struct *m2) +{ + return (vector3_equal(m1->epsilon_diag, m2->epsilon_diag) && + vector3_equal(m1->epsilon_offdiag, m2->epsilon_offdiag) && + vector3_equal(m1->mu_diag, m2->mu_diag) && + vector3_equal(m1->mu_offdiag, m2->mu_offdiag) && + vector3_equal(m1->E_chi2_diag, m2->E_chi2_diag) && + vector3_equal(m1->E_chi3_diag, m2->E_chi3_diag) && + vector3_equal(m1->H_chi2_diag, m2->H_chi2_diag) && + vector3_equal(m1->D_conductivity_diag, m2->D_conductivity_diag) && + vector3_equal(m1->B_conductivity_diag, m2->B_conductivity_diag) && + susceptibility_list_equal(m1->E_susceptibilities, m2->E_susceptibilities) && + susceptibility_list_equal(m1->H_susceptibilities, m2->H_susceptibilities)); +} + +static bool material_type_equal(const material_type m1, const material_type m2) +{ + if (m1 == m2) return true; + if (m1->which_subclass != m2->which_subclass) return false; + switch (m1->which_subclass) { + case material_data::MATERIAL_TYPE_SELF: case material_data::PERFECT_METAL: return true; + case material_data::MATERIAL_FUNCTION: return m1->user_func == m2->user_func && m1->user_data == m2->user_data; + case material_data::MEDIUM: return m1->medium == m2->medium || medium_struct_equal(m1->medium, m2->medium); + default: return false; + } } /***************************************************************/ /***************************************************************/ /***************************************************************/ + typedef struct { double m00, m01, m02, m11, m12, @@ -86,9 +135,9 @@ void sym_matrix_rotate(symmetric_matrix *RAR, A[0][1] = A[1][0] = A_->m01; A[0][2] = A[2][0] = A_->m02; A[1][2] = A[2][1] = A_->m12; - for (i = 0; i < 3; ++i) for (j = 0; j < 3; ++j) + for (i = 0; i < 3; ++i) for (j = 0; j < 3; ++j) AR[i][j] = A[i][0]*R[0][j] + A[i][1]*R[1][j] + A[i][2]*R[2][j]; - for (i = 0; i < 3; ++i) for (j = i; j < 3; ++j) + for (i = 0; i < 3; ++i) for (j = i; j < 3; ++j) A[i][j] = R[0][i]*AR[0][j] + R[1][i]*AR[1][j] + R[2][i]*AR[2][j]; RAR->m00 = A[0][0]; RAR->m11 = A[1][1]; @@ -99,7 +148,7 @@ void sym_matrix_rotate(symmetric_matrix *RAR, } /* Set Vinv = inverse of V, where both V and Vinv are real-symmetric matrices.*/ -void sym_matrix_invert(symmetric_matrix *Vinv, +void sym_matrix_invert(symmetric_matrix *Vinv, const symmetric_matrix *V) { double m00 = V->m00, m11 = V->m11, m22 = V->m22; @@ -114,19 +163,19 @@ void sym_matrix_invert(symmetric_matrix *Vinv, } else { double detinv; - + /* compute the determinant: */ detinv = m00*m11*m22 - m02*m11*m02 + 2.0 * m01*m12*m02 - m01*m01*m22 - m12*m12*m00; - + if (detinv == 0.0) meep::abort( "singular 3x3 matrix"); - + detinv = 1.0/detinv; - + Vinv->m00 = detinv * (m11*m22 - m12*m12); Vinv->m11 = detinv * (m00*m22 - m02*m02); Vinv->m22 = detinv * (m11*m00 - m01*m01); - + Vinv->m02 = detinv * (m01*m12 - m11*m02); Vinv->m01 = detinv * (m12*m02 - m01*m22); Vinv->m12 = detinv * (m01*m02 - m00*m12); @@ -152,7 +201,7 @@ int sym_matrix_positive_definite(symmetric_matrix *V) det2 = m00*m11 - m01*m01; det3 = det2*m22 - m02*m11*m02 + 2.0 * m01*m12*m02 - m12*m12*m00; #endif /* real matrix */ - + return (m00 > 0.0 && det2 > 0.0 && det3 > 0.0); } @@ -175,7 +224,7 @@ void set_dimensions(int dims) vector3 vec_to_vector3(const meep::vec &pt) { vector3 v3; - + switch (pt.dim) { case meep::D1: v3.x = 0; @@ -225,30 +274,30 @@ static geom_box gv2box(const meep::volume &v) return box; } -static bool is_variable(material_type mt) +static bool is_variable(material_type md) { - material_data *md = (material_data *)mt.data; return (md->which_subclass == material_data::MATERIAL_FUNCTION); } +static bool is_variable(void* md) { return is_variable((material_type) md); } -static bool is_self(material_type mt) +static bool is_self(material_type md) { - material_data *md = (material_data *)mt.data; return (md->which_subclass == material_data::MATERIAL_TYPE_SELF); } +static bool is_self(void* md) { return is_self((material_type) md); } -static bool is_medium(material_type mt, medium_struct **m) -{ - material_data *md = (material_data *)mt.data; +static bool is_medium(material_type md, medium_struct **m) +{ if (md->which_subclass == material_data::MEDIUM) { *m = md->medium; return true; }; return false; } +static bool is_medium(void* md, medium_struct **m) { return is_medium((material_type) md, m); } static bool is_metal(meep::field_type ft, const material_type *material) { - material_data *md=(material_data *)material->data; + material_data *md = *material; if (ft == meep::E_stuff) switch (md->which_subclass) { case material_data::MEDIUM: @@ -349,18 +398,18 @@ static meep::realnum linear_interpolate( // return material of the point p from the file (assumed already read) static void epsilon_file_material(material_data *md, vector3 p) { - default_material.data=(void *)md; + default_material = (void*) md; if (!epsilon_data) return; if (md->which_subclass != material_data::MEDIUM) meep::abort("epsilon-input-file only works with a type=medium default-material"); medium_struct *mm=md->medium; - double rx = geometry_lattice.size.x == 0 + double rx = geometry_lattice.size.x == 0 ? 0 : 0.5 + (p.x-geometry_center.x) / geometry_lattice.size.x; - double ry = geometry_lattice.size.y == 0 + double ry = geometry_lattice.size.y == 0 ? 0 : 0.5 + (p.y-geometry_center.y) / geometry_lattice.size.y; - double rz = geometry_lattice.size.z == 0 + double rz = geometry_lattice.size.z == 0 ? 0 : 0.5 + (p.z-geometry_center.z) / geometry_lattice.size.z; - mm->epsilon_diag.x = mm->epsilon_diag.y = mm->epsilon_diag.z = + mm->epsilon_diag.x = mm->epsilon_diag.y = mm->epsilon_diag.z = linear_interpolate(rx, ry, rz, epsilon_data, epsilon_dims[0], epsilon_dims[1], epsilon_dims[2], 1); mm->epsilon_offdiag.x = mm->epsilon_offdiag.y = mm->epsilon_offdiag.z = 0; @@ -383,9 +432,9 @@ class geom_epsilon : public meep::material_function { geometric_object_list geometry; geom_box_tree geometry_tree; geom_box_tree restricted_tree; - + cond_profile cond[5][2]; // [direction][side] - + public: geom_epsilon(geometric_object_list g, material_type_list mlist, const meep::volume &v); @@ -395,7 +444,7 @@ class geom_epsilon : public meep::material_function { double L, double dx, double (*prof)(int,double*,void*), void*, double R); - + virtual void set_volume(const meep::volume &v); virtual void unset_volume(void); @@ -411,7 +460,7 @@ class geom_epsilon : public meep::material_function { virtual double chi1p1(meep::field_type ft, const meep::vec &r); virtual void eff_chi1inv_row(meep::component c, double chi1inv_row[3], - const meep::volume &v, + const meep::volume &v, double tol, int maxeval); void fallback_chi1inv_row(meep::component c, double chi1inv_row[3], @@ -453,32 +502,32 @@ geom_epsilon::geom_epsilon(geometric_object_list g, for (int i = 0; i < geometry.num_items; ++i) { display_geometric_object_info(5, geometry.items[i]); - + medium_struct *mm; if ( is_medium(geometry.items[i].material, &mm) ) printf("%*sdielectric constant epsilon diagonal " "= (%g,%g,%g)\n", 5 + 5, "", - mm->epsilon_diag.x, - mm->epsilon_diag.y, + mm->epsilon_diag.x, + mm->epsilon_diag.y, mm->epsilon_diag.z ); } } - + geom_fix_objects0(geometry); geom_box box = gv2box(v); geometry_tree = create_geom_box_tree0(geometry, box); if (verbose && meep::am_master()) { printf("Geometric-object bounding-box tree:\n"); display_geom_box_tree(5, geometry_tree); - + int tree_depth, tree_nobjects; geom_box_tree_stats(geometry_tree, &tree_depth, &tree_nobjects); master_printf("Geometric object tree has depth %d " "and %d object nodes (vs. %d actual objects)\n", tree_depth, tree_nobjects, geometry.num_items); } - + restricted_tree = geometry_tree; } @@ -491,7 +540,7 @@ geom_epsilon::~geom_epsilon() } void geom_epsilon::set_cond_profile(meep::direction dir, - meep::boundary_side side, + meep::boundary_side side, double L, double dx, double (*P)(int,double*,void*), void *data, double R) @@ -526,11 +575,11 @@ void geom_epsilon::unset_volume(void) void geom_epsilon::set_volume(const meep::volume &v) { unset_volume(); - + geom_box box = gv2box(v); restricted_tree = create_geom_box_tree0(geometry, box); } - + #if 0 #TODO figure this out static material_type eval_material_func(function material_func, vector3 p) @@ -538,13 +587,13 @@ static material_type eval_material_func(function material_func, vector3 p) SCM pscm = ctl_convert_vector3_to_scm(p); material_type material; SCM mo; - + mo = gh_call1(material_func, pscm); material_type_input(mo, &material); - + while (material.which_subclass == MTS::MATERIAL_FUNCTION) { material_type m; - + mo = gh_call1(material.subclass. material_function_data->material_func, pscm); @@ -552,20 +601,20 @@ static material_type eval_material_func(function material_func, vector3 p) material_type_destroy(material); material = m; } - + if (material.which_subclass == MTS::MATERIAL_TYPE_SELF) { epsilon_file_material(material, p); } CK(material.which_subclass != MTS::MATERIAL_FUNCTION, "infinite loop in material functions"); - + return material; } #endif -static void material_epsmu(meep::field_type ft, material_type material, +static void material_epsmu(meep::field_type ft, material_type material, symmetric_matrix *epsmu, symmetric_matrix *epsmu_inv) { - material_data *md=(material_data *)material.data; + material_data *md = material; if (ft == meep::E_stuff) switch (md->which_subclass) { case material_data::MEDIUM: @@ -628,9 +677,9 @@ bool geom_epsilon::get_material_pt(material_type &material, const meep::vec &r) { vector3 p = vec_to_vector3(r); boolean inobject; - material = + material = (material_type) material_of_unshifted_point_in_tree_inobject(p, restricted_tree,&inobject); - material_data *md = (material_data *)(material.data); + material_data *md = material; bool destroy_material = false; if (!inobject && epsilon_data) { @@ -643,7 +692,7 @@ bool geom_epsilon::get_material_pt(material_type &material, const meep::vec &r) destroy_material = true; } else - material = default_material; + material = (material_type) default_material; } else if (md->which_subclass == material_data::MATERIAL_FUNCTION) { // TODO figure this out @@ -672,17 +721,17 @@ double geom_epsilon::chi1p1(meep::field_type ft, const meep::vec &r) material_type material; bool destroy_material = get_material_pt(material, r); - material_epsmu(ft, material, &chi1p1, &chi1p1_inv); - + material_epsmu(ft, material, &chi1p1, &chi1p1_inv); + if (destroy_material) material_type_destroy(material); - + return (chi1p1.m00 + chi1p1.m11 + chi1p1.m22)/3; } /* Find frontmost object in v, along with the constant material behind it. Returns false if material behind the object is not constant. - + Requires moderately horrifying logic to figure things out properly, stolen from MPB. */ static bool get_front_object(const meep::volume &v, @@ -708,7 +757,7 @@ static bool get_front_object(const meep::volume &v, { {0,0,0}, {1,1,1},{1,1,-1},{1,-1,1},{1,-1,-1}, {-1,1,1},{-1,1,-1},{-1,-1,1},{-1,-1,-1} } - }; + }; pixel = gv2box(v); pcenter = p = vec_to_vector3(v.center()); double d1, d2, d3; @@ -727,12 +776,12 @@ static bool get_front_object(const meep::volume &v, if ((id == id1 && vector3_equal(shiftby, shiftby1)) || (id == id2 && vector3_equal(shiftby, shiftby2))) continue; - - mat=default_material; + + mat = (material_type) default_material; if (o) - { material_data *md = (material_data *)o->material.data; + { material_data *md = (material_data *)o->material; if (md->which_subclass != material_data::MATERIAL_TYPE_SELF) - mat=o->material; + mat = md; } if (id1 == -1) { o1 = o; @@ -741,17 +790,17 @@ static bool get_front_object(const meep::volume &v, mat1 = mat; } else if (id2 == -1 || ( (id >= id1 && id >= id2) && - (id1 == id2 - || material_type_equal(&mat1,&mat2)))) { + (id1 == id2 + || material_type_equal(mat1,mat2)))) { o2 = o; shiftby2 = shiftby; id2 = id; mat2 = mat; } else if (!(id1 < id2 && - (id1 == id || material_type_equal(&mat1,&mat))) && + (id1 == id || material_type_equal(mat1,mat))) && !(id2 < id1 && - (id2 == id || material_type_equal(&mat2,&mat)))) + (id2 == id || material_type_equal(mat2,mat)))) return false; } @@ -763,7 +812,7 @@ static bool get_front_object(const meep::volume &v, shiftby2 = shiftby1; } - + if ( (o1 && is_variable(o1->material)) || (o2 && is_variable(o2->material)) || ( (is_variable(default_material) || epsilon_data) @@ -802,7 +851,7 @@ void geom_epsilon::eff_chi1inv_row(meep::component c, double chi1inv_row[3], p, &o, shiftby, mat, mat_behind)) { noavg: destroy_material = get_material_pt(mat, v.center()); - trivial: + trivial: material_epsmu(meep::type(c), mat, &meps, &meps_inv); switch (component_direction(c)) { case meep::X: case meep::R: @@ -835,7 +884,7 @@ void geom_epsilon::eff_chi1inv_row(meep::component c, double chi1inv_row[3], // } /* check for trivial case of only one object/material */ - if (material_type_equal(&mat, &mat_behind)) goto trivial; + if (material_type_equal(mat, mat_behind)) goto trivial; // it doesn't make sense to average metals (electric or magnetic) if (is_metal(meep::type(c), &mat) || is_metal(meep::type(c), &mat_behind)) @@ -848,7 +897,7 @@ void geom_epsilon::eff_chi1inv_row(meep::component c, double chi1inv_row[3], pixel.low = vector3_minus(pixel.low, shiftby); pixel.high = vector3_minus(pixel.high, shiftby); - double fill = box_overlap_with_object(pixel, *o, tol, maxeval); + double fill = box_overlap_with_object(pixel, *o, tol, maxeval); material_epsmu(meep::type(c), mat, &meps, &meps_inv); symmetric_matrix eps2, epsinv2; @@ -918,7 +967,7 @@ void geom_epsilon::eff_chi1inv_row(meep::component c, double chi1inv_row[3], #undef SQR -#define SWAP(a,b) { double xxx = a; a = b; b = xxx; } +#define SWAP(a,b) { double xxx = a; a = b; b = xxx; } /* invert rotation matrix = transpose */ SWAP(Rot[0][1], Rot[1][0]); SWAP(Rot[0][2], Rot[2][0]); @@ -1009,14 +1058,14 @@ void geom_epsilon::fallback_chi1inv_row(meep::component c, if (destroy_material) material_type_destroy(material); if (chi1p1.m01 != 0 || chi1p1.m02 != 0 || chi1p1.m12 != 0 - || chi1p1.m00 != chi1p1.m11 || chi1p1.m11 != chi1p1.m22 || + || chi1p1.m00 != chi1p1.m11 || chi1p1.m11 != chi1p1.m22 || chi1p1.m00 != chi1p1.m22) { int rownum = meep::component_direction(c) % 3; if (rownum == 0) { chi1inv_row[0] = chi1p1.m00; chi1inv_row[1] = chi1p1.m01; chi1inv_row[2] = chi1p1.m02; - } + } else if (rownum == 1) { chi1inv_row[0] = chi1p1.m01; chi1inv_row[1] = chi1p1.m11; @@ -1036,7 +1085,7 @@ void geom_epsilon::fallback_chi1inv_row(meep::component c, vector3 gvmin, gvmax; gvmin = vec_to_vector3(v.get_min_corner()); gvmax = vec_to_vector3(v.get_max_corner()); - xmin[0] = gvmin.x; xmax[0] = gvmax.x; + xmin[0] = gvmin.x; xmax[0] = gvmax.x; if (dim == meep::Dcyl) { xmin[1] = gvmin.z; xmin[2] = gvmin.y; xmax[1] = gvmax.z; xmax[2] = gvmax.y; } @@ -1099,7 +1148,7 @@ bool geom_epsilon::has_chi3(meep::component c) for (int i = 0; i < geometry.num_items; ++i) if ( is_medium(geometry.items[i].material, &mm) ) if (get_chi3(c, mm)!=0) - return true; + return true; for (int i = 0; i < extra_materials.num_items; ++i) if ( is_medium(extra_materials.items[i], &mm )) @@ -1114,10 +1163,10 @@ bool geom_epsilon::has_chi3(meep::component c) double geom_epsilon::chi3(meep::component c, const meep::vec &r) { material_type material; bool destroy_material = get_material_pt(material, r); - + double chi3_val; - - material_data *md=(material_data*)material.data; + + material_data *md = material; switch (md->which_subclass) { case material_data::MEDIUM: chi3_val = get_chi3(c, md->medium); @@ -1125,10 +1174,10 @@ double geom_epsilon::chi3(meep::component c, const meep::vec &r) { default: chi3_val = 0; } - + if (destroy_material) material_type_destroy(material); - + return chi3_val; } @@ -1149,16 +1198,16 @@ bool geom_epsilon::has_chi2(meep::component c) medium_struct *mm; for (int i = 0; i < geometry.num_items; ++i) - if ( is_medium(geometry.items[i].material, &mm) + if ( is_medium(geometry.items[i].material, &mm) && get_chi2(c, mm) != 0 - ) return true; + ) return true; for (int i = 0; i < extra_materials.num_items; ++i) - if ( is_medium(extra_materials.items[i], &mm) + if ( is_medium(extra_materials.items[i], &mm) && get_chi2(c, mm) != 0 ) return true; - return ( is_medium(default_material, &mm) + return ( is_medium(default_material, &mm) && get_chi2(c, mm) != 0 ); } @@ -1166,9 +1215,9 @@ bool geom_epsilon::has_chi2(meep::component c) double geom_epsilon::chi2(meep::component c, const meep::vec &r) { material_type material; bool destroy_material = get_material_pt(material, r); - + double chi2_val; - material_data *md=(material_data *)material.data; + material_data *md = material; switch (md->which_subclass) { case material_data::MEDIUM: chi2_val = get_chi2(c, md->medium); @@ -1176,14 +1225,14 @@ double geom_epsilon::chi2(meep::component c, const meep::vec &r) { default: chi2_val = 0; } - + if (destroy_material) material_type_destroy(material); - + return chi2_val; } -static bool mu_not_1(material_type &m) +static bool mu_not_1(material_type m) { medium_struct *mm; @@ -1197,12 +1246,13 @@ static bool mu_not_1(material_type &m) ) ); } +static bool mu_not_1(void *m) { return mu_not_1((material_type) m); } bool geom_epsilon::has_mu() { for (int i = 0; i < geometry.num_items; ++i) if (mu_not_1(geometry.items[i].material)) - return true; + return true; for (int i = 0; i < extra_materials.num_items; ++i) if (mu_not_1(extra_materials.items[i])) @@ -1238,7 +1288,7 @@ bool geom_epsilon::has_conductivity(meep::component c) for (int i = 0; i < geometry.num_items; ++i) if ( is_medium(geometry.items[i].material,&mm) && get_cnd(c,mm) - ) return true; + ) return true; for (int i = 0; i < extra_materials.num_items; ++i) if ( is_medium(extra_materials.items[i],&mm) @@ -1256,7 +1306,7 @@ double geom_epsilon::conductivity(meep::component c, const meep::vec &r) { bool destroy_material = get_material_pt(material, r); double cond_val; - material_data *md=(material_data *)material.data; + material_data *md = material; switch (md->which_subclass) { case material_data::MEDIUM: cond_val = get_cnd(c, md->medium); @@ -1264,7 +1314,7 @@ double geom_epsilon::conductivity(meep::component c, const meep::vec &r) { default: cond_val = 0; } - + if (destroy_material) material_type_destroy(material); @@ -1299,7 +1349,7 @@ double geom_epsilon::conductivity(meep::component c, const meep::vec &r) { } } } - + return cond_val; } @@ -1318,19 +1368,19 @@ static bool susceptibility_equiv(const susceptibility *o0, return 1; } -void geom_epsilon::sigma_row(meep::component c, double sigrow[3], +void geom_epsilon::sigma_row(meep::component c, double sigrow[3], const meep::vec &r) { vector3 p = vec_to_vector3(r); boolean inobject; - material_type material = + material_type material = (material_type) material_of_unshifted_point_in_tree_inobject(p, restricted_tree, &inobject); - material_data *md=(material_data *)material.data; - + material_data *md = material; + int destroy_material = 0; if (md->which_subclass == material_data::MATERIAL_TYPE_SELF) { - material = default_material; + material = (material_type) default_material; } if (md->which_subclass == material_data::MATERIAL_FUNCTION) { // TODO figure this out @@ -1339,10 +1389,10 @@ void geom_epsilon::sigma_row(meep::component c, double sigrow[3], // p); destroy_material = 1; } - + sigrow[0] = sigrow[1] = sigrow[2] = 0.0; if (md->which_subclass == material_data::MEDIUM) { - susceptibility_list slist = + susceptibility_list slist = type(c) == meep::E_stuff ? md->medium->E_susceptibilities : md->medium->H_susceptibilities; @@ -1369,7 +1419,7 @@ void geom_epsilon::sigma_row(meep::component c, double sigrow[3], break; } } - + if (destroy_material) material_type_destroy(material); } @@ -1399,7 +1449,7 @@ void geom_epsilon::add_susceptibilities(meep::structure *s) { add_susceptibilities(meep::H_stuff, s); } -void geom_epsilon::add_susceptibilities(meep::field_type ft, +void geom_epsilon::add_susceptibilities(meep::field_type ft, meep::structure *s) { pol *pols = 0; medium_struct *mm; @@ -1411,7 +1461,7 @@ void geom_epsilon::add_susceptibilities(meep::field_type ft, ? mm->E_susceptibilities : mm->H_susceptibilities ); - + for (int i = 0; i < extra_materials.num_items; ++i) if ( is_medium(extra_materials.items[i], &mm ) ) pols = add_pols(pols, ft == meep::E_stuff @@ -1424,7 +1474,7 @@ void geom_epsilon::add_susceptibilities(meep::field_type ft, ? mm->E_susceptibilities : mm->H_susceptibilities ); - + for (struct pol *p = pols; p; p = p->next) { susceptibility *ss=&(p->user_s); @@ -1432,7 +1482,7 @@ void geom_epsilon::add_susceptibilities(meep::field_type ft, meep::abort("unknown susceptibility"); bool noisy = (ss->noise_amp!=0.0); - meep::susceptibility *sus = + meep::susceptibility *sus = noisy ? new meep::noisy_lorentzian_susceptibility(ss->noise_amp, ss->frequency, ss->gamma, @@ -1457,7 +1507,7 @@ void geom_epsilon::add_susceptibilities(meep::field_type ft, } } current_pol = NULL; - + while (pols) { struct pol *p = pols; pols = pols->next; @@ -1566,13 +1616,13 @@ void set_materials_from_geometry(meep::structure *s, master_printf("Working in %s dimensions.\n", meep::dimension_name(s->gv.dim)); master_printf("Computational cell is %g x %g x %g with resolution %g\n", - size.x, size.y, size.z, resolution); - + size.x, size.y, size.z, resolution); + material_type_list extra_materials; extra_materials.items=0; extra_materials.num_items=0; geom_epsilon geps(g, extra_materials, gv.pad().surroundings()); - + /***************************************************************/ /***************************************************************/ /***************************************************************/ diff --git a/libmeepgeom/ring-ll.cpp b/libmeepgeom/ring-ll.cpp index bfbbc3deb..af4f39a80 100644 --- a/libmeepgeom/ring-ll.cpp +++ b/libmeepgeom/ring-ll.cpp @@ -52,7 +52,7 @@ int main(int argc, char *argv[]) double dpml=2; // thickness of PML double sxy = 2.0*(r+w+pad+dpml); // cell size - double resolution = 10.0; + double resolution = 10.0; // (set-param! resolution 10) // (set! geometry-lattice (make lattice (size sxy sxy no-size))) @@ -78,13 +78,13 @@ int main(int argc, char *argv[]) // (radius (+ r w)) (material (make dielectric (index n)))) // (make cylinder (center 0 0) (height infinity) // (radius r) (material air)))) - material_type dielectric = meep_geom::make_dielectric(n*n); + meep_geom::material_type dielectric = meep_geom::make_dielectric(n*n); geometric_object objects[2]; vector3 v3zero = {0.0, 0.0, 0.0}; vector3 zaxis = {0.0, 0.0, 1.0}; objects[0] = make_cylinder(dielectric, v3zero, r+w, ENORMOUS, zaxis); objects[1] = make_cylinder(meep_geom::vacuum, v3zero, r, ENORMOUS, zaxis); - geometric_object_list g={ 2, objects }; + geometric_object_list g={ 2, objects }; meep_geom::set_materials_from_geometry(&the_structure, g); fields f(&the_structure); @@ -101,7 +101,7 @@ int main(int argc, char *argv[]) volume v( vec(r+0.1,0.0), vec(0.0, 0.0) ); f.add_volume_source(Ez, src, v); - // (run-sources+ 300 + // (run-sources+ 300 // (at-beginning output-epsilon) // (after-sources (harminv Ez (vector3 (+ r 0.1)) fcen df))) while( f.round_time() < f.last_source_time() ) @@ -115,7 +115,7 @@ int main(int argc, char *argv[]) { f.step(); fieldData.push_back( f.get_field(Ez, eval_pt) ); }; - + #define MAXBANDS 100 cdouble amp[MAXBANDS]; double freq_re[MAXBANDS]; @@ -129,7 +129,7 @@ int main(int argc, char *argv[]) master_printf("--------------------------------------------------------------------------------------------\n"); for(int nb=0; nb 1.0e-2*fabs(ref_freq_re[nb]) || fabs(freq_im[nb]-ref_freq_im[nb]) > 1.0e-2*fabs(ref_freq_im[nb]) || abs( amp[nb]- ref_amp[nb]) > 1.0e-2* abs( ref_amp[nb]) - - ) + + ) abort("harminv band %i disagrees with ref: {re f, im f, re A, im A}={%e,%e,%e,%e}!= {%e,%e,%e,%e}\n", nb, freq_re[nb], freq_im[nb], real(amp[nb]), imag(amp[nb]), ref_freq_re[nb], ref_freq_im[nb], real(ref_amp[nb]), imag(ref_amp[nb])); From eb81b1b06b5f89c81dc3ea3cef12c16414e01875 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Mon, 30 Oct 2017 11:30:30 -0400 Subject: [PATCH 2/2] update python usage of material_type --- python/typemap_utils.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/python/typemap_utils.cpp b/python/typemap_utils.cpp index 31ad71e66..eef35b355 100644 --- a/python/typemap_utils.cpp +++ b/python/typemap_utils.cpp @@ -1,4 +1,4 @@ -/* Copyright (C) 2005-2017 Massachusetts Institute of Technology +/* Copyright (C) 2005-2017 Massachusetts Institute of Technology * * 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 @@ -268,7 +268,7 @@ static int pymaterial_to_material(PyObject *po, material_type *mt) { return 0; } - mt->data = (void *)md; + *mt = md; return 1; } @@ -278,7 +278,7 @@ static int pysphere_to_sphere(PyObject *py_sphere, geometric_object *go) { material_type material; vector3 center; double radius; - + if (!get_attr_v3(py_sphere, ¢er, "center") || !get_attr_dbl(py_sphere, &radius, "radius") || !get_attr_material(py_sphere, &material)) { @@ -318,7 +318,7 @@ static int pywedge_to_wedge(PyObject *py_wedge, geometric_object *wedge) { return 0; } - double wedge_angle; + double wedge_angle; vector3 wedge_start; if (!get_attr_dbl(py_wedge, &wedge_angle, "wedge_angle") || @@ -389,7 +389,7 @@ static int pyellipsoid_to_ellipsoid(PyObject *py_ell, geometric_object *e) { return 0; } - material_type material = blk.material; + material_type material = (material_type) blk.material; vector3 center = blk.center; vector3 e1 = blk.subclass.block_data->e1; vector3 e2 = blk.subclass.block_data->e2;