From ca82d45e6dbc08c049377e5a9adb4bfb2457abc2 Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Sun, 1 Mar 2020 16:22:13 -0700 Subject: [PATCH 01/35] added some TODOs, need more --- src/meep.hpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/meep.hpp b/src/meep.hpp index 4d95e8289..df566312b 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1009,6 +1009,7 @@ class dft_chunk { void operator-=(const dft_chunk &chunk); // the frequencies to loop_in_chunks + // TODO: change here double omega_min, domega; int Nomega; @@ -1103,6 +1104,7 @@ class dft_flux { void remove(); + // TODO: change here double freq_min, dfreq; int Nfreq; dft_chunk *E, *H; @@ -1115,6 +1117,7 @@ class dft_flux { // dft.cpp (normally created with fields::add_dft_energy) class dft_energy { public: + // TODO: change here dft_energy(dft_chunk *E_, dft_chunk *H_, dft_chunk *D_, dft_chunk *B_, double fmin, double fmax, int Nf, const volume &where_); dft_energy(const dft_energy &f); @@ -1140,6 +1143,7 @@ class dft_energy { void remove(); + // TODO: change here double freq_min, dfreq; int Nfreq; dft_chunk *E, *H, *D, *B; @@ -1149,6 +1153,7 @@ class dft_energy { // stress.cpp (normally created with fields::add_dft_force) class dft_force { public: + // TODO: change here dft_force(dft_chunk *offdiag1_, dft_chunk *offdiag2_, dft_chunk *diag_, double fmin, double fmax, int Nf, const volume &where_); dft_force(const dft_force &f); @@ -1167,6 +1172,7 @@ class dft_force { void remove(); + // TODO: change here double freq_min, dfreq; int Nfreq; dft_chunk *offdiag1, *offdiag2, *diag; @@ -1178,6 +1184,7 @@ class dft_near2far { public: /* fourier tranforms of tangential E and H field components in a medium with the given scalar eps and mu */ + // TODO: change here dft_near2far(dft_chunk *F, double fmin, double fmax, int Nf, double eps, double mu, const volume &where_, const direction periodic_d_[2], const int periodic_n_[2], const double periodic_k_[2], const double period_[2]); @@ -1215,6 +1222,7 @@ class dft_near2far { void remove(); + // TODO: change here double freq_min, dfreq; int Nfreq; dft_chunk *F; @@ -1233,6 +1241,7 @@ class dft_near2far { store the Fourier transform per point or per current. */ class dft_ldos { public: + // TODO: change here dft_ldos(double freq_min, double freq_max, int Nfreq); ~dft_ldos() { delete[] Fdft; @@ -1249,6 +1258,7 @@ class dft_ldos { std::complex *Jdft; // Nomega array of J(t) DFT values double Jsum; // sum of |J| over all points public: + // TODO: change here double omega_min, domega; int Nomega; }; From 7a21f0e1848693dcc57fdfd07c8b6ad1039054bd Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Sun, 1 Mar 2020 16:44:49 -0700 Subject: [PATCH 02/35] Should be all TODOs for src/meep.hpp --- src/meep.hpp | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/meep.hpp b/src/meep.hpp index df566312b..fd558ce6d 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1266,12 +1266,14 @@ class dft_ldos { // dft.cpp (normally created with fields::add_dft_fields) class dft_fields { public: + // TODO: change here dft_fields(dft_chunk *chunks, double freq_min, double freq_max, int Nfreq, const volume &where); void scale_dfts(std::complex scale); void remove(); + // TODO: change here double freq_min, dfreq; int Nfreq; dft_chunk *chunks; @@ -1731,27 +1733,36 @@ class fields { double max_abs(derived_component c, const volume &where); // dft.cpp + // TODO: change here dft_chunk *add_dft(component c, const volume &where, double freq_min, double freq_max, int Nfreq, bool include_dV_and_interp_weights = true, std::complex stored_weight = 1.0, dft_chunk *chunk_next = 0, bool sqrt_dV_and_interp_weights = false, std::complex extra_weight = 1.0, bool use_centered_grid = true, int vc = 0); + // TODO: change here dft_chunk *add_dft_pt(component c, const vec &where, double freq_min, double freq_max, int Nfreq); + // TODO: change here dft_chunk *add_dft(const volume_list *where, double freq_min, double freq_max, int Nfreq, bool include_dV = true); void update_dfts(); + // TODO: change here dft_flux add_dft_flux(const volume_list *where, double freq_min, double freq_max, int Nfreq, bool use_symmetry = true); + // TODO: change here dft_flux add_dft_flux(direction d, const volume &where, double freq_min, double freq_max, int Nfreq, bool use_symmetry = true); + // TODO: change here dft_flux add_dft_flux_box(const volume &where, double freq_min, double freq_max, int Nfreq); + // TODO: change here dft_flux add_dft_flux_plane(const volume &where, double freq_min, double freq_max, int Nfreq); // a "mode monitor" is just a dft_flux with symmetry reduction turned off. + // TODO: change here dft_flux add_mode_monitor(direction d, const volume &where, double freq_min, double freq_max, int Nfreq); + // TODO: change here dft_fields add_dft_fields(component *components, int num_components, const volume where, double freq_min, double freq_max, int Nfreq, bool use_centered_grid=true); @@ -1799,12 +1810,15 @@ class fields { void get_mode_mode_overlap(void *mode1_data, void *mode2_data, dft_flux flux, std::complex overlaps[2]); + // TODO: change here dft_energy add_dft_energy(const volume_list *where, double freq_min, double freq_max, int Nfreq); // stress.cpp + // TODO: change here dft_force add_dft_force(const volume_list *where, double freq_min, double freq_max, int Nfreq); // near2far.cpp + // TODO: change here dft_near2far add_dft_near2far(const volume_list *where, double freq_min, double freq_max, int Nfreq, int Nperiods = 1); // monitor.cpp From 0ba7651320179cc2bb38975588c92161266ddf21 Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Sun, 1 Mar 2020 17:36:23 -0700 Subject: [PATCH 03/35] Should have added TODOs for all relevant function/class/struct defintions using dft. Need to hunt down any possible uses of them... --- src/dft.cpp | 43 +++++++++++++++++++++++++++++++++++++++++++ src/meep.hpp | 1 + src/near2far.cpp | 7 +++++++ src/stress.cpp | 12 +++++++++++- 4 files changed, 62 insertions(+), 1 deletion(-) diff --git a/src/dft.cpp b/src/dft.cpp index af5452fd1..744191537 100644 --- a/src/dft.cpp +++ b/src/dft.cpp @@ -33,6 +33,7 @@ namespace meep { struct dft_chunk_data { // for passing to field::loop_in_chunks as void* component c; int vc; + // TODO: change here double omega_min, domega; int Nomega; complex stored_weight, extra_weight; @@ -84,6 +85,7 @@ dft_chunk::dft_chunk(fields_chunk *fc_, ivec is_, ivec ie_, vec s0_, vec s1_, ve sn = sn_; vc = data->vc; + // TODO: change here omega_min = data->omega_min; domega = data->domega; Nomega = data->Nomega; @@ -145,6 +147,7 @@ static void add_dft_chunkloop(fields_chunk *fc, int ichunk, component cgrid, ive shift_phase * S.phase_shift(c, sn), shift, S, sn, chunkloop_data); } +// TODO: change here dft_chunk *fields::add_dft(component c, const volume &where, double freq_min, double freq_max, int Nfreq, bool include_dV_and_interp_weights, complex stored_weight, dft_chunk *chunk_next, @@ -164,6 +167,7 @@ dft_chunk *fields::add_dft(component c, const volume &where, double freq_min, do dft_chunk_data data; data.c = c; data.vc = vc; + // TODO: change here if (Nfreq <= 1) freq_min = freq_max = (freq_min + freq_max) * 0.5; data.omega_min = freq_min * 2 * pi; data.domega = Nfreq <= 1 ? 0.0 : (freq_max * 2 * pi - data.omega_min) / (Nfreq - 1); @@ -182,12 +186,14 @@ dft_chunk *fields::add_dft(component c, const volume &where, double freq_min, do return data.dft_chunks; } +// TODO: change here dft_chunk *fields::add_dft(const volume_list *where, double freq_min, double freq_max, int Nfreq, bool include_dV_and_interp_weights) { dft_chunk *chunks = 0; while (where) { if (is_derived(where->c)) abort("derived_component invalid for dft"); cdouble stored_weight = where->weight; + // TODO: change here chunks = add_dft(component(where->c), where->v, freq_min, freq_max, Nfreq, include_dV_and_interp_weights, stored_weight, chunks); where = where->next; @@ -195,8 +201,10 @@ dft_chunk *fields::add_dft(const volume_list *where, double freq_min, double fre return chunks; } +// TODO: change here dft_chunk *fields::add_dft_pt(component c, const vec &where, double freq_min, double freq_max, int Nfreq) { + // TODO: change here return add_dft(c, where, freq_min, freq_max, Nfreq, false); } @@ -218,6 +226,7 @@ void dft_chunk::update_dft(double time) { if (!fc->f[c][0]) return; for (int i = 0; i < Nomega; ++i) + // TODO: change here, check this whole function dft_phase[i] = polar(1.0, (omega_min + i * domega) * time) * scale; int numcmp = fc->f[c][1] ? 2 : 1; @@ -246,11 +255,13 @@ void dft_chunk::update_dft(double time) { if (numcmp == 2) { complex fc(f[0], f[1]); for (int i = 0; i < Nomega; ++i) + // TODO: check this dft[Nomega * idx_dft + i] += dft_phase[i] * fc; } else { realnum fr = f[0]; for (int i = 0; i < Nomega; ++i) + // TODO: check this dft[Nomega * idx_dft + i] += dft_phase[i] * fr; } idx_dft++; @@ -335,17 +346,20 @@ void load_dft_hdf5(dft_chunk *dft_chunks, component c, h5file *file, const char load_dft_hdf5(dft_chunks, component_name(c), file, dprefix); } +// TODO: change here dft_flux::dft_flux(const component cE_, const component cH_, dft_chunk *E_, dft_chunk *H_, double fmin, double fmax, int Nf, const volume &where_, direction normal_direction_, bool use_symmetry_) : Nfreq(Nf), E(E_), H(H_), cE(cE_), cH(cH_), where(where_), normal_direction(normal_direction_), use_symmetry(use_symmetry_) { + // TODO: change here if (Nf <= 1) fmin = fmax = (fmin + fmax) * 0.5; freq_min = fmin; dfreq = Nf <= 1 ? 0.0 : (fmax - fmin) / (Nf - 1); } dft_flux::dft_flux(const dft_flux &f) : where(f.where) { + // TODO: change here freq_min = f.freq_min; Nfreq = f.Nfreq; dfreq = f.dfreq; @@ -401,9 +415,11 @@ void dft_flux::scale_dfts(complex scale) { if (H) H->scale_dft(scale); } +// TODO: change here dft_flux fields::add_dft_flux(const volume_list *where_, double freq_min, double freq_max, int Nfreq, bool use_symmetry) { if (!where_) // handle empty list of volumes + // TODO: change here return dft_flux(Ex, Hy, NULL, NULL, freq_min, freq_max, Nfreq, v, NO_DIRECTION, use_symmetry); dft_chunk *E = 0, *H = 0; @@ -437,8 +453,10 @@ dft_flux fields::add_dft_flux(const volume_list *where_, double freq_min, double } for (int i = 0; i < 2; ++i) { + // TODO: change here E = add_dft(cE[i], where->v, freq_min, freq_max, Nfreq, true, where->weight * double(1 - 2 * i), E); + // TODO: change here H = add_dft(cH[i], where->v, freq_min, freq_max, Nfreq, false, 1.0, H); } @@ -449,18 +467,23 @@ dft_flux fields::add_dft_flux(const volume_list *where_, double freq_min, double // if the volume list has only one entry, store its component's direction. // if the volume list has > 1 entry, store NO_DIRECTION. direction flux_dir = (where_->next ? NO_DIRECTION : component_direction(where_->c)); + // TODO: change here return dft_flux(cE[0], cH[0], E, H, freq_min, freq_max, Nfreq, firstvol, flux_dir, use_symmetry); } +// TODO: change here dft_energy::dft_energy(dft_chunk *E_, dft_chunk *H_, dft_chunk *D_, dft_chunk *B_, double fmin, double fmax, int Nf, const volume &where_) : Nfreq(Nf), E(E_), H(H_), D(D_), B(B_), where(where_) { + // TODO: change here if (Nf <= 1) fmin = fmax = (fmin + fmax) * 0.5; freq_min = fmin; dfreq = Nf <= 1 ? 0.0 : (fmax - fmin) / (Nf - 1); } +// TODO: change here dft_energy::dft_energy(const dft_energy &f) : where(f.where) { + // TODO: change here freq_min = f.freq_min; Nfreq = f.Nfreq; dfreq = f.dfreq; @@ -511,10 +534,12 @@ double *dft_energy::total() { return F; } +// TODO: change here dft_energy fields::add_dft_energy(const volume_list *where_, double freq_min, double freq_max, int Nfreq) { if (!where_) // handle empty list of volumes + // TODO: change here return dft_energy(NULL, NULL, NULL, NULL, freq_min, freq_max, Nfreq, v); dft_chunk *E = 0, *D = 0, *H = 0, *B = 0; @@ -523,15 +548,20 @@ dft_energy fields::add_dft_energy(const volume_list *where_, double freq_min, do volume_list *where_save = where; while (where) { LOOP_OVER_FIELD_DIRECTIONS(gv.dim, d) { + // TODO: change here E = add_dft(direction_component(Ex, d), where->v, freq_min, freq_max, Nfreq, true, 1.0, E); + // TODO: change here D = add_dft(direction_component(Dx, d), where->v, freq_min, freq_max, Nfreq, false, 1.0, D); + // TODO: change here H = add_dft(direction_component(Hx, d), where->v, freq_min, freq_max, Nfreq, true, 1.0, H); + // TODO: change here B = add_dft(direction_component(Bx, d), where->v, freq_min, freq_max, Nfreq, false, 1.0, B); } where = where->next; } delete where_save; + // TODO: change here return dft_energy(E, H, D, B, freq_min, freq_max, Nfreq, firstvol); } @@ -613,20 +643,25 @@ direction fields::normal_direction(const volume &where) const { return d; } +// TODO: change here dft_flux fields::add_dft_flux(direction d, const volume &where, double freq_min, double freq_max, int Nfreq, bool use_symmetry) { if (d == NO_DIRECTION) d = normal_direction(where); volume_list vl(where, direction_component(Sx, d)); + // TODO: change here dft_flux flux = add_dft_flux(&vl, freq_min, freq_max, Nfreq, use_symmetry); flux.normal_direction = d; return flux; } +// TODO: change here dft_flux fields::add_mode_monitor(direction d, const volume &where, double freq_min, double freq_max, int Nfreq) { + // TODO: change here return add_dft_flux(d, where, freq_min, freq_max, Nfreq, /*use_symmetry=*/false); } +// TODO: change here dft_flux fields::add_dft_flux_box(const volume &where, double freq_min, double freq_max, int Nfreq) { volume_list *faces = 0; @@ -642,20 +677,25 @@ dft_flux fields::add_dft_flux_box(const volume &where, double freq_min, double f } } + // TODO: change here dft_flux flux = add_dft_flux(faces, freq_min, freq_max, Nfreq); delete faces; return flux; } +// TODO: change here dft_flux fields::add_dft_flux_plane(const volume &where, double freq_min, double freq_max, int Nfreq) { + // TODO: change here return add_dft_flux(NO_DIRECTION, where, freq_min, freq_max, Nfreq); } +// TODO: change here dft_fields::dft_fields(dft_chunk *chunks_, double freq_min_, double freq_max_, int Nfreq_, const volume &where_) : where(where_) { chunks = chunks_; + // TODO: change here freq_min = freq_min_; dfreq = Nfreq_ <= 1 ? 0.0 : (freq_max_ - freq_min_) / (Nfreq_ - 1); Nfreq = Nfreq_; @@ -671,6 +711,7 @@ void dft_fields::remove() { } } +// TODO: change here dft_fields fields::add_dft_fields(component *components, int num_components, const volume where, double freq_min, double freq_max, int Nfreq, bool use_centered_grid) { bool include_dV_and_interp_weights = false; @@ -679,10 +720,12 @@ dft_fields fields::add_dft_fields(component *components, int num_components, con cdouble stored_weight = 1.0; dft_chunk *chunks = 0; for (int nc = 0; nc < num_components; nc++) + // TODO: change here chunks = add_dft(components[nc], where, freq_min, freq_max, Nfreq, include_dV_and_interp_weights, stored_weight, chunks, sqrt_dV_and_interp_weights,extra_weight,use_centered_grid); + // TODO: change here return dft_fields(chunks, freq_min, freq_max, Nfreq, where); } diff --git a/src/meep.hpp b/src/meep.hpp index fd558ce6d..0780b27a0 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1082,6 +1082,7 @@ void load_dft_hdf5(dft_chunk *dft_chunks, const char *name, h5file *file, const // dft.cpp (normally created with fields::add_dft_flux) class dft_flux { public: + // TODO: change here dft_flux(const component cE_, const component cH_, dft_chunk *E_, dft_chunk *H_, double fmin, double fmax, int Nf, const volume &where_, direction normal_direction_, bool use_symmetry_); diff --git a/src/near2far.cpp b/src/near2far.cpp index 1628d77d7..540ed388d 100644 --- a/src/near2far.cpp +++ b/src/near2far.cpp @@ -29,11 +29,13 @@ using namespace std; namespace meep { +// TODO: change here dft_near2far::dft_near2far(dft_chunk *F_, double fmin, double fmax, int Nf, double eps_, double mu_, const volume &where_, const direction periodic_d_[2], const int periodic_n_[2], const double periodic_k_[2], const double period_[2]) : Nfreq(Nf), F(F_), eps(eps_), mu(mu_), where(where_) { + // TODO: change here if (Nf <= 1) fmin = fmax = (fmin + fmax) * 0.5; freq_min = fmin; dfreq = Nf <= 1 ? 0.0 : (fmax - fmin) / (Nf - 1); @@ -45,6 +47,7 @@ dft_near2far::dft_near2far(dft_chunk *F_, double fmin, double fmax, int Nf, doub } } +// TODO: change here dft_near2far::dft_near2far(const dft_near2far &f) : freq_min(f.freq_min), dfreq(f.dfreq), Nfreq(f.Nfreq), F(f.F), eps(f.eps), mu(f.mu), where(f.where) { @@ -329,6 +332,7 @@ void dft_near2far::farfield_lowlevel(std::complex *EH, const vec &x) { #endif for (int i = 0; i < Nfreq; ++i) { std::complex EH6[6]; + // TODO: change here??? double freq = freq_min + i * dfreq; size_t idx_dft = 0; LOOP_OVER_IVECS(f->fc->gv, f->is, f->ie, idx) { @@ -523,6 +527,7 @@ double *dft_near2far::flux(direction df, const volume &where, double resolution) static double approxeq(double a, double b) { return fabs(a - b) < 0.5e-11 * (fabs(a) + fabs(b)); } +// TODO: change here dft_near2far fields::add_dft_near2far(const volume_list *where, double freq_min, double freq_max, int Nfreq, int Nperiods) { dft_chunk *F = 0; /* E and H chunks*/ @@ -598,12 +603,14 @@ dft_near2far fields::add_dft_near2far(const volume_list *where, double freq_min, double s = j == 0 ? 1 : -1; /* sign of n x c */ if (is_electric(c)) s = -s; + // TODO: change here F = add_dft(c, w->v, freq_min, freq_max, Nfreq, true, s * w->weight, F, false, 1.0, false, c0); } } } + // TODO: change here return dft_near2far(F, freq_min, freq_max, Nfreq, eps, mu, everywhere, periodic_d, periodic_n, periodic_k, period); } diff --git a/src/stress.cpp b/src/stress.cpp index 4dc05eb76..c478309d6 100644 --- a/src/stress.cpp +++ b/src/stress.cpp @@ -23,10 +23,11 @@ using namespace std; namespace meep { - +// TODO: change here dft_force::dft_force(dft_chunk *offdiag1_, dft_chunk *offdiag2_, dft_chunk *diag_, double fmin, double fmax, int Nf, const volume &where_) : where(where_) { + // TODO: change here if (Nf <= 1) fmin = fmax = (fmin + fmax) * 0.5; freq_min = fmin; Nfreq = Nf; @@ -38,6 +39,7 @@ dft_force::dft_force(dft_chunk *offdiag1_, dft_chunk *offdiag2_, dft_chunk *diag } dft_force::dft_force(const dft_force &f) : where(f.where) { + // TODO: change here freq_min = f.freq_min; Nfreq = f.Nfreq; dfreq = f.dfreq; @@ -132,6 +134,7 @@ void dft_force::scale_dfts(complex scale) { /* note that the components where->c indicate the direction of the force to be computed, so they should be vector components (such as Ex, Ey, ... or Sx, ...) rather than pseudovectors (like Hx, ...). */ +// TODO: change here dft_force fields::add_dft_force(const volume_list *where_, double freq_min, double freq_max, int Nfreq) { dft_chunk *offdiag1 = 0, *offdiag2 = 0, *diag = 0; @@ -148,20 +151,26 @@ dft_force fields::add_dft_force(const volume_list *where_, double freq_min, doub if (coordinate_mismatch(gv.dim, fd)) abort("coordinate-type mismatch in add_dft_force"); if (fd != nd) { // off-diagaonal stress-tensor terms + // TODO: change here offdiag1 = add_dft(direction_component(Ex, fd), where->v, freq_min, freq_max, Nfreq, true, where->weight, offdiag1); + // TODO: change here offdiag2 = add_dft(direction_component(Ex, nd), where->v, freq_min, freq_max, Nfreq, false, 1.0, offdiag2); + // TODO: change here offdiag1 = add_dft(direction_component(Hx, fd), where->v, freq_min, freq_max, Nfreq, true, where->weight, offdiag1); + // TODO: change here offdiag2 = add_dft(direction_component(Hx, nd), where->v, freq_min, freq_max, Nfreq, false, 1.0, offdiag2); } else // diagonal stress-tensor terms LOOP_OVER_FIELD_DIRECTIONS(gv.dim, d) { complex weight1 = where->weight * (d == fd ? +0.5 : -0.5); + // TODO: change here diag = add_dft(direction_component(Ex, d), where->v, freq_min, freq_max, Nfreq, true, 1.0, diag, true, weight1, false); + // TODO: change here diag = add_dft(direction_component(Hx, d), where->v, freq_min, freq_max, Nfreq, true, 1.0, diag, true, weight1, false); } @@ -169,6 +178,7 @@ dft_force fields::add_dft_force(const volume_list *where_, double freq_min, doub } delete where_save; + // TODO: change here return dft_force(offdiag1, offdiag2, diag, freq_min, freq_max, Nfreq, everywhere); } From 354b4606ae0b96abf8a3bb8d8c628ae05a64c1dc Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Mon, 2 Mar 2020 14:03:37 -0700 Subject: [PATCH 04/35] Changed freq_min/freq_max/dfreq to freqs --- src/meep.hpp | 67 +++++++++++++++++++++++++--------------------------- 1 file changed, 32 insertions(+), 35 deletions(-) diff --git a/src/meep.hpp b/src/meep.hpp index 0780b27a0..1fc205685 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1010,8 +1010,8 @@ class dft_chunk { // the frequencies to loop_in_chunks // TODO: change here - double omega_min, domega; int Nomega; + double omegas[Nomega]; component c; // component to DFT (possibly transformed by symmetry) @@ -1083,10 +1083,9 @@ void load_dft_hdf5(dft_chunk *dft_chunks, const char *name, h5file *file, const class dft_flux { public: // TODO: change here - dft_flux(const component cE_, const component cH_, dft_chunk *E_, dft_chunk *H_, double fmin, - double fmax, int Nf, const volume &where_, direction normal_direction_, - bool use_symmetry_); - dft_flux(const dft_flux &f); + dft_flux(const component cE_, const component cH_, dft_chunk *E_, dft_chunk *H_, double *fs, + int Nf, const volume &where_, direction normal_direction_, bool use_symmetry_); + dft_flux(const dft_flux &f); double *flux(); @@ -1106,8 +1105,8 @@ class dft_flux { void remove(); // TODO: change here - double freq_min, dfreq; int Nfreq; + double freqs[Nfreq]; dft_chunk *E, *H; component cE, cH; volume where; @@ -1119,8 +1118,8 @@ class dft_flux { class dft_energy { public: // TODO: change here - dft_energy(dft_chunk *E_, dft_chunk *H_, dft_chunk *D_, dft_chunk *B_, double fmin, double fmax, - int Nf, const volume &where_); + dft_energy(dft_chunk *E_, dft_chunk *H_, dft_chunk *D_, dft_chunk *B_, double *fs, int Nf, + const volume &where_); dft_energy(const dft_energy &f); double *electric(); @@ -1145,8 +1144,8 @@ class dft_energy { void remove(); // TODO: change here - double freq_min, dfreq; int Nfreq; + double freqs[Nfreq]; dft_chunk *E, *H, *D, *B; volume where; }; @@ -1155,8 +1154,8 @@ class dft_energy { class dft_force { public: // TODO: change here - dft_force(dft_chunk *offdiag1_, dft_chunk *offdiag2_, dft_chunk *diag_, double fmin, double fmax, - int Nf, const volume &where_); + dft_force(dft_chunk *offdiag1_, dft_chunk *offdiag2_, dft_chunk *diag_, double *fs, int Nf, + const volume &where_); dft_force(const dft_force &f); double *force(); @@ -1174,8 +1173,8 @@ class dft_force { void remove(); // TODO: change here - double freq_min, dfreq; int Nfreq; + double freqs[Nfreq]; dft_chunk *offdiag1, *offdiag2, *diag; volume where; }; @@ -1186,8 +1185,8 @@ class dft_near2far { /* fourier tranforms of tangential E and H field components in a medium with the given scalar eps and mu */ // TODO: change here - dft_near2far(dft_chunk *F, double fmin, double fmax, int Nf, double eps, double mu, - const volume &where_, const direction periodic_d_[2], const int periodic_n_[2], + dft_near2far(dft_chunk *F, double *fs, int Nf, double eps, double mu, const volume &where_, + const direction periodic_d_[2], const int periodic_n_[2], const double periodic_k_[2], const double period_[2]); dft_near2far(const dft_near2far &f); @@ -1224,8 +1223,8 @@ class dft_near2far { void remove(); // TODO: change here - double freq_min, dfreq; int Nfreq; + double freqs[Nfreq]; dft_chunk *F; double eps, mu; volume where; @@ -1243,7 +1242,7 @@ class dft_near2far { class dft_ldos { public: // TODO: change here - dft_ldos(double freq_min, double freq_max, int Nfreq); + dft_ldos(double *freqs, int Nfreq); ~dft_ldos() { delete[] Fdft; delete[] Jdft; @@ -1260,23 +1259,23 @@ class dft_ldos { double Jsum; // sum of |J| over all points public: // TODO: change here - double omega_min, domega; int Nomega; + double omegas[Nomega]; }; // dft.cpp (normally created with fields::add_dft_fields) class dft_fields { public: // TODO: change here - dft_fields(dft_chunk *chunks, double freq_min, double freq_max, int Nfreq, const volume &where); + dft_fields(dft_chunk *chunks, double *freqs, int Nfreq, const volume &where); void scale_dfts(std::complex scale); void remove(); // TODO: change here - double freq_min, dfreq; int Nfreq; + double freqs[Nfreq]; dft_chunk *chunks; volume where; }; @@ -1735,37 +1734,35 @@ class fields { // dft.cpp // TODO: change here - dft_chunk *add_dft(component c, const volume &where, double freq_min, double freq_max, int Nfreq, + dft_chunk *add_dft(component c, const volume &where, double *freqs, int Nfreq, bool include_dV_and_interp_weights = true, std::complex stored_weight = 1.0, dft_chunk *chunk_next = 0, bool sqrt_dV_and_interp_weights = false, std::complex extra_weight = 1.0, bool use_centered_grid = true, int vc = 0); // TODO: change here - dft_chunk *add_dft_pt(component c, const vec &where, double freq_min, double freq_max, int Nfreq); + dft_chunk *add_dft_pt(component c, const vec &where, double *freqs, int Nfreq); // TODO: change here - dft_chunk *add_dft(const volume_list *where, double freq_min, double freq_max, int Nfreq, - bool include_dV = true); + dft_chunk *add_dft(const volume_list *where, double *freqs, int Nfreq, bool include_dV = true); void update_dfts(); // TODO: change here - dft_flux add_dft_flux(const volume_list *where, double freq_min, double freq_max, int Nfreq, + dft_flux add_dft_flux(const volume_list *where, double *freqs, int Nfreq, bool use_symmetry = true); // TODO: change here - dft_flux add_dft_flux(direction d, const volume &where, double freq_min, double freq_max, - int Nfreq, bool use_symmetry = true); + dft_flux add_dft_flux(direction d, const volume &where, double *freqs, int Nfreq, + bool use_symmetry = true); // TODO: change here - dft_flux add_dft_flux_box(const volume &where, double freq_min, double freq_max, int Nfreq); + dft_flux add_dft_flux_box(const volume &where, double *freqs, int Nfreq); // TODO: change here - dft_flux add_dft_flux_plane(const volume &where, double freq_min, double freq_max, int Nfreq); + dft_flux add_dft_flux_plane(const volume &where, double *freqs, int Nfreq); // a "mode monitor" is just a dft_flux with symmetry reduction turned off. // TODO: change here - dft_flux add_mode_monitor(direction d, const volume &where, double freq_min, double freq_max, - int Nfreq); + dft_flux add_mode_monitor(direction d, const volume &where, double *freqs, int Nfreq); // TODO: change here dft_fields add_dft_fields(component *components, int num_components, const volume where, - double freq_min, double freq_max, int Nfreq, bool use_centered_grid=true); + double *freqs, int Nfreq, bool use_centered_grid=true); /********************************************************/ /* process_dft_component is an intermediate-level */ @@ -1812,16 +1809,16 @@ class fields { std::complex overlaps[2]); // TODO: change here - dft_energy add_dft_energy(const volume_list *where, double freq_min, double freq_max, int Nfreq); + dft_energy add_dft_energy(const volume_list *where, double *freqs, int Nfreq); // stress.cpp // TODO: change here - dft_force add_dft_force(const volume_list *where, double freq_min, double freq_max, int Nfreq); + dft_force add_dft_force(const volume_list *where, double *freqs, int Nfreq); // near2far.cpp // TODO: change here - dft_near2far add_dft_near2far(const volume_list *where, double freq_min, double freq_max, - int Nfreq, int Nperiods = 1); + dft_near2far add_dft_near2far(const volume_list *where, double *freqs, int Nfreq, + int Nperiods = 1); // monitor.cpp std::complex get_chi1inv(component, direction, const vec &loc, double omega = 0, bool parallel = true) const; From 612eb9813a979ae1e0ff4897ea2c8d7d61b188ef Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Mon, 2 Mar 2020 14:18:23 -0700 Subject: [PATCH 05/35] Changed freq_min/freq_max/dfreq to freqs --- src/dft.cpp | 126 ++++++++++++++++++++++++---------------------------- 1 file changed, 57 insertions(+), 69 deletions(-) diff --git a/src/dft.cpp b/src/dft.cpp index 744191537..2ca77041b 100644 --- a/src/dft.cpp +++ b/src/dft.cpp @@ -34,8 +34,8 @@ struct dft_chunk_data { // for passing to field::loop_in_chunks as void* component c; int vc; // TODO: change here - double omega_min, domega; int Nomega; + double omegas[Nomega]; complex stored_weight, extra_weight; double dt_factor; bool include_dV_and_interp_weights; @@ -86,9 +86,8 @@ dft_chunk::dft_chunk(fields_chunk *fc_, ivec is_, ivec ie_, vec s0_, vec s1_, ve vc = data->vc; // TODO: change here - omega_min = data->omega_min; - domega = data->domega; Nomega = data->Nomega; + omegas = data->omegas; dft_phase = new complex[Nomega]; N = 1; @@ -148,11 +147,10 @@ static void add_dft_chunkloop(fields_chunk *fc, int ichunk, component cgrid, ive } // TODO: change here -dft_chunk *fields::add_dft(component c, const volume &where, double freq_min, double freq_max, - int Nfreq, bool include_dV_and_interp_weights, - complex stored_weight, dft_chunk *chunk_next, - bool sqrt_dV_and_interp_weights, complex extra_weight, - bool use_centered_grid, int vc) { +dft_chunk *fields::add_dft(component c, const volume &where, double *freqs, int Nfreq, + bool include_dV_and_interp_weights, complex stored_weight, + dft_chunk *chunk_next, bool sqrt_dV_and_interp_weights, + complex extra_weight, bool use_centered_grid, int vc) { if (coordinate_mismatch(gv.dim, c)) return NULL; /* If you call add_dft before adding sources, it will do nothing @@ -168,10 +166,12 @@ dft_chunk *fields::add_dft(component c, const volume &where, double freq_min, do data.c = c; data.vc = vc; // TODO: change here - if (Nfreq <= 1) freq_min = freq_max = (freq_min + freq_max) * 0.5; - data.omega_min = freq_min * 2 * pi; - data.domega = Nfreq <= 1 ? 0.0 : (freq_max * 2 * pi - data.omega_min) / (Nfreq - 1); + if (Nfreq < 1) abort("Nfreq must be at least 1"); data.Nomega = Nfreq; + double omegas[Nfreq]; + int i; + for (i = 0; i < Nfreq; i++) { omegas[i] = freqs[i] * 2 * pi; } + data.omegas = omegas; data.stored_weight = stored_weight; data.extra_weight = extra_weight; data.dt_factor = dt / sqrt(2.0 * pi); @@ -187,25 +187,24 @@ dft_chunk *fields::add_dft(component c, const volume &where, double freq_min, do } // TODO: change here -dft_chunk *fields::add_dft(const volume_list *where, double freq_min, double freq_max, int Nfreq, +dft_chunk *fields::add_dft(const volume_list *where, double *freqs, int Nfreq, bool include_dV_and_interp_weights) { dft_chunk *chunks = 0; while (where) { if (is_derived(where->c)) abort("derived_component invalid for dft"); cdouble stored_weight = where->weight; // TODO: change here - chunks = add_dft(component(where->c), where->v, freq_min, freq_max, Nfreq, - include_dV_and_interp_weights, stored_weight, chunks); + chunks = add_dft(component(where->c), where->v, freqs, Nfreq, include_dV_and_interp_weights, + stored_weight, chunks); where = where->next; } return chunks; } // TODO: change here -dft_chunk *fields::add_dft_pt(component c, const vec &where, double freq_min, double freq_max, - int Nfreq) { +dft_chunk *fields::add_dft_pt(component c, const vec &where, double *freqs, int Nfreq) { // TODO: change here - return add_dft(c, where, freq_min, freq_max, Nfreq, false); + return add_dft(c, where, freqs, Nfreq, false); } void fields::update_dfts() { @@ -227,7 +226,7 @@ void dft_chunk::update_dft(double time) { for (int i = 0; i < Nomega; ++i) // TODO: change here, check this whole function - dft_phase[i] = polar(1.0, (omega_min + i * domega) * time) * scale; + dft_phase[i] = polar(1.0, omegas[i] * time) * scale; int numcmp = fc->f[c][1] ? 2 : 1; @@ -348,21 +347,19 @@ void load_dft_hdf5(dft_chunk *dft_chunks, component c, h5file *file, const char // TODO: change here dft_flux::dft_flux(const component cE_, const component cH_, dft_chunk *E_, dft_chunk *H_, - double fmin, double fmax, int Nf, const volume &where_, - direction normal_direction_, bool use_symmetry_) + double *fs, int Nf, const volume &where_, direction normal_direction_, + bool use_symmetry_) : Nfreq(Nf), E(E_), H(H_), cE(cE_), cH(cH_), where(where_), normal_direction(normal_direction_), use_symmetry(use_symmetry_) { // TODO: change here - if (Nf <= 1) fmin = fmax = (fmin + fmax) * 0.5; - freq_min = fmin; - dfreq = Nf <= 1 ? 0.0 : (fmax - fmin) / (Nf - 1); + if (Nf < 1) abort("Nf must be at least 1"); + freqs = fs; } dft_flux::dft_flux(const dft_flux &f) : where(f.where) { // TODO: change here - freq_min = f.freq_min; Nfreq = f.Nfreq; - dfreq = f.dfreq; + freqs = f.freqs; E = f.E; H = f.H; cE = f.cE; @@ -416,11 +413,11 @@ void dft_flux::scale_dfts(complex scale) { } // TODO: change here -dft_flux fields::add_dft_flux(const volume_list *where_, double freq_min, double freq_max, - int Nfreq, bool use_symmetry) { +dft_flux fields::add_dft_flux(const volume_list *where_, double *freqs, int Nfreq, + bool use_symmetry) { if (!where_) // handle empty list of volumes // TODO: change here - return dft_flux(Ex, Hy, NULL, NULL, freq_min, freq_max, Nfreq, v, NO_DIRECTION, use_symmetry); + return dft_flux(Ex, Hy, NULL, NULL, freqs, Nfreq, v, NO_DIRECTION, use_symmetry); dft_chunk *E = 0, *H = 0; component cE[2] = {Ex, Ey}, cH[2] = {Hy, Hx}; @@ -454,10 +451,9 @@ dft_flux fields::add_dft_flux(const volume_list *where_, double freq_min, double for (int i = 0; i < 2; ++i) { // TODO: change here - E = add_dft(cE[i], where->v, freq_min, freq_max, Nfreq, true, - where->weight * double(1 - 2 * i), E); + E = add_dft(cE[i], where->v, freqs, Nfreq, true, where->weight * double(1 - 2 * i), E); // TODO: change here - H = add_dft(cH[i], where->v, freq_min, freq_max, Nfreq, false, 1.0, H); + H = add_dft(cH[i], where->v, freqs, Nfreq, false, 1.0, H); } where = where->next; @@ -468,25 +464,23 @@ dft_flux fields::add_dft_flux(const volume_list *where_, double freq_min, double // if the volume list has > 1 entry, store NO_DIRECTION. direction flux_dir = (where_->next ? NO_DIRECTION : component_direction(where_->c)); // TODO: change here - return dft_flux(cE[0], cH[0], E, H, freq_min, freq_max, Nfreq, firstvol, flux_dir, use_symmetry); + return dft_flux(cE[0], cH[0], E, H, freqs, Nfreq, firstvol, flux_dir, use_symmetry); } // TODO: change here -dft_energy::dft_energy(dft_chunk *E_, dft_chunk *H_, dft_chunk *D_, dft_chunk *B_, double fmin, - double fmax, int Nf, const volume &where_) +dft_energy::dft_energy(dft_chunk *E_, dft_chunk *H_, dft_chunk *D_, dft_chunk *B_, double *fs, + int Nf, const volume &where_) : Nfreq(Nf), E(E_), H(H_), D(D_), B(B_), where(where_) { - // TODO: change here - if (Nf <= 1) fmin = fmax = (fmin + fmax) * 0.5; - freq_min = fmin; - dfreq = Nf <= 1 ? 0.0 : (fmax - fmin) / (Nf - 1); + // TODO: change here + if (Nf < 1) abort("Nf must be at least 1"); + freqs = fs; } // TODO: change here dft_energy::dft_energy(const dft_energy &f) : where(f.where) { // TODO: change here - freq_min = f.freq_min; Nfreq = f.Nfreq; - dfreq = f.dfreq; + freqs = f.freqs; E = f.E; H = f.H; D = f.D; @@ -535,12 +529,11 @@ double *dft_energy::total() { } // TODO: change here -dft_energy fields::add_dft_energy(const volume_list *where_, double freq_min, double freq_max, - int Nfreq) { +dft_energy fields::add_dft_energy(const volume_list *where_, double *freqs, int Nfreq) { if (!where_) // handle empty list of volumes // TODO: change here - return dft_energy(NULL, NULL, NULL, NULL, freq_min, freq_max, Nfreq, v); + return dft_energy(NULL, NULL, NULL, NULL, freqs, Nfreq, v); dft_chunk *E = 0, *D = 0, *H = 0, *B = 0; volume firstvol(where_->v); @@ -549,20 +542,20 @@ dft_energy fields::add_dft_energy(const volume_list *where_, double freq_min, do while (where) { LOOP_OVER_FIELD_DIRECTIONS(gv.dim, d) { // TODO: change here - E = add_dft(direction_component(Ex, d), where->v, freq_min, freq_max, Nfreq, true, 1.0, E); + E = add_dft(direction_component(Ex, d), where->v, freqs, Nfreq, true, 1.0, E); // TODO: change here - D = add_dft(direction_component(Dx, d), where->v, freq_min, freq_max, Nfreq, false, 1.0, D); + D = add_dft(direction_component(Dx, d), where->v, freqs, Nfreq, false, 1.0, D); // TODO: change here - H = add_dft(direction_component(Hx, d), where->v, freq_min, freq_max, Nfreq, true, 1.0, H); + H = add_dft(direction_component(Hx, d), where->v, freqs, Nfreq, true, 1.0, H); // TODO: change here - B = add_dft(direction_component(Bx, d), where->v, freq_min, freq_max, Nfreq, false, 1.0, B); + B = add_dft(direction_component(Bx, d), where->v, freqs, Nfreq, false, 1.0, B); } where = where->next; } delete where_save; // TODO: change here - return dft_energy(E, H, D, B, freq_min, freq_max, Nfreq, firstvol); + return dft_energy(E, H, D, B, freqs, Nfreq, firstvol); } void dft_energy::save_hdf5(h5file *file, const char *dprefix) { @@ -644,26 +637,24 @@ direction fields::normal_direction(const volume &where) const { } // TODO: change here -dft_flux fields::add_dft_flux(direction d, const volume &where, double freq_min, double freq_max, - int Nfreq, bool use_symmetry) { +dft_flux fields::add_dft_flux(direction d, const volume &where, double *freqs, int Nfreq, + bool use_symmetry) { if (d == NO_DIRECTION) d = normal_direction(where); volume_list vl(where, direction_component(Sx, d)); // TODO: change here - dft_flux flux = add_dft_flux(&vl, freq_min, freq_max, Nfreq, use_symmetry); + dft_flux flux = add_dft_flux(&vl, freqs, Nfreq, use_symmetry); flux.normal_direction = d; return flux; } // TODO: change here -dft_flux fields::add_mode_monitor(direction d, const volume &where, double freq_min, - double freq_max, int Nfreq) { +dft_flux fields::add_mode_monitor(direction d, const volume &where, double *freqs, int Nfreq) { // TODO: change here - return add_dft_flux(d, where, freq_min, freq_max, Nfreq, /*use_symmetry=*/false); + return add_dft_flux(d, where, freqs, Nfreq, /*use_symmetry=*/false); } // TODO: change here -dft_flux fields::add_dft_flux_box(const volume &where, double freq_min, double freq_max, - int Nfreq) { +dft_flux fields::add_dft_flux_box(const volume &where, double *freqs, int Nfreq) { volume_list *faces = 0; LOOP_OVER_DIRECTIONS(where.dim, d) { if (where.in_direction(d) > 0) { @@ -678,27 +669,24 @@ dft_flux fields::add_dft_flux_box(const volume &where, double freq_min, double f } // TODO: change here - dft_flux flux = add_dft_flux(faces, freq_min, freq_max, Nfreq); + dft_flux flux = add_dft_flux(faces, freqs, Nfreq); delete faces; return flux; } // TODO: change here -dft_flux fields::add_dft_flux_plane(const volume &where, double freq_min, double freq_max, - int Nfreq) { +dft_flux fields::add_dft_flux_plane(const volume &where, double *freqs, int Nfreq) { // TODO: change here - return add_dft_flux(NO_DIRECTION, where, freq_min, freq_max, Nfreq); + return add_dft_flux(NO_DIRECTION, where, freqs, Nfreq); } // TODO: change here -dft_fields::dft_fields(dft_chunk *chunks_, double freq_min_, double freq_max_, int Nfreq_, - const volume &where_) +dft_fields::dft_fields(dft_chunk *chunks_, double *freqs_, int Nfreq_, const volume &where_) : where(where_) { chunks = chunks_; // TODO: change here - freq_min = freq_min_; - dfreq = Nfreq_ <= 1 ? 0.0 : (freq_max_ - freq_min_) / (Nfreq_ - 1); Nfreq = Nfreq_; + freqs = freqs_; } void dft_fields::scale_dfts(cdouble scale) { chunks->scale_dft(scale); } @@ -713,7 +701,7 @@ void dft_fields::remove() { // TODO: change here dft_fields fields::add_dft_fields(component *components, int num_components, const volume where, - double freq_min, double freq_max, int Nfreq, bool use_centered_grid) { + double *freqs, int Nfreq, bool use_centered_grid) { bool include_dV_and_interp_weights = false; bool sqrt_dV_and_interp_weights = false; // default option from meep.hpp (expose to user?) std::complex extra_weight = 1.0; // default option from meep.hpp (expose to user?) @@ -721,12 +709,12 @@ dft_fields fields::add_dft_fields(component *components, int num_components, con dft_chunk *chunks = 0; for (int nc = 0; nc < num_components; nc++) // TODO: change here - chunks = add_dft(components[nc], where, freq_min, freq_max, Nfreq, - include_dV_and_interp_weights, stored_weight, chunks, - sqrt_dV_and_interp_weights,extra_weight,use_centered_grid); + chunks = add_dft(components[nc], where, freqs, Nfreq, include_dV_and_interp_weights, + stored_weight, chunks, sqrt_dV_and_interp_weights, extra_weight, + use_centered_grid); // TODO: change here - return dft_fields(chunks, freq_min, freq_max, Nfreq, where); + return dft_fields(chunks, freqs, Nfreq, where); } /***************************************************************/ From 286e22fb61960491a00e59d8a1116711d8303985 Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Mon, 2 Mar 2020 17:13:21 -0700 Subject: [PATCH 06/35] greater detail in TODOs; add overloaded versions of constructors, helper functions, etc. to header --- src/meep.hpp | 87 +++++++++++++++++++++++++++++++++++++--------------- 1 file changed, 62 insertions(+), 25 deletions(-) diff --git a/src/meep.hpp b/src/meep.hpp index 1fc205685..1b993167d 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1009,7 +1009,7 @@ class dft_chunk { void operator-=(const dft_chunk &chunk); // the frequencies to loop_in_chunks - // TODO: change here + // TODO: change dft_chunk properties int Nomega; double omegas[Nomega]; @@ -1082,9 +1082,12 @@ void load_dft_hdf5(dft_chunk *dft_chunks, const char *name, h5file *file, const // dft.cpp (normally created with fields::add_dft_flux) class dft_flux { public: - // TODO: change here + // TODO: add dft_flux constructor dft_flux(const component cE_, const component cH_, dft_chunk *E_, dft_chunk *H_, double *fs, int Nf, const volume &where_, direction normal_direction_, bool use_symmetry_); + dft_flux(const component cE_, const component cH_, dft_chunk *E_, dft_chunk *H_, double fmin, + double fmax, int Nf, const volume &where_, direction normal_direction_, + bool use_symmetry_); dft_flux(const dft_flux &f); double *flux(); @@ -1104,7 +1107,7 @@ class dft_flux { void remove(); - // TODO: change here + // TODO: change dft_flux properties int Nfreq; double freqs[Nfreq]; dft_chunk *E, *H; @@ -1117,9 +1120,11 @@ class dft_flux { // dft.cpp (normally created with fields::add_dft_energy) class dft_energy { public: - // TODO: change here + // TODO: add dft_energy constructor dft_energy(dft_chunk *E_, dft_chunk *H_, dft_chunk *D_, dft_chunk *B_, double *fs, int Nf, const volume &where_); + dft_energy(dft_chunk *E_, dft_chunk *H_, dft_chunk *D_, dft_chunk *B_, double fmin, double fmax, + int Nf, const volume &where_); dft_energy(const dft_energy &f); double *electric(); @@ -1143,7 +1148,7 @@ class dft_energy { void remove(); - // TODO: change here + // TODO: change dft_energy properties int Nfreq; double freqs[Nfreq]; dft_chunk *E, *H, *D, *B; @@ -1153,9 +1158,11 @@ class dft_energy { // stress.cpp (normally created with fields::add_dft_force) class dft_force { public: - // TODO: change here + // TODO: add dft_force constructor dft_force(dft_chunk *offdiag1_, dft_chunk *offdiag2_, dft_chunk *diag_, double *fs, int Nf, const volume &where_); + dft_force::dft_force(dft_chunk *offdiag1_, dft_chunk *offdiag2_, dft_chunk *diag_, double fmin, + double fmax, int Nf, const volume &where_) dft_force(const dft_force &f); double *force(); @@ -1172,7 +1179,7 @@ class dft_force { void remove(); - // TODO: change here + // TODO: change dft_force properties int Nfreq; double freqs[Nfreq]; dft_chunk *offdiag1, *offdiag2, *diag; @@ -1184,10 +1191,13 @@ class dft_near2far { public: /* fourier tranforms of tangential E and H field components in a medium with the given scalar eps and mu */ - // TODO: change here + // TODO: add dft_near2far constructor dft_near2far(dft_chunk *F, double *fs, int Nf, double eps, double mu, const volume &where_, const direction periodic_d_[2], const int periodic_n_[2], const double periodic_k_[2], const double period_[2]); + dft_near2far(dft_chunk *F, double fmin, double fmax, int Nf, double eps, double mu, + const volume &where_, const direction periodic_d_[2], const int periodic_n_[2], + const double periodic_k_[2], const double period_[2]); dft_near2far(const dft_near2far &f); /* return an array (Ex,Ey,Ez,Hx,Hy,Hz) x Nfreq of the far fields at x */ @@ -1222,7 +1232,7 @@ class dft_near2far { void remove(); - // TODO: change here + // TODO: change dft_near2far properties int Nfreq; double freqs[Nfreq]; dft_chunk *F; @@ -1241,8 +1251,9 @@ class dft_near2far { store the Fourier transform per point or per current. */ class dft_ldos { public: - // TODO: change here + // TODO: add dft_lods constructor dft_ldos(double *freqs, int Nfreq); + dft_ldos(double freq_min, double freq_max, int Nfreq); ~dft_ldos() { delete[] Fdft; delete[] Jdft; @@ -1258,7 +1269,7 @@ class dft_ldos { std::complex *Jdft; // Nomega array of J(t) DFT values double Jsum; // sum of |J| over all points public: - // TODO: change here + // TODO: change dft_ldos properties int Nomega; double omegas[Nomega]; }; @@ -1266,14 +1277,15 @@ class dft_ldos { // dft.cpp (normally created with fields::add_dft_fields) class dft_fields { public: - // TODO: change here + // TODO: add dft_fields constructor dft_fields(dft_chunk *chunks, double *freqs, int Nfreq, const volume &where); + dft_fields(dft_chunk *chunks, double freq_min, double freq_max, int Nfreq, const volume &where); void scale_dfts(std::complex scale); void remove(); - // TODO: change here + // TODO: change dft_fields properties int Nfreq; double freqs[Nfreq]; dft_chunk *chunks; @@ -1733,36 +1745,57 @@ class fields { double max_abs(derived_component c, const volume &where); // dft.cpp - // TODO: change here + // TODO: overload add_dft dft_chunk *add_dft(component c, const volume &where, double *freqs, int Nfreq, bool include_dV_and_interp_weights = true, std::complex stored_weight = 1.0, dft_chunk *chunk_next = 0, bool sqrt_dV_and_interp_weights = false, std::complex extra_weight = 1.0, bool use_centered_grid = true, int vc = 0); - // TODO: change here + dft_chunk *add_dft(component c, const volume &where, double freq_min, double freq_max, int Nfreq, + bool include_dV_and_interp_weights = true, + std::complex stored_weight = 1.0, dft_chunk *chunk_next = 0, + bool sqrt_dV_and_interp_weights = false, + std::complex extra_weight = 1.0, bool use_centered_grid = true, + int vc = 0); + // TODO: overload add_dft_pt dft_chunk *add_dft_pt(component c, const vec &where, double *freqs, int Nfreq); - // TODO: change here + dft_chunk *add_dft_pt(component c, const vec &where, double freq_min, double freq_max, int Nfreq); + // TODO: overload add_dft dft_chunk *add_dft(const volume_list *where, double *freqs, int Nfreq, bool include_dV = true); + dft_chunk *add_dft(const volume_list *where, double freq_min, double freq_max, int Nfreq, + bool include_dV = true); void update_dfts(); - // TODO: change here + // TODO: overload add_dft_flux dft_flux add_dft_flux(const volume_list *where, double *freqs, int Nfreq, bool use_symmetry = true); - // TODO: change here + dft_flux add_dft_flux(const volume_list *where, double freq_min, double freq_max, int Nfreq, + bool use_symmetry = true); + // TODO: overload add_dft_flux dft_flux add_dft_flux(direction d, const volume &where, double *freqs, int Nfreq, bool use_symmetry = true); - // TODO: change here + dft_flux add_dft_flux(direction d, const volume &where, double freq_min, double freq_max, + int Nfreq, bool use_symmetry = true); + // TODO: overload add_dft_flux_box dft_flux add_dft_flux_box(const volume &where, double *freqs, int Nfreq); - // TODO: change here + dft_flux add_dft_flux_box(const volume &where, double freq_min, double freq_max, int Nfreq); + // TODO: overload add_dft_flux_plane dft_flux add_dft_flux_plane(const volume &where, double *freqs, int Nfreq); + dft_flux add_dft_flux_plane(const volume &where, double freq_min, double freq_max, int Nfreq); // a "mode monitor" is just a dft_flux with symmetry reduction turned off. - // TODO: change here + // TODO: overload add_mode_monitor dft_flux add_mode_monitor(direction d, const volume &where, double *freqs, int Nfreq); + dft_flux add_mode_monitor(direction d, const volume &where, double freq_min, double freq_max, + int Nfreq); - // TODO: change here + + // TODO: overload add_dft_fields dft_fields add_dft_fields(component *components, int num_components, const volume where, double *freqs, int Nfreq, bool use_centered_grid=true); + dft_fields add_dft_fields(component *components, int num_components, const volume where, + double freq_min, double freq_max, int Nfreq, + bool use_centered_grid=true); /********************************************************/ /* process_dft_component is an intermediate-level */ @@ -1808,17 +1841,21 @@ class fields { void get_mode_mode_overlap(void *mode1_data, void *mode2_data, dft_flux flux, std::complex overlaps[2]); - // TODO: change here + // TODO: overload add_dft_energy dft_energy add_dft_energy(const volume_list *where, double *freqs, int Nfreq); + dft_energy add_dft_energy(const volume_list *where, double freq_min, double freq_max, int Nfreq); // stress.cpp - // TODO: change here + // TODO: overload add_dft_force dft_force add_dft_force(const volume_list *where, double *freqs, int Nfreq); + dft_force add_dft_force(const volume_list *where, double freq_min, double freq_max, int Nfreq); // near2far.cpp - // TODO: change here + // TODO: overload add_dft_near2far dft_near2far add_dft_near2far(const volume_list *where, double *freqs, int Nfreq, int Nperiods = 1); + dft_near2far add_dft_near2far(const volume_list *where, double freq_min, double freq_max, + int Nfreq, int Nperiods = 1); // monitor.cpp std::complex get_chi1inv(component, direction, const vec &loc, double omega = 0, bool parallel = true) const; From ea7ba082bd1df03e83cd228d2290ca6d51226399 Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Mon, 2 Mar 2020 17:18:44 -0700 Subject: [PATCH 07/35] greater detail in TODOs; add overloaded versions of constructors, helper functions, etc. to header --- src/meep.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/meep.hpp b/src/meep.hpp index 1b993167d..054750286 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1088,7 +1088,7 @@ class dft_flux { dft_flux(const component cE_, const component cH_, dft_chunk *E_, dft_chunk *H_, double fmin, double fmax, int Nf, const volume &where_, direction normal_direction_, bool use_symmetry_); - dft_flux(const dft_flux &f); + dft_flux(const dft_flux &f); double *flux(); From 92895cfed6d9b915e2d495a43ac861d9d64cc290 Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Mon, 2 Mar 2020 20:30:28 -0700 Subject: [PATCH 08/35] typo --- src/meep.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/meep.hpp b/src/meep.hpp index 054750286..d579227ae 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1161,8 +1161,8 @@ class dft_force { // TODO: add dft_force constructor dft_force(dft_chunk *offdiag1_, dft_chunk *offdiag2_, dft_chunk *diag_, double *fs, int Nf, const volume &where_); - dft_force::dft_force(dft_chunk *offdiag1_, dft_chunk *offdiag2_, dft_chunk *diag_, double fmin, - double fmax, int Nf, const volume &where_) + dft_force(dft_chunk *offdiag1_, dft_chunk *offdiag2_, dft_chunk *diag_, double fmin, double fmax, + int Nf, const volume &where_) dft_force(const dft_force &f); double *force(); From 4cfc13b2c02ac4fcc4d2cb89ffa9439291b49a0e Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Mon, 2 Mar 2020 20:31:19 -0700 Subject: [PATCH 09/35] details on TODOs, added some overloaded constructors --- src/dft.cpp | 24 ++++++++++++++++++++---- 1 file changed, 20 insertions(+), 4 deletions(-) diff --git a/src/dft.cpp b/src/dft.cpp index 2ca77041b..464de81d6 100644 --- a/src/dft.cpp +++ b/src/dft.cpp @@ -33,7 +33,7 @@ namespace meep { struct dft_chunk_data { // for passing to field::loop_in_chunks as void* component c; int vc; - // TODO: change here + // TODO: change dft_chunk_data properties int Nomega; double omegas[Nomega]; complex stored_weight, extra_weight; @@ -85,7 +85,7 @@ dft_chunk::dft_chunk(fields_chunk *fc_, ivec is_, ivec ie_, vec s0_, vec s1_, ve sn = sn_; vc = data->vc; - // TODO: change here + // TODO: change dft_chunk properties Nomega = data->Nomega; omegas = data->omegas; dft_phase = new complex[Nomega]; @@ -146,7 +146,7 @@ static void add_dft_chunkloop(fields_chunk *fc, int ichunk, component cgrid, ive shift_phase * S.phase_shift(c, sn), shift, S, sn, chunkloop_data); } -// TODO: change here +// TODO: overload fields::add_dft dft_chunk *fields::add_dft(component c, const volume &where, double *freqs, int Nfreq, bool include_dV_and_interp_weights, complex stored_weight, dft_chunk *chunk_next, bool sqrt_dV_and_interp_weights, @@ -165,7 +165,7 @@ dft_chunk *fields::add_dft(component c, const volume &where, double *freqs, int dft_chunk_data data; data.c = c; data.vc = vc; - // TODO: change here + // TODO: change how freqs becomes omegas if (Nfreq < 1) abort("Nfreq must be at least 1"); data.Nomega = Nfreq; double omegas[Nfreq]; @@ -186,6 +186,22 @@ dft_chunk *fields::add_dft(component c, const volume &where, double *freqs, int return data.dft_chunks; } +dft_chunk *fields::add_dft(component c, const volume &where, double freq_min, double freq_max, + int Nfreq, bool include_dV_and_interp_weights, + complex stored_weight, dft_chunk *chunk_next, + bool sqrt_dV_and_interp_weights, complex extra_weight, + bool use_centered_grid, int vc) { + if (Nfreq < 1) abort("Nfreq must be at least 1"); + double dfreq = (freq_max - freq_min) / (Nfreq - 1); + double freqs[Nfreq]; + int i; + for (i = 0; i < Nfreq; i++) { + freqs[i] = freq_min + dfreq * i; + } + return add_dft(c, where, freqs, Nfreq, include_dV_and_interp_weights, stored_weight, chunk_next, + sqrt_dV_and_interp_weights, extra_weight, use_centered_grid, vc); +} + // TODO: change here dft_chunk *fields::add_dft(const volume_list *where, double *freqs, int Nfreq, bool include_dV_and_interp_weights) { From 157638f6860b24f06bdf67b8ab069acf1227fd42 Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Mon, 2 Mar 2020 21:23:48 -0700 Subject: [PATCH 10/35] greater detail in TODOs; add overloaded versions of constructors, helper functions, etc. --- src/dft.cpp | 160 ++++++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 148 insertions(+), 12 deletions(-) diff --git a/src/dft.cpp b/src/dft.cpp index 464de81d6..e25193201 100644 --- a/src/dft.cpp +++ b/src/dft.cpp @@ -202,7 +202,7 @@ dft_chunk *fields::add_dft(component c, const volume &where, double freq_min, do sqrt_dV_and_interp_weights, extra_weight, use_centered_grid, vc); } -// TODO: change here +// TODO: overload add_dft dft_chunk *fields::add_dft(const volume_list *where, double *freqs, int Nfreq, bool include_dV_and_interp_weights) { dft_chunk *chunks = 0; @@ -217,12 +217,30 @@ dft_chunk *fields::add_dft(const volume_list *where, double *freqs, int Nfreq, return chunks; } -// TODO: change here +dft_chunk *fields::add_dft(const volume_list *where, double freq_min, double freq_max, int Nfreq, + bool include_dV_and_interp_weights) { + dft_chunk *chunks = 0; + while (where) { + if (is_derived(where->c)) abort("derived_component invalid for dft"); + cdouble stored_weight = where->weight; + chunks = add_dft(component(where->c), where->v, freq_min, freq_max, Nfreq, + include_dV_and_interp_weights, stored_weight, chunks); + where = where->next; + } + return chunks; +} + +// TODO: overload add_dft_pt dft_chunk *fields::add_dft_pt(component c, const vec &where, double *freqs, int Nfreq) { // TODO: change here return add_dft(c, where, freqs, Nfreq, false); } +dft_chunk *fields::add_dft_pt(component c, const vec &where, double freq_min, double freq_max, + int Nfreq) { + return add_dft(c, where, freq_min, freq_max, Nfreq, false); +} + void fields::update_dfts() { am_now_working_on(FourierTransforming); for (int i = 0; i < num_chunks; i++) @@ -361,7 +379,7 @@ void load_dft_hdf5(dft_chunk *dft_chunks, component c, h5file *file, const char load_dft_hdf5(dft_chunks, component_name(c), file, dprefix); } -// TODO: change here +// TODO: overload dft_flux constructor dft_flux::dft_flux(const component cE_, const component cH_, dft_chunk *E_, dft_chunk *H_, double *fs, int Nf, const volume &where_, direction normal_direction_, bool use_symmetry_) @@ -372,6 +390,21 @@ dft_flux::dft_flux(const component cE_, const component cH_, dft_chunk *E_, dft_ freqs = fs; } +dft_flux::dft_flux(const component cE_, const component cH_, dft_chunk *E_, dft_chunk *H_, + double fmin, double fmax, int Nf, const volume &where_, + direction normal_direction_, bool use_symmetry_) + : Nfreq(Nf), E(E_), H(H_), cE(cE_), cH(cH_), where(where_), normal_direction(normal_direction_), + use_symmetry(use_symmetry_) { + if (Nf < 1) abort("Nf must be at least 1"); + double dfreqs = (fmax - fmin) / (Nf - 1); + double fs[Nf]; + int i; + for (i = 0; i < Nf; i++) { + fs[i] = fmin + dfreqs * i; + } + freqs = fs; +} + dft_flux::dft_flux(const dft_flux &f) : where(f.where) { // TODO: change here Nfreq = f.Nfreq; @@ -428,7 +461,7 @@ void dft_flux::scale_dfts(complex scale) { if (H) H->scale_dft(scale); } -// TODO: change here +// TODO: overload fields::add_dft_flux dft_flux fields::add_dft_flux(const volume_list *where_, double *freqs, int Nfreq, bool use_symmetry) { if (!where_) // handle empty list of volumes @@ -483,7 +516,19 @@ dft_flux fields::add_dft_flux(const volume_list *where_, double *freqs, int Nfre return dft_flux(cE[0], cH[0], E, H, freqs, Nfreq, firstvol, flux_dir, use_symmetry); } -// TODO: change here +dft_flux fields::add_dft_flux(const volume_list *where_, double freq_min, double freq_max, + int Nfreq, bool use_symmetry) { + if (Nfreq < 1) abort("Nfreq must be at least 1"); + double dfreqs = (fmax - fmin) / (Nfreq - 1); + double freqs[Nfreq]; + int i; + for (i = 0; i < Nfreq; i++) { + freqs[i] = fmin + dfreqs * i; + } + return add_dft_flux(where_, freqs, Nfreq, use_symmetry); +} + +// TODO: overload dft_energy constructor dft_energy::dft_energy(dft_chunk *E_, dft_chunk *H_, dft_chunk *D_, dft_chunk *B_, double *fs, int Nf, const volume &where_) : Nfreq(Nf), E(E_), H(H_), D(D_), B(B_), where(where_) { @@ -492,6 +537,19 @@ dft_energy::dft_energy(dft_chunk *E_, dft_chunk *H_, dft_chunk *D_, dft_chunk *B freqs = fs; } +dft_energy::dft_energy(dft_chunk *E_, dft_chunk *H_, dft_chunk *D_, dft_chunk *B_, double fmin, + double fmax, int Nf, const volume &where_) + : Nfreq(Nf), E(E_), H(H_), D(D_), B(B_), where(where_) { + if (Nf < 1) abort("Nf must be at least 1"); + dfreq = (fmax - fmin) / (Nf - 1); + double fs[Nf]; + int i; + for (i = 0; i < Nf; i++) { + fs[i] = fmin + dfreqs * i; + } + freqs = fs; +} + // TODO: change here dft_energy::dft_energy(const dft_energy &f) : where(f.where) { // TODO: change here @@ -544,7 +602,7 @@ double *dft_energy::total() { return F; } -// TODO: change here +// TODO: overload fields::add_dft_energy dft_energy fields::add_dft_energy(const volume_list *where_, double *freqs, int Nfreq) { if (!where_) // handle empty list of volumes @@ -574,6 +632,18 @@ dft_energy fields::add_dft_energy(const volume_list *where_, double *freqs, int return dft_energy(E, H, D, B, freqs, Nfreq, firstvol); } +dft_energy fields::add_dft_energy(const volume_list *where_, double freq_min, double freq_max, + int Nfreq) { + if (Nfreq < 1) abort("Nfreq must be at least 1"); + dfreq = (fmax - fmin) / (Nfreq - 1); + double freqs[Nfreq]; + int i; + for (i = 0; i < Nfreq; i++) { + freqs[i] = fmin + dfreqs * i; + } + return add_dft_energy(where_, freqs, Nfreq); +} + void dft_energy::save_hdf5(h5file *file, const char *dprefix) { save_dft_hdf5(E, "E", file, dprefix); file->prevent_deadlock(); // hackery @@ -652,7 +722,7 @@ direction fields::normal_direction(const volume &where) const { return d; } -// TODO: change here +// TODO: overload fields::add_dft_flux dft_flux fields::add_dft_flux(direction d, const volume &where, double *freqs, int Nfreq, bool use_symmetry) { if (d == NO_DIRECTION) d = normal_direction(where); @@ -663,13 +733,37 @@ dft_flux fields::add_dft_flux(direction d, const volume &where, double *freqs, i return flux; } -// TODO: change here +dft_flux fields::add_dft_flux(direction d, const volume &where, double freq_min, double freq_max, + int Nfreq, bool use_symmetry) { + if (Nfreq < 1) abort("Nfreq must be at least 1"); + dfreq = (fmax - fmin) / (Nfreq - 1); + double freqs[Nfreq]; + int i; + for (i = 0; i < Nfreq; i++) { + freqs[i] = fmin + dfreqs * i; + } + return add_dft_flux(d, where, freqs, Nfreq, use_symmetry); +} + +// TODO: overload fields::add_mode_monitor dft_flux fields::add_mode_monitor(direction d, const volume &where, double *freqs, int Nfreq) { // TODO: change here return add_dft_flux(d, where, freqs, Nfreq, /*use_symmetry=*/false); } -// TODO: change here +dft_flux fields::add_mode_monitor(direction d, const volume &where, double freq_min, + double freq_max, int Nfreq) { + if (Nfreq < 1) abort("Nfreq must be at least 1"); + dfreq = (fmax - fmin) / (Nfreq - 1); + double freqs[Nfreq]; + int i; + for (i = 0; i < Nfreq; i++) { + freqs[i] = fmin + dfreqs * i; + } + return add_dft_flux(d, where, freqs, Nfreq, /*use_symmetry=*/false); +} + +// TODO: overload fields::add_dft_flux_box dft_flux fields::add_dft_flux_box(const volume &where, double *freqs, int Nfreq) { volume_list *faces = 0; LOOP_OVER_DIRECTIONS(where.dim, d) { @@ -690,13 +784,30 @@ dft_flux fields::add_dft_flux_box(const volume &where, double *freqs, int Nfreq) return flux; } -// TODO: change here +dft_flux fields::add_dft_flux_box(const volume &where, double freq_min, double freq_max, + int Nfreq) { + if (Nfreq < 1) abort("Nfreq must be at least 1"); + dfreq = (fmax - fmin) / (Nfreq - 1); + double freqs[Nfreq]; + int i; + for (i = 0; i < Nfreq; i++) { + freqs[i] = fmin + dfreqs * i; + } + return add_dft_flux_box(where, freqs, Nfreq); +} + +// TODO: overload fields::add_dft_flux_plane dft_flux fields::add_dft_flux_plane(const volume &where, double *freqs, int Nfreq) { // TODO: change here return add_dft_flux(NO_DIRECTION, where, freqs, Nfreq); } -// TODO: change here +dft_flux fields::add_dft_flux_plane(const volume &where, double freq_min, double freq_max, + int Nfreq) { + return add_dft_flux(NO_DIRECTION, where, freq_min, freq_max, Nfreq); +} + +// TODO: overload dft_fields constructor dft_fields::dft_fields(dft_chunk *chunks_, double *freqs_, int Nfreq_, const volume &where_) : where(where_) { chunks = chunks_; @@ -705,6 +816,20 @@ dft_fields::dft_fields(dft_chunk *chunks_, double *freqs_, int Nfreq_, const vol freqs = freqs_; } +dft_fields::dft_fields(dft_chunk *chunks_, double freq_min_, double freq_max_, int Nfreq_, + const volume &where_) + : where(where_) { + chunks = chunks_; + Nfreq = Nfreq_; + double dfreqs = (freq_max_ - freq_min_) / (Nfreq_ - 1); + double freqs_[Nfreq_]; + int i; + for (i = 0; i < Nfreq_; i++) { + freqs_[i] = freq_min_ + dfreqs * i; + } + freqs = freqs_; +} + void dft_fields::scale_dfts(cdouble scale) { chunks->scale_dft(scale); } void dft_fields::remove() { @@ -715,7 +840,7 @@ void dft_fields::remove() { } } -// TODO: change here +// TODO: overload fields::add_dft_fields dft_fields fields::add_dft_fields(component *components, int num_components, const volume where, double *freqs, int Nfreq, bool use_centered_grid) { bool include_dV_and_interp_weights = false; @@ -733,6 +858,17 @@ dft_fields fields::add_dft_fields(component *components, int num_components, con return dft_fields(chunks, freqs, Nfreq, where); } +dft_fields fields::add_dft_fields(component *components, int num_components, const volume where, + double freq_min, double freq_max, int Nfreq, bool use_centered_grid) { + dfreq = (fmax - fmin) / (Nfreq - 1); + double freqs[Nfreq]; + int i; + for (i = 0; i < Nfreq; i++) { + freqs[i] = fmin + dfreqs * i; + } + return add_dft_fields(components, num_components, volume where, freqs, Nfreq, use_centered_grid); +} + /***************************************************************/ /* chunk-level processing for fields::process_dft_component. */ /***************************************************************/ From c24fa0aa047f1de0fd623b0e266c536a722c4f84 Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Mon, 2 Mar 2020 21:38:41 -0700 Subject: [PATCH 11/35] greater detail in TODOs; add overloaded versions of constructors, helper functions, etc. to header --- src/dft.cpp | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/dft.cpp b/src/dft.cpp index e25193201..00c863cd5 100644 --- a/src/dft.cpp +++ b/src/dft.cpp @@ -804,7 +804,14 @@ dft_flux fields::add_dft_flux_plane(const volume &where, double *freqs, int Nfre dft_flux fields::add_dft_flux_plane(const volume &where, double freq_min, double freq_max, int Nfreq) { - return add_dft_flux(NO_DIRECTION, where, freq_min, freq_max, Nfreq); + if (Nfreq < 1) abort("Nfreq must be at least 1"); + dfreq = (fmax - fmin) / (Nfreq - 1); + double freqs[Nfreq]; + int i; + for (i = 0; i < Nfreq; i++) { + freqs[i] = fmin + dfreqs * i; + } + return add_dft_flux(NO_DIRECTION, where, freqs, Nfreq); } // TODO: overload dft_fields constructor From d89d799084332c2c889e4f3cf676065e9529e380 Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Tue, 3 Mar 2020 10:57:29 -0700 Subject: [PATCH 12/35] greater detail in TODOs; add overloaded versions of constructors, helper functions, etc. --- src/stress.cpp | 61 ++++++++++++++++++++++++++++++++++++-------------- 1 file changed, 44 insertions(+), 17 deletions(-) diff --git a/src/stress.cpp b/src/stress.cpp index c478309d6..3ff8a4454 100644 --- a/src/stress.cpp +++ b/src/stress.cpp @@ -23,14 +23,31 @@ using namespace std; namespace meep { -// TODO: change here +// TODO: overload dft_force constructor +dft_force::dft_force(dft_chunk *offdiag1_, dft_chunk *offdiag2_, dft_chunk *diag_, double *fs, + int Nf, const volume &where_) + : where(where_) { + if (Nf < 1) abort("Nf must be at least 1"); + freqs = fs; + Nfreq = Nf; + offdiag1 = offdiag1_; + offdiag2 = offdiag2_; + diag = diag_; +} + dft_force::dft_force(dft_chunk *offdiag1_, dft_chunk *offdiag2_, dft_chunk *diag_, double fmin, double fmax, int Nf, const volume &where_) : where(where_) { // TODO: change here - if (Nf <= 1) fmin = fmax = (fmin + fmax) * 0.5; - freq_min = fmin; + if (Nf < 1) abort("Nf must be at least 1"); + double dfreq = (fmax - fmin) / (Nf - 1); + double fs[Nf]; + int i; + for (i = 0; i < Nf; i++) { + fs[i] = fmin + dfreq * i; + } Nfreq = Nf; + freqs = fs; dfreq = Nf <= 1 ? 0.0 : (fmax - fmin) / (Nf - 1); offdiag1 = offdiag1_; offdiag2 = offdiag2_; @@ -40,9 +57,8 @@ dft_force::dft_force(dft_chunk *offdiag1_, dft_chunk *offdiag2_, dft_chunk *diag dft_force::dft_force(const dft_force &f) : where(f.where) { // TODO: change here - freq_min = f.freq_min; Nfreq = f.Nfreq; - dfreq = f.dfreq; + freqs = f.freqs; offdiag1 = f.offdiag1; offdiag2 = f.offdiag2; diag = f.diag; @@ -134,9 +150,8 @@ void dft_force::scale_dfts(complex scale) { /* note that the components where->c indicate the direction of the force to be computed, so they should be vector components (such as Ex, Ey, ... or Sx, ...) rather than pseudovectors (like Hx, ...). */ -// TODO: change here -dft_force fields::add_dft_force(const volume_list *where_, double freq_min, double freq_max, - int Nfreq) { +// TODO: overload fields::add_dft_force +dft_force fields::add_dft_force(const volume_list *where_, double *freqs, int Nfreq) { dft_chunk *offdiag1 = 0, *offdiag2 = 0, *diag = 0; volume_list *where = S.reduce(where_); @@ -152,34 +167,46 @@ dft_force fields::add_dft_force(const volume_list *where_, double freq_min, doub if (fd != nd) { // off-diagaonal stress-tensor terms // TODO: change here - offdiag1 = add_dft(direction_component(Ex, fd), where->v, freq_min, freq_max, Nfreq, true, + offdiag1 = add_dft(direction_component(Ex, fd), where->v, freqs, Nfreq, true, where->weight, offdiag1); // TODO: change here - offdiag2 = add_dft(direction_component(Ex, nd), where->v, freq_min, freq_max, Nfreq, false, + offdiag2 = add_dft(direction_component(Ex, nd), where->v, freqs, Nfreq, false, 1.0, offdiag2); // TODO: change here - offdiag1 = add_dft(direction_component(Hx, fd), where->v, freq_min, freq_max, Nfreq, true, + offdiag1 = add_dft(direction_component(Hx, fd), where->v, freqs, Nfreq, true, where->weight, offdiag1); // TODO: change here - offdiag2 = add_dft(direction_component(Hx, nd), where->v, freq_min, freq_max, Nfreq, false, + offdiag2 = add_dft(direction_component(Hx, nd), where->v, freqs, Nfreq, false, 1.0, offdiag2); } else // diagonal stress-tensor terms LOOP_OVER_FIELD_DIRECTIONS(gv.dim, d) { complex weight1 = where->weight * (d == fd ? +0.5 : -0.5); // TODO: change here - diag = add_dft(direction_component(Ex, d), where->v, freq_min, freq_max, Nfreq, true, 1.0, - diag, true, weight1, false); + diag = add_dft(direction_component(Ex, d), where->v, freqs, Nfreq, true, 1.0, diag, true, + weight1, false); // TODO: change here - diag = add_dft(direction_component(Hx, d), where->v, freq_min, freq_max, Nfreq, true, 1.0, - diag, true, weight1, false); + diag = add_dft(direction_component(Hx, d), where->v, freqs, Nfreq, true, 1.0, diag, true, + weight1, false); } everywhere = everywhere | where->v; } delete where_save; // TODO: change here - return dft_force(offdiag1, offdiag2, diag, freq_min, freq_max, Nfreq, everywhere); + return dft_force(offdiag1, offdiag2, diag, freqs, Nfreq, everywhere); +} + +dft_force fields::add_dft_force(const volume_list *where_, double freq_min, double freq_max, + int Nfreq) { + if (Nfreq < 1) abort("Nfreq must be at least 1"); + double dfreq = (freq_max - freq_min) / (Nfreq - 1); + double freqs[Nfreq]; + int i; + for (i = 0; i < Nfreq; i++) { + freqs[i] = freq_min + dfreq * i; + } + return add_dft_force(where_, freqs, Nfreq); } } // namespace meep From 1698f2d83c37aeba83639901050fb40b3e2ca990 Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Tue, 3 Mar 2020 11:45:08 -0700 Subject: [PATCH 13/35] greater detail in TODOs; add overloaded versions of constructors, helper functions, etc. --- src/dft.cpp | 36 +++++++++---------------- src/near2far.cpp | 68 ++++++++++++++++++++++++++++++++++-------------- src/stress.cpp | 6 ++--- 3 files changed, 63 insertions(+), 47 deletions(-) diff --git a/src/dft.cpp b/src/dft.cpp index 00c863cd5..ab3f903c6 100644 --- a/src/dft.cpp +++ b/src/dft.cpp @@ -169,8 +169,7 @@ dft_chunk *fields::add_dft(component c, const volume &where, double *freqs, int if (Nfreq < 1) abort("Nfreq must be at least 1"); data.Nomega = Nfreq; double omegas[Nfreq]; - int i; - for (i = 0; i < Nfreq; i++) { omegas[i] = freqs[i] * 2 * pi; } + for (int i = 0; i < Nfreq; i++) { omegas[i] = freqs[i] * 2 * pi; } data.omegas = omegas; data.stored_weight = stored_weight; data.extra_weight = extra_weight; @@ -194,8 +193,7 @@ dft_chunk *fields::add_dft(component c, const volume &where, double freq_min, do if (Nfreq < 1) abort("Nfreq must be at least 1"); double dfreq = (freq_max - freq_min) / (Nfreq - 1); double freqs[Nfreq]; - int i; - for (i = 0; i < Nfreq; i++) { + for (int i = 0; i < Nfreq; i++) { freqs[i] = freq_min + dfreq * i; } return add_dft(c, where, freqs, Nfreq, include_dV_and_interp_weights, stored_weight, chunk_next, @@ -398,8 +396,7 @@ dft_flux::dft_flux(const component cE_, const component cH_, dft_chunk *E_, dft_ if (Nf < 1) abort("Nf must be at least 1"); double dfreqs = (fmax - fmin) / (Nf - 1); double fs[Nf]; - int i; - for (i = 0; i < Nf; i++) { + for (int i = 0; i < Nf; i++) { fs[i] = fmin + dfreqs * i; } freqs = fs; @@ -521,8 +518,7 @@ dft_flux fields::add_dft_flux(const volume_list *where_, double freq_min, double if (Nfreq < 1) abort("Nfreq must be at least 1"); double dfreqs = (fmax - fmin) / (Nfreq - 1); double freqs[Nfreq]; - int i; - for (i = 0; i < Nfreq; i++) { + for (int i = 0; i < Nfreq; i++) { freqs[i] = fmin + dfreqs * i; } return add_dft_flux(where_, freqs, Nfreq, use_symmetry); @@ -543,8 +539,7 @@ dft_energy::dft_energy(dft_chunk *E_, dft_chunk *H_, dft_chunk *D_, dft_chunk *B if (Nf < 1) abort("Nf must be at least 1"); dfreq = (fmax - fmin) / (Nf - 1); double fs[Nf]; - int i; - for (i = 0; i < Nf; i++) { + for (int i = 0; i < Nf; i++) { fs[i] = fmin + dfreqs * i; } freqs = fs; @@ -637,8 +632,7 @@ dft_energy fields::add_dft_energy(const volume_list *where_, double freq_min, do if (Nfreq < 1) abort("Nfreq must be at least 1"); dfreq = (fmax - fmin) / (Nfreq - 1); double freqs[Nfreq]; - int i; - for (i = 0; i < Nfreq; i++) { + for (int i = 0; i < Nfreq; i++) { freqs[i] = fmin + dfreqs * i; } return add_dft_energy(where_, freqs, Nfreq); @@ -738,8 +732,7 @@ dft_flux fields::add_dft_flux(direction d, const volume &where, double freq_min, if (Nfreq < 1) abort("Nfreq must be at least 1"); dfreq = (fmax - fmin) / (Nfreq - 1); double freqs[Nfreq]; - int i; - for (i = 0; i < Nfreq; i++) { + for (int i = 0; i < Nfreq; i++) { freqs[i] = fmin + dfreqs * i; } return add_dft_flux(d, where, freqs, Nfreq, use_symmetry); @@ -756,8 +749,7 @@ dft_flux fields::add_mode_monitor(direction d, const volume &where, double freq_ if (Nfreq < 1) abort("Nfreq must be at least 1"); dfreq = (fmax - fmin) / (Nfreq - 1); double freqs[Nfreq]; - int i; - for (i = 0; i < Nfreq; i++) { + for (int i = 0; i < Nfreq; i++) { freqs[i] = fmin + dfreqs * i; } return add_dft_flux(d, where, freqs, Nfreq, /*use_symmetry=*/false); @@ -789,8 +781,7 @@ dft_flux fields::add_dft_flux_box(const volume &where, double freq_min, double f if (Nfreq < 1) abort("Nfreq must be at least 1"); dfreq = (fmax - fmin) / (Nfreq - 1); double freqs[Nfreq]; - int i; - for (i = 0; i < Nfreq; i++) { + for (int i = 0; i < Nfreq; i++) { freqs[i] = fmin + dfreqs * i; } return add_dft_flux_box(where, freqs, Nfreq); @@ -807,8 +798,7 @@ dft_flux fields::add_dft_flux_plane(const volume &where, double freq_min, double if (Nfreq < 1) abort("Nfreq must be at least 1"); dfreq = (fmax - fmin) / (Nfreq - 1); double freqs[Nfreq]; - int i; - for (i = 0; i < Nfreq; i++) { + for (int i = 0; i < Nfreq; i++) { freqs[i] = fmin + dfreqs * i; } return add_dft_flux(NO_DIRECTION, where, freqs, Nfreq); @@ -830,8 +820,7 @@ dft_fields::dft_fields(dft_chunk *chunks_, double freq_min_, double freq_max_, i Nfreq = Nfreq_; double dfreqs = (freq_max_ - freq_min_) / (Nfreq_ - 1); double freqs_[Nfreq_]; - int i; - for (i = 0; i < Nfreq_; i++) { + for (int i = 0; i < Nfreq_; i++) { freqs_[i] = freq_min_ + dfreqs * i; } freqs = freqs_; @@ -869,8 +858,7 @@ dft_fields fields::add_dft_fields(component *components, int num_components, con double freq_min, double freq_max, int Nfreq, bool use_centered_grid) { dfreq = (fmax - fmin) / (Nfreq - 1); double freqs[Nfreq]; - int i; - for (i = 0; i < Nfreq; i++) { + for (int i = 0; i < Nfreq; i++) { freqs[i] = fmin + dfreqs * i; } return add_dft_fields(components, num_components, volume where, freqs, Nfreq, use_centered_grid); diff --git a/src/near2far.cpp b/src/near2far.cpp index 540ed388d..dfa49d1f9 100644 --- a/src/near2far.cpp +++ b/src/near2far.cpp @@ -29,16 +29,36 @@ using namespace std; namespace meep { -// TODO: change here +// TODO: overload dft_near2far constructor +dft_near2far::dft_near2far(dft_chunk *F_, double *fs, int Nf, double eps_, double mu_, + const volume &where_, const direction periodic_d_[2], + const int periodic_n_[2], const double periodic_k_[2], + const double period_[2]) + : Nfreq(Nf), F(F_), eps(eps_), mu(mu_), where(where_) { + // TODO: initialize freqs with fs if Nf >= 1 + if (Nf < 1) abort("Nf must be at least 1"); + freqs = fs; + for (int i = 0; i < 2; ++i) { + periodic_d[i] = periodic_d_[i]; + periodic_n[i] = periodic_n_[i]; + periodic_k[i] = periodic_k_[i]; + period[i] = period_[i]; + } +} + dft_near2far::dft_near2far(dft_chunk *F_, double fmin, double fmax, int Nf, double eps_, double mu_, const volume &where_, const direction periodic_d_[2], const int periodic_n_[2], const double periodic_k_[2], const double period_[2]) : Nfreq(Nf), F(F_), eps(eps_), mu(mu_), where(where_) { - // TODO: change here - if (Nf <= 1) fmin = fmax = (fmin + fmax) * 0.5; - freq_min = fmin; - dfreq = Nf <= 1 ? 0.0 : (fmax - fmin) / (Nf - 1); + // TODO: create fs array with which to initalize freqs + if (Nf < 1) abort("Nf must be at least 1"); + double dfreq = (fmax - fmin) / (Nf - 1); + double fs[Nf]; + for (int i = 0; i < Nf; i++) { + fs[i] = fmin + dfreq * i; + } + freqs = fs; for (int i = 0; i < 2; ++i) { periodic_d[i] = periodic_d_[i]; periodic_n[i] = periodic_n_[i]; @@ -47,10 +67,9 @@ dft_near2far::dft_near2far(dft_chunk *F_, double fmin, double fmax, int Nf, doub } } -// TODO: change here +// TODO: change here for updated properties of dft_near2far dft_near2far::dft_near2far(const dft_near2far &f) - : freq_min(f.freq_min), dfreq(f.dfreq), Nfreq(f.Nfreq), F(f.F), eps(f.eps), mu(f.mu), - where(f.where) { + : freqs(f.freqs), Nfreq(f.Nfreq), F(f.F), eps(f.eps), mu(f.mu), where(f.where) { for (int i = 0; i < 2; ++i) { periodic_d[i] = f.periodic_d[i]; periodic_n[i] = f.periodic_n[i]; @@ -332,8 +351,8 @@ void dft_near2far::farfield_lowlevel(std::complex *EH, const vec &x) { #endif for (int i = 0; i < Nfreq; ++i) { std::complex EH6[6]; - // TODO: change here??? - double freq = freq_min + i * dfreq; + // TODO: replace freq_min + dfreq * i with freqs[i] + double freq = freqs[i]; size_t idx_dft = 0; LOOP_OVER_IVECS(f->fc->gv, f->is, f->ie, idx) { IVEC_LOOP_LOC(f->fc->gv, x0); @@ -527,9 +546,9 @@ double *dft_near2far::flux(direction df, const volume &where, double resolution) static double approxeq(double a, double b) { return fabs(a - b) < 0.5e-11 * (fabs(a) + fabs(b)); } -// TODO: change here -dft_near2far fields::add_dft_near2far(const volume_list *where, double freq_min, double freq_max, - int Nfreq, int Nperiods) { +// TODO: overload fields::add_dft_near2far +dft_near2far fields::add_dft_near2far(const volume_list *where, double *freqs, int Nfreq, + int Nperiods) { dft_chunk *F = 0; /* E and H chunks*/ double eps = 0, mu = 0; volume everywhere = where->v; @@ -603,16 +622,27 @@ dft_near2far fields::add_dft_near2far(const volume_list *where, double freq_min, double s = j == 0 ? 1 : -1; /* sign of n x c */ if (is_electric(c)) s = -s; - // TODO: change here - F = add_dft(c, w->v, freq_min, freq_max, Nfreq, true, s * w->weight, F, false, 1.0, false, - c0); + // TODO: replace freq_min and freq_max with freqs + F = add_dft(c, w->v, freqs, Nfreq, true, s * w->weight, F, false, 1.0, false, c0); } } } - // TODO: change here - return dft_near2far(F, freq_min, freq_max, Nfreq, eps, mu, everywhere, periodic_d, periodic_n, - periodic_k, period); + // TODO: replace freq_min and freq_max with freqs + return dft_near2far(F, freqs, Nfreq, eps, mu, everywhere, periodic_d, periodic_n, periodic_k, + period); + +} + +dft_near2far fields::add_dft_near2far(const volume_list *where, double freq_min, double freq_max, + int Nfreq, int Nperiods) { + if (Nfreq < 1) abort("Nfreq must be at least 1"); + double dfreq = (freq_max - freq_min) / (Nfreq - 1); + double freqs[Nfreq]; + for (int i = 0; i < Nfreq; i++) { + freqs[i] = freq_min + dfreq * i; + } + return add_dft_near2far(where, freqs, Nfreq, Nperiods); } } // namespace meep diff --git a/src/stress.cpp b/src/stress.cpp index 3ff8a4454..9ea04ed66 100644 --- a/src/stress.cpp +++ b/src/stress.cpp @@ -42,8 +42,7 @@ dft_force::dft_force(dft_chunk *offdiag1_, dft_chunk *offdiag2_, dft_chunk *diag if (Nf < 1) abort("Nf must be at least 1"); double dfreq = (fmax - fmin) / (Nf - 1); double fs[Nf]; - int i; - for (i = 0; i < Nf; i++) { + for (int i = 0; i < Nf; i++) { fs[i] = fmin + dfreq * i; } Nfreq = Nf; @@ -202,8 +201,7 @@ dft_force fields::add_dft_force(const volume_list *where_, double freq_min, doub if (Nfreq < 1) abort("Nfreq must be at least 1"); double dfreq = (freq_max - freq_min) / (Nfreq - 1); double freqs[Nfreq]; - int i; - for (i = 0; i < Nfreq; i++) { + for (int i = 0; i < Nfreq; i++) { freqs[i] = freq_min + dfreq * i; } return add_dft_force(where_, freqs, Nfreq); From b4f09268763b37587a0f829aa57afef5dd861216 Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Tue, 3 Mar 2020 15:09:57 -0700 Subject: [PATCH 14/35] Fixed how freqs is constructed and added freqs deletion to the remove functions --- src/dft.cpp | 10 +++++++++- src/meep.hpp | 17 +++++++++-------- src/near2far.cpp | 2 ++ src/stress.cpp | 2 ++ 4 files changed, 22 insertions(+), 9 deletions(-) diff --git a/src/dft.cpp b/src/dft.cpp index ab3f903c6..4bd75db62 100644 --- a/src/dft.cpp +++ b/src/dft.cpp @@ -35,7 +35,7 @@ struct dft_chunk_data { // for passing to field::loop_in_chunks as void* int vc; // TODO: change dft_chunk_data properties int Nomega; - double omegas[Nomega]; + double *omegas = new double[Nomega]; complex stored_weight, extra_weight; double dt_factor; bool include_dV_and_interp_weights; @@ -106,6 +106,8 @@ dft_chunk::dft_chunk(fields_chunk *fc_, ivec is_, ivec ie_, vec s0_, vec s1_, ve dft_chunk::~dft_chunk() { delete[] dft; delete[] dft_phase; + // TODO: add omegas to destructor? + delete[] omegas; // delete from fields_chunk list dft_chunk *cur = fc->dft_chunks; @@ -119,6 +121,8 @@ dft_chunk::~dft_chunk() { } void dft_flux::remove() { + // TODO: add freqs to remove? + delete[] freqs; while (E) { dft_chunk *nxt = E->next_in_dft; delete E; @@ -678,6 +682,8 @@ void dft_energy::scale_dfts(complex scale) { } void dft_energy::remove() { + // TODO: add freqs to remove? + delete[] freqs; while (E) { dft_chunk *nxt = E->next_in_dft; delete E; @@ -829,6 +835,8 @@ dft_fields::dft_fields(dft_chunk *chunks_, double freq_min_, double freq_max_, i void dft_fields::scale_dfts(cdouble scale) { chunks->scale_dft(scale); } void dft_fields::remove() { + // TODO: add freqs to remove? + delete[] freqs; while (chunks) { dft_chunk *nxt = chunks->next_in_dft; delete chunks; diff --git a/src/meep.hpp b/src/meep.hpp index d579227ae..ce3e25c02 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1011,7 +1011,7 @@ class dft_chunk { // the frequencies to loop_in_chunks // TODO: change dft_chunk properties int Nomega; - double omegas[Nomega]; + double *omegas = new double[Nomega]; component c; // component to DFT (possibly transformed by symmetry) @@ -1109,7 +1109,7 @@ class dft_flux { // TODO: change dft_flux properties int Nfreq; - double freqs[Nfreq]; + double *freqs = new double[Nfreq]; dft_chunk *E, *H; component cE, cH; volume where; @@ -1150,7 +1150,7 @@ class dft_energy { // TODO: change dft_energy properties int Nfreq; - double freqs[Nfreq]; + double *freqs = new double[Nfreq]; dft_chunk *E, *H, *D, *B; volume where; }; @@ -1162,7 +1162,7 @@ class dft_force { dft_force(dft_chunk *offdiag1_, dft_chunk *offdiag2_, dft_chunk *diag_, double *fs, int Nf, const volume &where_); dft_force(dft_chunk *offdiag1_, dft_chunk *offdiag2_, dft_chunk *diag_, double fmin, double fmax, - int Nf, const volume &where_) + int Nf, const volume &where_); dft_force(const dft_force &f); double *force(); @@ -1181,7 +1181,7 @@ class dft_force { // TODO: change dft_force properties int Nfreq; - double freqs[Nfreq]; + double *freqs = new double[Nfreq]; dft_chunk *offdiag1, *offdiag2, *diag; volume where; }; @@ -1234,7 +1234,7 @@ class dft_near2far { // TODO: change dft_near2far properties int Nfreq; - double freqs[Nfreq]; + double *freqs = new double[Nfreq]; dft_chunk *F; double eps, mu; volume where; @@ -1257,6 +1257,7 @@ class dft_ldos { ~dft_ldos() { delete[] Fdft; delete[] Jdft; + delete[] omegas; } void update(fields &f); // to be called after each timestep @@ -1271,7 +1272,7 @@ class dft_ldos { public: // TODO: change dft_ldos properties int Nomega; - double omegas[Nomega]; + double *omegas = new double[Nomega]; }; // dft.cpp (normally created with fields::add_dft_fields) @@ -1287,7 +1288,7 @@ class dft_fields { // TODO: change dft_fields properties int Nfreq; - double freqs[Nfreq]; + double *freqs = new double[Nfreq]; dft_chunk *chunks; volume where; }; diff --git a/src/near2far.cpp b/src/near2far.cpp index dfa49d1f9..a0e2583be 100644 --- a/src/near2far.cpp +++ b/src/near2far.cpp @@ -79,6 +79,8 @@ dft_near2far::dft_near2far(const dft_near2far &f) } void dft_near2far::remove() { + // TODO: add freqs to remove? + delete[] freqs; while (F) { dft_chunk *nxt = F->next_in_dft; delete F; diff --git a/src/stress.cpp b/src/stress.cpp index 9ea04ed66..73ca3aada 100644 --- a/src/stress.cpp +++ b/src/stress.cpp @@ -65,6 +65,8 @@ dft_force::dft_force(const dft_force &f) : where(f.where) { } void dft_force::remove() { + // TODO: add freqs to remove? + delete[] freqs; while (offdiag1) { dft_chunk *nxt = offdiag1->next_in_dft; delete offdiag1; From a429fd816dfd5b5d5b5b6542c2acf1beb752140d Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Tue, 3 Mar 2020 17:37:58 -0700 Subject: [PATCH 15/35] changed freqs definitions in classes to be just pointers. compiles until swig comes into play --- src/dft.cpp | 39 ++++++++++++++++++++------------------- src/dft_ldos.cpp | 33 ++++++++++++++++++++++++--------- src/meep.hpp | 15 ++++++++------- src/mpb.cpp | 5 ++--- src/stress.cpp | 2 +- 5 files changed, 55 insertions(+), 39 deletions(-) diff --git a/src/dft.cpp b/src/dft.cpp index 4bd75db62..6986858b3 100644 --- a/src/dft.cpp +++ b/src/dft.cpp @@ -87,7 +87,7 @@ dft_chunk::dft_chunk(fields_chunk *fc_, ivec is_, ivec ie_, vec s0_, vec s1_, ve // TODO: change dft_chunk properties Nomega = data->Nomega; - omegas = data->omegas; + omegas = data->omegas = new double [Nomega]; dft_phase = new complex[Nomega]; N = 1; @@ -520,10 +520,10 @@ dft_flux fields::add_dft_flux(const volume_list *where_, double *freqs, int Nfre dft_flux fields::add_dft_flux(const volume_list *where_, double freq_min, double freq_max, int Nfreq, bool use_symmetry) { if (Nfreq < 1) abort("Nfreq must be at least 1"); - double dfreqs = (fmax - fmin) / (Nfreq - 1); + double dfreq = (freq_max - freq_min) / (Nfreq - 1); double freqs[Nfreq]; for (int i = 0; i < Nfreq; i++) { - freqs[i] = fmin + dfreqs * i; + freqs[i] = freq_min + dfreq * i; } return add_dft_flux(where_, freqs, Nfreq, use_symmetry); } @@ -541,10 +541,10 @@ dft_energy::dft_energy(dft_chunk *E_, dft_chunk *H_, dft_chunk *D_, dft_chunk *B double fmax, int Nf, const volume &where_) : Nfreq(Nf), E(E_), H(H_), D(D_), B(B_), where(where_) { if (Nf < 1) abort("Nf must be at least 1"); - dfreq = (fmax - fmin) / (Nf - 1); + double dfreq = (fmax - fmin) / (Nf - 1); double fs[Nf]; for (int i = 0; i < Nf; i++) { - fs[i] = fmin + dfreqs * i; + fs[i] = fmin + dfreq * i; } freqs = fs; } @@ -634,10 +634,10 @@ dft_energy fields::add_dft_energy(const volume_list *where_, double *freqs, int dft_energy fields::add_dft_energy(const volume_list *where_, double freq_min, double freq_max, int Nfreq) { if (Nfreq < 1) abort("Nfreq must be at least 1"); - dfreq = (fmax - fmin) / (Nfreq - 1); + double dfreq = (freq_max - freq_min) / (Nfreq - 1); double freqs[Nfreq]; for (int i = 0; i < Nfreq; i++) { - freqs[i] = fmin + dfreqs * i; + freqs[i] = freq_min + dfreq * i; } return add_dft_energy(where_, freqs, Nfreq); } @@ -736,10 +736,10 @@ dft_flux fields::add_dft_flux(direction d, const volume &where, double *freqs, i dft_flux fields::add_dft_flux(direction d, const volume &where, double freq_min, double freq_max, int Nfreq, bool use_symmetry) { if (Nfreq < 1) abort("Nfreq must be at least 1"); - dfreq = (fmax - fmin) / (Nfreq - 1); + double dfreq = (freq_max - freq_min) / (Nfreq - 1); double freqs[Nfreq]; for (int i = 0; i < Nfreq; i++) { - freqs[i] = fmin + dfreqs * i; + freqs[i] = freq_min + dfreq * i; } return add_dft_flux(d, where, freqs, Nfreq, use_symmetry); } @@ -753,10 +753,10 @@ dft_flux fields::add_mode_monitor(direction d, const volume &where, double *freq dft_flux fields::add_mode_monitor(direction d, const volume &where, double freq_min, double freq_max, int Nfreq) { if (Nfreq < 1) abort("Nfreq must be at least 1"); - dfreq = (fmax - fmin) / (Nfreq - 1); + double dfreq = (freq_max - freq_min) / (Nfreq - 1); double freqs[Nfreq]; for (int i = 0; i < Nfreq; i++) { - freqs[i] = fmin + dfreqs * i; + freqs[i] = freq_min + dfreq * i; } return add_dft_flux(d, where, freqs, Nfreq, /*use_symmetry=*/false); } @@ -785,10 +785,10 @@ dft_flux fields::add_dft_flux_box(const volume &where, double *freqs, int Nfreq) dft_flux fields::add_dft_flux_box(const volume &where, double freq_min, double freq_max, int Nfreq) { if (Nfreq < 1) abort("Nfreq must be at least 1"); - dfreq = (fmax - fmin) / (Nfreq - 1); + double dfreq = (freq_max - freq_min) / (Nfreq - 1); double freqs[Nfreq]; for (int i = 0; i < Nfreq; i++) { - freqs[i] = fmin + dfreqs * i; + freqs[i] = freq_min + dfreq * i; } return add_dft_flux_box(where, freqs, Nfreq); } @@ -802,10 +802,10 @@ dft_flux fields::add_dft_flux_plane(const volume &where, double *freqs, int Nfre dft_flux fields::add_dft_flux_plane(const volume &where, double freq_min, double freq_max, int Nfreq) { if (Nfreq < 1) abort("Nfreq must be at least 1"); - dfreq = (fmax - fmin) / (Nfreq - 1); + double dfreq = (freq_max - freq_min) / (Nfreq - 1); double freqs[Nfreq]; for (int i = 0; i < Nfreq; i++) { - freqs[i] = fmin + dfreqs * i; + freqs[i] = freq_min + dfreq * i; } return add_dft_flux(NO_DIRECTION, where, freqs, Nfreq); } @@ -863,13 +863,14 @@ dft_fields fields::add_dft_fields(component *components, int num_components, con } dft_fields fields::add_dft_fields(component *components, int num_components, const volume where, - double freq_min, double freq_max, int Nfreq, bool use_centered_grid) { - dfreq = (fmax - fmin) / (Nfreq - 1); + double freq_min, double freq_max, int Nfreq, + bool use_centered_grid) { + double dfreq = (freq_max - freq_min) / (Nfreq - 1); double freqs[Nfreq]; for (int i = 0; i < Nfreq; i++) { - freqs[i] = fmin + dfreqs * i; + freqs[i] = freq_min + dfreq * i; } - return add_dft_fields(components, num_components, volume where, freqs, Nfreq, use_centered_grid); + return add_dft_fields(components, num_components, where, freqs, Nfreq, use_centered_grid); } /***************************************************************/ diff --git a/src/dft_ldos.cpp b/src/dft_ldos.cpp index add49c9b8..ec780f231 100644 --- a/src/dft_ldos.cpp +++ b/src/dft_ldos.cpp @@ -22,16 +22,31 @@ using namespace std; namespace meep { -dft_ldos::dft_ldos(double freq_min, double freq_max, int Nfreq) { - if (Nfreq <= 1) { - omega_min = (freq_min + freq_max) * pi; - domega = 0; - Nomega = 1; +// TODO: overload dft_ldos constructor +dft_ldos::dft_ldos(double *freqs, int Nfreq) { + if (Nfreq < 1) { abort("Nfreq must be at least 1"); } + else { + Nomega = Nfreq; + double om[Nomega]; + for (int i = 0; i < Nfreq; i++) { om[i] = freqs[i] * 2 * pi; } + omegas = om; } + Fdft = new complex[Nomega]; + Jdft = new complex[Nomega]; + for (int i = 0; i < Nomega; ++i) + Fdft[i] = Jdft[i] = 0.0; + Jsum = 1.0; +} + +// TODO: overload dft_ldos constructor +dft_ldos::dft_ldos(double freq_min, double freq_max, int Nfreq) { + if (Nfreq < 1) { abort("Nfreq must be at least 1"); } else { - omega_min = freq_min * 2 * pi; - domega = (freq_max - freq_min) * 2 * pi / (Nfreq - 1); + double dfreq = (freq_max - freq_min) / (Nfreq - 1); Nomega = Nfreq; + double om[Nomega]; + for (int i = 0; i < Nfreq; i++) { om[i] = (freq_min + dfreq * i) * 2 * pi; } + omegas = om; } Fdft = new complex[Nomega]; Jdft = new complex[Nomega]; @@ -130,8 +145,8 @@ void dft_ldos::update(fields &f) { } } for (int i = 0; i < Nomega; ++i) { - complex Ephase = polar(1.0, (omega_min + i * domega) * f.time()) * scale; - complex Hphase = polar(1.0, (omega_min + i * domega) * (f.time() - f.dt / 2)) * scale; + complex Ephase = polar(1.0, (omegas[i]) * f.time()) * scale; + complex Hphase = polar(1.0, (omegas[i]) * (f.time() - f.dt / 2)) * scale; Fdft[i] += Ephase * EJ + Hphase * HJ; // NOTE: take only 1st time dependence: assumes all sources have same J(t) diff --git a/src/meep.hpp b/src/meep.hpp index ce3e25c02..7de9b112c 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1011,7 +1011,7 @@ class dft_chunk { // the frequencies to loop_in_chunks // TODO: change dft_chunk properties int Nomega; - double *omegas = new double[Nomega]; + double *omegas; component c; // component to DFT (possibly transformed by symmetry) @@ -1109,7 +1109,7 @@ class dft_flux { // TODO: change dft_flux properties int Nfreq; - double *freqs = new double[Nfreq]; + double *freqs; dft_chunk *E, *H; component cE, cH; volume where; @@ -1150,7 +1150,7 @@ class dft_energy { // TODO: change dft_energy properties int Nfreq; - double *freqs = new double[Nfreq]; + double *freqs; dft_chunk *E, *H, *D, *B; volume where; }; @@ -1181,7 +1181,7 @@ class dft_force { // TODO: change dft_force properties int Nfreq; - double *freqs = new double[Nfreq]; + double *freqs; dft_chunk *offdiag1, *offdiag2, *diag; volume where; }; @@ -1234,7 +1234,7 @@ class dft_near2far { // TODO: change dft_near2far properties int Nfreq; - double *freqs = new double[Nfreq]; + double *freqs; dft_chunk *F; double eps, mu; volume where; @@ -1243,6 +1243,7 @@ class dft_near2far { double periodic_k[2], period[2]; }; +// dft_ldos.cpp /* Class to compute local-density-of-states spectra: the power spectrum P(omega) of the work done by the sources. Specialized to handle only the case where all sources have the same time dependence, which greatly @@ -1272,7 +1273,7 @@ class dft_ldos { public: // TODO: change dft_ldos properties int Nomega; - double *omegas = new double[Nomega]; + double *omegas; }; // dft.cpp (normally created with fields::add_dft_fields) @@ -1288,7 +1289,7 @@ class dft_fields { // TODO: change dft_fields properties int Nfreq; - double *freqs = new double[Nfreq]; + double *freqs; dft_chunk *chunks; volume where; }; diff --git a/src/mpb.cpp b/src/mpb.cpp index ea1637801..234298d6b 100644 --- a/src/mpb.cpp +++ b/src/mpb.cpp @@ -758,9 +758,8 @@ void fields::get_eigenmode_coefficients(dft_flux flux, const volume &eig_vol, in double *vgrp, kpoint_func user_kpoint_func, void *user_kpoint_data, vec *kpoints, vec *kdom_list, double *cscale, direction d) { - double freq_min = flux.freq_min; - double dfreq = flux.dfreq; int num_freqs = flux.Nfreq; + double *freqs = flux.freqs; bool match_frequency = true; if (flux.use_symmetry && S.multiplicity() > 1 && parity == 0) @@ -779,7 +778,7 @@ void fields::get_eigenmode_coefficients(dft_flux flux, const volume &eig_vol, in /*- call mpb to compute the eigenmode --------------------------*/ /*--------------------------------------------------------------*/ int band_num = bands[nb]; - double freq = freq_min + nf * dfreq; + double freq = freqs[nf]; double kdom[3]; if (user_kpoint_func) kpoint = user_kpoint_func(freq, band_num, user_kpoint_data); am_now_working_on(MPBTime); diff --git a/src/stress.cpp b/src/stress.cpp index 73ca3aada..364232f3f 100644 --- a/src/stress.cpp +++ b/src/stress.cpp @@ -40,12 +40,12 @@ dft_force::dft_force(dft_chunk *offdiag1_, dft_chunk *offdiag2_, dft_chunk *diag : where(where_) { // TODO: change here if (Nf < 1) abort("Nf must be at least 1"); + Nfreq = Nf; double dfreq = (fmax - fmin) / (Nf - 1); double fs[Nf]; for (int i = 0; i < Nf; i++) { fs[i] = fmin + dfreq * i; } - Nfreq = Nf; freqs = fs; dfreq = Nf <= 1 ? 0.0 : (fmax - fmin) / (Nf - 1); offdiag1 = offdiag1_; From 67c32050443e52b4e2e2cccb118e1405b6579f88 Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Sun, 8 Mar 2020 23:10:09 -0600 Subject: [PATCH 16/35] * Change properties of DftObj class. * Add TODOs for functions that need to be overloaded --- python/simulation.py | 40 ++++++++++++++++++++++++++++------------ 1 file changed, 28 insertions(+), 12 deletions(-) diff --git a/python/simulation.py b/python/simulation.py index a58e42430..56a2f0186 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -316,12 +316,8 @@ def remove(self): return self.swigobj_attr('remove') @property - def freq_min(self): - return self.swigobj_attr('freq_min') - - @property - def dfreq(self): - return self.swigobj_attr('dfreq') + def freqs(self): + return self.swigobj_attr('freqs') @property def Nfreq(self): @@ -721,8 +717,8 @@ def _check_material_frequencies(self): dft_freqs = [] for dftf in self.dft_objects: - dft_freqs.append(dftf.freq_min) - dft_freqs.append(dftf.freq_min + dftf.Nfreq * dftf.dfreq) + dft_freqs.append(dftf.freqs[0]) + dft_freqs.append(dftf.freqs[-1]) warn_src = ('Note: your sources include frequencies outside the range of validity of the ' + 'material models. This is fine as long as you eventually only look at outputs ' + @@ -1630,22 +1626,32 @@ def _evaluate_dft_objects(self): if dft.swigobj is None: dft.swigobj = dft.func(*dft.args) - def add_dft_fields(self, components, freq_min, freq_max, nfreq, where=None, center=None, size=None, yee_grid=False): + # TODO: overload add_dft_fields, _add_dft_fields + def add_dft_fields(self, components, freq_min=None, freq_max=None, freqs=None, nfreq, where=None, center=None, size=None, yee_grid=False): + if freqs is None: + if freq_min is None or freq_max is None: + raise ValueError("Either freqs or freq_min and freq_max is required for add_dft_fields") + else: + freqs = np.linspace(start=freq_min, end=freq_max, num=nfreq) + else: + freqs.sort() + freq_min = freqs[0] + freq_max = freqs[-1] center_v3 = Vector3(*center) if center is not None else None size_v3 = Vector3(*size) if size is not None else None use_centered_grid = not yee_grid - dftf = DftFields(self._add_dft_fields, [components, where, center_v3, size_v3, freq_min, freq_max, nfreq, use_centered_grid]) + dftf = DftFields(self._add_dft_fields, [components, where, center_v3, size_v3, freqs, nfreq, use_centered_grid]) self.dft_objects.append(dftf) return dftf - def _add_dft_fields(self, components, where, center, size, freq_min, freq_max, nfreq, use_centered_grid): + def _add_dft_fields(self, components, where, center, size, freq_min, freq_max, freqs, nfreq, use_centered_grid): if self.fields is None: self.init_sim() try: where = self._volume_from_kwargs(where, center, size) except ValueError: where = self.fields.total_volume() - return self.fields.add_dft_fields(components, where, freq_min, freq_max, nfreq, use_centered_grid) + return self.fields.add_dft_fields(components, where, freq_min, freq_max, freqs, nfreq, use_centered_grid) def output_dft(self, dft_fields, fname): if self.fields is None: @@ -1664,6 +1670,7 @@ def get_dft_data(self, dft_chunk): mp._get_dft_data(dft_chunk, arr) return arr + # TODO: overload add_near2far, _add_near2far def add_near2far(self, fcen, df, nfreq, *near2fars, **kwargs): nperiods = kwargs.get('nperiods', 1) n2f = DftNear2Far(self._add_near2far, [fcen, df, nfreq, nperiods, near2fars]) @@ -1675,6 +1682,7 @@ def _add_near2far(self, fcen, df, nfreq, nperiods, near2fars): self.init_sim() return self._add_fluxish_stuff(self.fields.add_dft_near2far, fcen, df, nfreq, near2fars, nperiods) + # TODO: overload add_energy, _add_energy def add_energy(self, fcen, df, nfreq, *energys): en = DftEnergy(self._add_energy, [fcen, df, nfreq, energys]) self.dft_objects.append(en) @@ -1685,6 +1693,7 @@ def _add_energy(self, fcen, df, nfreq, energys): self.init_sim() return self._add_fluxish_stuff(self.fields.add_dft_energy, fcen, df, nfreq, energys) + # TODO: check if this needs to be changed def _display_energy(self, name, func, energys): if energys: freqs = get_energy_freqs(energys[0]) @@ -1770,6 +1779,7 @@ def load_minus_near2far_data(self, n2f, n2fdata): self.load_near2far_data(n2f, n2fdata) n2f.scale_dfts(complex(-1.0)) + # TODO: overload add_force, _add_force def add_force(self, fcen, df, nfreq, *forces): force = DftForce(self._add_force, [fcen, df, nfreq, forces]) self.dft_objects.append(force) @@ -1780,6 +1790,7 @@ def _add_force(self, fcen, df, nfreq, forces): self.init_sim() return self._add_fluxish_stuff(self.fields.add_dft_force, fcen, df, nfreq, forces) + # TODO: check if display_forces needs to be changed def display_forces(self, *forces): force_freqs = get_force_freqs(forces[0]) display_csv(self, 'force', zip(force_freqs, *[get_forces(f) for f in forces])) @@ -1812,6 +1823,7 @@ def load_minus_force_data(self, force, fdata): self.load_force_data(force, fdata) force.scale_dfts(complex(-1.0)) + # TODO: overload add_flux, _add_flux def add_flux(self, fcen, df, nfreq, *fluxes): flux = DftFlux(self._add_flux, [fcen, df, nfreq, fluxes]) self.dft_objects.append(flux) @@ -1822,6 +1834,7 @@ def _add_flux(self, fcen, df, nfreq, fluxes): self.init_sim() return self._add_fluxish_stuff(self.fields.add_dft_flux, fcen, df, nfreq, fluxes) + # TODO: overload add_mode_monitor, _add_mode_monitor def add_mode_monitor(self, fcen, df, nfreq, *fluxes): flux = DftFlux(self._add_mode_monitor, [fcen, df, nfreq, fluxes]) self.dft_objects.append(flux) @@ -1841,6 +1854,7 @@ def _add_mode_monitor(self, fcen, df, nfreq, fluxes): return self.fields.add_mode_monitor(d, v.swigobj, fcen - df / 2, fcen + df / 2, nfreq) + # TODO: overload add_eigenmode def add_eigenmode(self, fcen, df, nfreq, *fluxes): warnings.warn('add_eigenmode is deprecated. Please use add_mode_monitor instead.', DeprecationWarning) return self.add_mode_monitor(fcen, df, nfreq, *fluxes) @@ -2945,6 +2959,7 @@ def Ldos(fcen, df, nfreq): return mp._dft_ldos(fcen - df / 2, fcen + df / 2, nfreq) +# TODO: overload dft_ldos def dft_ldos(fcen=None, df=None, nfreq=None, ldos=None): if ldos is None: if fcen is None or df is None or nfreq is None: @@ -2999,6 +3014,7 @@ def get_near2far_freqs(f): def scale_energy_fields(s, ef): + # TODO: check if this is a typo df.scale_dfts(s) From 082247da1d72ab72094c75407172f90596e7530d Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Mon, 9 Mar 2020 16:14:48 -0600 Subject: [PATCH 17/35] * Created add_*_uneven functions * Wrapped add_* functions around add_*_uneven functions * Changed parameters of _add_* functions --- python/simulation.py | 133 ++++++++++++++++++++++++++----------------- 1 file changed, 81 insertions(+), 52 deletions(-) diff --git a/python/simulation.py b/python/simulation.py index 56a2f0186..bbf122d77 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -1626,17 +1626,9 @@ def _evaluate_dft_objects(self): if dft.swigobj is None: dft.swigobj = dft.func(*dft.args) - # TODO: overload add_dft_fields, _add_dft_fields - def add_dft_fields(self, components, freq_min=None, freq_max=None, freqs=None, nfreq, where=None, center=None, size=None, yee_grid=False): - if freqs is None: - if freq_min is None or freq_max is None: - raise ValueError("Either freqs or freq_min and freq_max is required for add_dft_fields") - else: - freqs = np.linspace(start=freq_min, end=freq_max, num=nfreq) - else: - freqs.sort() - freq_min = freqs[0] - freq_max = freqs[-1] + # TODO: make add_dft_fields_uneven, wrap add_dft_fields around it, change parameters of _add_dft_fields + def add_dft_fields_uneven(self, components, freqs, nfreq, where=None, center=None, size=None, yee_grid=False): + freqs.sort() center_v3 = Vector3(*center) if center is not None else None size_v3 = Vector3(*size) if size is not None else None use_centered_grid = not yee_grid @@ -1644,14 +1636,18 @@ def add_dft_fields(self, components, freq_min=None, freq_max=None, freqs=None, n self.dft_objects.append(dftf) return dftf - def _add_dft_fields(self, components, where, center, size, freq_min, freq_max, freqs, nfreq, use_centered_grid): + def add_dft_fields(self, components, freq_min, freq_max, nfreq, where=None, center=None, size=None, yee_grid=False): + freqs = np.linspace(start=freq_min, end=freq_max, num=nfreq).tolist() + return self.add_dft_fields_uneven(components, freqs, nfreq, where, center, size, yee_grid) + + def _add_dft_fields(self, components, where, center, size, freqs, nfreq, use_centered_grid): if self.fields is None: self.init_sim() try: where = self._volume_from_kwargs(where, center, size) except ValueError: where = self.fields.total_volume() - return self.fields.add_dft_fields(components, where, freq_min, freq_max, freqs, nfreq, use_centered_grid) + return self.fields.add_dft_fields(components, where, freqs, nfreq, use_centered_grid) def output_dft(self, dft_fields, fname): if self.fields is None: @@ -1670,28 +1666,36 @@ def get_dft_data(self, dft_chunk): mp._get_dft_data(dft_chunk, arr) return arr - # TODO: overload add_near2far, _add_near2far - def add_near2far(self, fcen, df, nfreq, *near2fars, **kwargs): + # TODO: make add_near2far_uneven, wrap add_near2far around it, change parameters of _add_near2far + def add_near2far_uneven(self, freqs, nfreq, *near2fars, **kwargs): nperiods = kwargs.get('nperiods', 1) - n2f = DftNear2Far(self._add_near2far, [fcen, df, nfreq, nperiods, near2fars]) + n2f = DftNear2Far(self._add_near2far, [freqs, nfreq, nperiods, near2fars]) self.dft_objects.append(n2f) return n2f - def _add_near2far(self, fcen, df, nfreq, nperiods, near2fars): + def add_near2far(self, fcen, df, nfreq, *near2fars, **kwargs): + freqs = np.linspace(fcen - df * (nfreq - 1) / 2, fcen + df * (nfreq - 1) / 2, nfreq).tolist() + return self.add_near2far_uneven(freqs, nfreq, *near2fars, **kwargs) + + def _add_near2far(self, freqs, nfreq, nperiods, near2fars): if self.fields is None: self.init_sim() - return self._add_fluxish_stuff(self.fields.add_dft_near2far, fcen, df, nfreq, near2fars, nperiods) + return self._add_fluxish_stuff(self.fields.add_dft_near2far, freqs, nfreq, near2fars, nperiods) - # TODO: overload add_energy, _add_energy - def add_energy(self, fcen, df, nfreq, *energys): - en = DftEnergy(self._add_energy, [fcen, df, nfreq, energys]) + # TODO: make add_energy_uneven, wrap add_energy around it, change parameters of _add_energy + def add_energy_uneven(self, freqs, nfreq, *energys): + en = DftEnergy(self._add_energy, [freqs, nfreq, energys]) self.dft_objects.append(en) return en - def _add_energy(self, fcen, df, nfreq, energys): + def add_energy(self, fcen, df, nfreq, *energys): + freqs = np.linspace(fcen - df * (nfreq - 1) / 2, fcen + df * (nfreq - 1) / 2, nfreq).tolist() + return self.add_energy_uneven(freqs, nfreq, *energys) + + def _add_energy(self, freqs, nfreq, energys): if self.fields is None: self.init_sim() - return self._add_fluxish_stuff(self.fields.add_dft_energy, fcen, df, nfreq, energys) + return self._add_fluxish_stuff(self.fields.add_dft_energy, freqs, nfreq, energys) # TODO: check if this needs to be changed def _display_energy(self, name, func, energys): @@ -1779,16 +1783,20 @@ def load_minus_near2far_data(self, n2f, n2fdata): self.load_near2far_data(n2f, n2fdata) n2f.scale_dfts(complex(-1.0)) - # TODO: overload add_force, _add_force - def add_force(self, fcen, df, nfreq, *forces): - force = DftForce(self._add_force, [fcen, df, nfreq, forces]) + # TODO: make add_force_uneven, wrap add_force around it, change parameters of _add_force + def add_force_uneven(self, freqs, nfreq, *forces): + force = DftForce(self._add_force, [freqs, nfreq, forces]) self.dft_objects.append(force) return force - def _add_force(self, fcen, df, nfreq, forces): + def add_force(self, fcen, df, nfreq, *forces): + freqs = np.linspace(fcen - df * (nfreq - 1) / 2, fcen + df * (nfreq - 1) / 2, nfreq).tolist() + return self.add_force_uneven(freqs, nfreq, *forces) + + def _add_force(self, freqs, nfreq, forces): if self.fields is None: self.init_sim() - return self._add_fluxish_stuff(self.fields.add_dft_force, fcen, df, nfreq, forces) + return self._add_fluxish_stuff(self.fields.add_dft_force, freqs, nfreq, forces) # TODO: check if display_forces needs to be changed def display_forces(self, *forces): @@ -1823,24 +1831,32 @@ def load_minus_force_data(self, force, fdata): self.load_force_data(force, fdata) force.scale_dfts(complex(-1.0)) - # TODO: overload add_flux, _add_flux - def add_flux(self, fcen, df, nfreq, *fluxes): - flux = DftFlux(self._add_flux, [fcen, df, nfreq, fluxes]) + # TODO: make add_flux_uneven, wrap add_flux around it, change parameters of _add_flux + def add_flux_uneven(self, freqs, nfreq, *fluxes): + flux = DftFlux(self._add_flux, [freqs, nfreq, fluxes]) self.dft_objects.append(flux) return flux - def _add_flux(self, fcen, df, nfreq, fluxes): + def add_flux(self, fcen, df, nfreq, *fluxes): + freqs = np.linspace(fcen - df * (nfreq - 1) / 2, fcen + df * (nfreq - 1) / 2, nfreq).tolist() + return self.add_flux_uneven(freqs, nfreq, *fluxes) + + def _add_flux(self, freqs, nfreq, fluxes): if self.fields is None: self.init_sim() - return self._add_fluxish_stuff(self.fields.add_dft_flux, fcen, df, nfreq, fluxes) + return self._add_fluxish_stuff(self.fields.add_dft_flux, freqs, nfreq, fluxes) - # TODO: overload add_mode_monitor, _add_mode_monitor - def add_mode_monitor(self, fcen, df, nfreq, *fluxes): - flux = DftFlux(self._add_mode_monitor, [fcen, df, nfreq, fluxes]) + # TODO: make add_mode_monitor_uneven, wrap add_mode_monitor around it, change parameters of _add_mode_monitor + def add_mode_monitor_uneven(self, freqs, nfreq, *fluxes): + flux = DftFlux(self._add_mode_monitor, [freqs, nfreq, fluxes]) self.dft_objects.append(flux) return flux - def _add_mode_monitor(self, fcen, df, nfreq, fluxes): + def add_mode_monitor(self, fcen, df, nfreq, *fluxes): + freqs = np.linspace(fcen - df * (nfreq - 1) / 2, fcen + df * (nfreq - 1) / 2, nfreq).tolist() + return self.add_mode_monitor_uneven(freqs, nfreq, *fluxes) + + def _add_mode_monitor(self, freqs, nfreq, fluxes): if self.fields is None: self.init_sim() @@ -1852,9 +1868,8 @@ def _add_mode_monitor(self, fcen, df, nfreq, fluxes): d0 = region.direction d = self.fields.normal_direction(v.swigobj) if d0 < 0 else d0 - return self.fields.add_mode_monitor(d, v.swigobj, fcen - df / 2, fcen + df / 2, nfreq) + return self.fields.add_mode_monitor(d, v.swigobj, freqs, nfreq) - # TODO: overload add_eigenmode def add_eigenmode(self, fcen, df, nfreq, *fluxes): warnings.warn('add_eigenmode is deprecated. Please use add_mode_monitor instead.', DeprecationWarning) return self.add_mode_monitor(fcen, df, nfreq, *fluxes) @@ -1950,7 +1965,8 @@ def solve_cw(self, tol=1e-8, maxiters=10000, L=2): self._evaluate_dft_objects() return self.fields.solve_cw(tol, maxiters, L) - def _add_fluxish_stuff(self, add_dft_stuff, fcen, df, nfreq, stufflist, *args): + # TODO: change parameters of _add_fluxish_stuff + def _add_fluxish_stuff(self, add_dft_stuff, freqs, nfreq, stufflist, *args): vol_list = None for s in stufflist: @@ -1963,7 +1979,7 @@ def _add_fluxish_stuff(self, add_dft_stuff, fcen, df, nfreq, stufflist, *args): is_cylindrical=self.is_cylindrical).swigobj vol_list = mp.make_volume_list(v2, c, s.weight, vol_list) - stuff = add_dft_stuff(vol_list, fcen - df / 2, fcen + df / 2, nfreq, *args) + stuff = add_dft_stuff(vol_list, freqs, nfreq, *args) vol_list.__swig_destroy__(vol_list) return stuff @@ -2959,12 +2975,12 @@ def Ldos(fcen, df, nfreq): return mp._dft_ldos(fcen - df / 2, fcen + df / 2, nfreq) -# TODO: overload dft_ldos -def dft_ldos(fcen=None, df=None, nfreq=None, ldos=None): +# TODO: make dft_ldos_uneven, wrap dft_ldos wrap around it +def dft_ldos_uneven(freqs=None, nfreq=None, ldos=None): if ldos is None: - if fcen is None or df is None or nfreq is None: - raise ValueError("Either fcen, df, and nfreq, or an Ldos is required for dft_ldos") - ldos = mp._dft_ldos(fcen - df / 2, fcen + df / 2, nfreq) + if freqs is None or nfreq is None: + raise ValueError("Either freqs and nfreq or an Ldos is required for dft_ldos_uneven") + ldos = mp._dft_ldos(freqs, nfreq) def _ldos(sim, todo): if todo == 'step': @@ -2977,12 +2993,21 @@ def _ldos(sim, todo): return _ldos +def dft_ldos(fcen=None, df=None, nfreq=None, ldos=None): + if ldos is None: + if fcen is None or df is None or nfreq is None: + raise ValueError("Either fcen, df, and nfreq, or an Ldos is required for dft_ldos") + freqs = np.linspace(fcen - df * (nfreq - 1) / 2, fcen + df * (nfreq - 1) / 2, nfreq).tolist() + return dft_ldos_uneven(freqs, nfreq, ldos) + + def scale_flux_fields(s, flux): flux.scale_dfts(s) +# TODO: change this def get_flux_freqs(f): - return np.linspace(f.freq_min, f.freq_min + f.dfreq * f.Nfreq, num=f.Nfreq, endpoint=False).tolist() + return f.freqs def get_fluxes(f): @@ -2993,12 +3018,14 @@ def scale_force_fields(s, force): force.scale_dfts(s) +# TODO: change this def get_eigenmode_freqs(f): - return np.linspace(f.freq_min, f.freq_min + f.dfreq * f.Nfreq, num=f.Nfreq, endpoint=False).tolist() + return f.freqs +# TODO: change this def get_force_freqs(f): - return np.linspace(f.freq_min, f.freq_min + f.dfreq * f.Nfreq, num=f.Nfreq, endpoint=False).tolist() + return f.freqs def get_forces(f): @@ -3009,17 +3036,19 @@ def scale_near2far_fields(s, n2f): n2f.scale_dfts(s) +# TODO: change this def get_near2far_freqs(f): - return np.linspace(f.freq_min, f.freq_min + f.dfreq * f.Nfreq, num=f.Nfreq, endpoint=False).tolist() + return f.freqs def scale_energy_fields(s, ef): # TODO: check if this is a typo - df.scale_dfts(s) + ef.scale_dfts(s) +# TODO: change this def get_energy_freqs(f): - return np.linspace(f.freq_min, f.freq_min + f.dfreq * f.Nfreq, num=f.Nfreq, endpoint=False).tolist() + return f.freqs def get_electric_energy(f): From 9f986bf8ca87eff58a094decd81c4a9eba1be58d Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Mon, 9 Mar 2020 17:15:54 -0600 Subject: [PATCH 18/35] Make sure freqs is sorted --- python/simulation.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/python/simulation.py b/python/simulation.py index bbf122d77..fbcd55fa9 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -1668,6 +1668,7 @@ def get_dft_data(self, dft_chunk): # TODO: make add_near2far_uneven, wrap add_near2far around it, change parameters of _add_near2far def add_near2far_uneven(self, freqs, nfreq, *near2fars, **kwargs): + freqs.sort() nperiods = kwargs.get('nperiods', 1) n2f = DftNear2Far(self._add_near2far, [freqs, nfreq, nperiods, near2fars]) self.dft_objects.append(n2f) @@ -1684,6 +1685,7 @@ def _add_near2far(self, freqs, nfreq, nperiods, near2fars): # TODO: make add_energy_uneven, wrap add_energy around it, change parameters of _add_energy def add_energy_uneven(self, freqs, nfreq, *energys): + freqs.sort() en = DftEnergy(self._add_energy, [freqs, nfreq, energys]) self.dft_objects.append(en) return en @@ -1785,6 +1787,7 @@ def load_minus_near2far_data(self, n2f, n2fdata): # TODO: make add_force_uneven, wrap add_force around it, change parameters of _add_force def add_force_uneven(self, freqs, nfreq, *forces): + freqs.sort() force = DftForce(self._add_force, [freqs, nfreq, forces]) self.dft_objects.append(force) return force @@ -1833,6 +1836,7 @@ def load_minus_force_data(self, force, fdata): # TODO: make add_flux_uneven, wrap add_flux around it, change parameters of _add_flux def add_flux_uneven(self, freqs, nfreq, *fluxes): + freqs.sort() flux = DftFlux(self._add_flux, [freqs, nfreq, fluxes]) self.dft_objects.append(flux) return flux @@ -1848,6 +1852,7 @@ def _add_flux(self, freqs, nfreq, fluxes): # TODO: make add_mode_monitor_uneven, wrap add_mode_monitor around it, change parameters of _add_mode_monitor def add_mode_monitor_uneven(self, freqs, nfreq, *fluxes): + freqs.sort() flux = DftFlux(self._add_mode_monitor, [freqs, nfreq, fluxes]) self.dft_objects.append(flux) return flux @@ -2980,6 +2985,7 @@ def dft_ldos_uneven(freqs=None, nfreq=None, ldos=None): if ldos is None: if freqs is None or nfreq is None: raise ValueError("Either freqs and nfreq or an Ldos is required for dft_ldos_uneven") + freqs.sort() ldos = mp._dft_ldos(freqs, nfreq) def _ldos(sim, todo): From 1f972bf0c9d55ef2276d737193c1c8aebedf3eba Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Mon, 9 Mar 2020 19:20:20 -0600 Subject: [PATCH 19/35] Removed (nfreq - 1) term from linspace functions; had misunderstood what quantity df was. --- python/simulation.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/python/simulation.py b/python/simulation.py index fbcd55fa9..2a037fb66 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -1675,7 +1675,7 @@ def add_near2far_uneven(self, freqs, nfreq, *near2fars, **kwargs): return n2f def add_near2far(self, fcen, df, nfreq, *near2fars, **kwargs): - freqs = np.linspace(fcen - df * (nfreq - 1) / 2, fcen + df * (nfreq - 1) / 2, nfreq).tolist() + freqs = np.linspace(fcen - df / 2, fcen + df / 2, nfreq).tolist() return self.add_near2far_uneven(freqs, nfreq, *near2fars, **kwargs) def _add_near2far(self, freqs, nfreq, nperiods, near2fars): @@ -1691,7 +1691,7 @@ def add_energy_uneven(self, freqs, nfreq, *energys): return en def add_energy(self, fcen, df, nfreq, *energys): - freqs = np.linspace(fcen - df * (nfreq - 1) / 2, fcen + df * (nfreq - 1) / 2, nfreq).tolist() + freqs = np.linspace(fcen - df / 2, fcen + df / 2, nfreq).tolist() return self.add_energy_uneven(freqs, nfreq, *energys) def _add_energy(self, freqs, nfreq, energys): @@ -1793,7 +1793,7 @@ def add_force_uneven(self, freqs, nfreq, *forces): return force def add_force(self, fcen, df, nfreq, *forces): - freqs = np.linspace(fcen - df * (nfreq - 1) / 2, fcen + df * (nfreq - 1) / 2, nfreq).tolist() + freqs = np.linspace(fcen - df / 2, fcen + df / 2, nfreq).tolist() return self.add_force_uneven(freqs, nfreq, *forces) def _add_force(self, freqs, nfreq, forces): @@ -1842,7 +1842,7 @@ def add_flux_uneven(self, freqs, nfreq, *fluxes): return flux def add_flux(self, fcen, df, nfreq, *fluxes): - freqs = np.linspace(fcen - df * (nfreq - 1) / 2, fcen + df * (nfreq - 1) / 2, nfreq).tolist() + freqs = np.linspace(fcen - df / 2, fcen + df / 2, nfreq).tolist() return self.add_flux_uneven(freqs, nfreq, *fluxes) def _add_flux(self, freqs, nfreq, fluxes): @@ -1858,7 +1858,7 @@ def add_mode_monitor_uneven(self, freqs, nfreq, *fluxes): return flux def add_mode_monitor(self, fcen, df, nfreq, *fluxes): - freqs = np.linspace(fcen - df * (nfreq - 1) / 2, fcen + df * (nfreq - 1) / 2, nfreq).tolist() + freqs = np.linspace(fcen - df / 2, fcen + df / 2, nfreq).tolist() return self.add_mode_monitor_uneven(freqs, nfreq, *fluxes) def _add_mode_monitor(self, freqs, nfreq, fluxes): @@ -3003,7 +3003,7 @@ def dft_ldos(fcen=None, df=None, nfreq=None, ldos=None): if ldos is None: if fcen is None or df is None or nfreq is None: raise ValueError("Either fcen, df, and nfreq, or an Ldos is required for dft_ldos") - freqs = np.linspace(fcen - df * (nfreq - 1) / 2, fcen + df * (nfreq - 1) / 2, nfreq).tolist() + freqs = np.linspace(fcen - df / 2, fcen + df / 2, nfreq).tolist() return dft_ldos_uneven(freqs, nfreq, ldos) From b10366c97d1d7024e41aa1be1fdcf1530585bdd2 Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Mon, 9 Mar 2020 19:22:20 -0600 Subject: [PATCH 20/35] Replacement of get_omega_min and get_domega with get_omegas --- python/meep.i | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/python/meep.i b/python/meep.i index 3c8498c1e..86e433d38 100644 --- a/python/meep.i +++ b/python/meep.i @@ -1346,11 +1346,8 @@ void _get_eigenmode(meep::fields *f, double omega_src, meep::direction d, const // Make omega members of meep::dft_ldos available as 'freq' in python %extend meep::dft_ldos { - double get_omega_min() { - return $self->omega_min; - } - double get_domega() { - return $self->domega; + double* get_omegas() { + return $self->omegas; } int get_Nomega() { return $self->Nomega; From e075a3401507a2e37f7a265e0b8a0e422c2270be Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Mon, 9 Mar 2020 22:16:02 -0600 Subject: [PATCH 21/35] Alterations for dft_ldos in python/meep.i --- python/meep.i | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/python/meep.i b/python/meep.i index 86e433d38..b9300bd3f 100644 --- a/python/meep.i +++ b/python/meep.i @@ -1360,12 +1360,10 @@ void _get_eigenmode(meep::fields *f, double omega_src, meep::direction d, const stop = start + (self.domega / (2 * math.pi)) * self.Nomega return np.linspace(start, stop, num=self.Nomega, endpoint=False).tolist() - __swig_getmethods__["freq_min"] = get_omega_min + __swig_getmethods__["freqs"] = get_omegas __swig_getmethods__["nfreq"] = get_Nomega - __swig_getmethods__["dfreq"] = get_domega - if _newclass: freq_min = property(get_omega_min) + if _newclass: freqs = property(get_omegas) if _newclass: nfreq = property(get_Nomega) - if _newclass: dfreq = property(get_domega) %} } From d7c8a7d611476a6a002087dd43b5443b78a480d2 Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Tue, 10 Mar 2020 21:38:02 -0600 Subject: [PATCH 22/35] Change because of differing number of arguments --- python/simulation.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/simulation.py b/python/simulation.py index 2a037fb66..8ff31ef08 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -332,8 +332,8 @@ class DftFlux(DftObj): def __init__(self, func, args): super(DftFlux, self).__init__(func, args) - self.nfreqs = args[2] - self.regions = args[3] + self.nfreqs = args[1] + self.regions = args[2] self.num_components = 4 @property From bee4c0bb69d92d6d8878b96815f59f8e05742168 Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Wed, 11 Mar 2020 18:44:52 -0600 Subject: [PATCH 23/35] Additional attempt #1 at typemaps in SWIG --- python/meep.i | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/python/meep.i b/python/meep.i index ea566a2be..762b75b85 100644 --- a/python/meep.i +++ b/python/meep.i @@ -966,6 +966,24 @@ meep::volume_list *make_volume_list(const meep::volume &v, int c, // typemaps needed for add_dft_fields //-------------------------------------------------- +%typecheck(SWIG_TYPECHECK_POINTER, fragment="NumPy_Fragments") double* freqs { + $1 = is_array($input); +} + +%typemap(in, fragment="NumPy_Macros") double* freqs { + $1 = (double *)array_data($input); +} + +/* +%typecheck(SWIG_TYPECHECK_POINTER, fragment="NumPy_Fragments") double* fs { + $1 = is_array($input); +} + +%typemap(in, fragment="NumPy_Macros") double* fs { + $1 = (double *)array_data($input); +} +*/ + %typecheck(SWIG_TYPECHECK_POINTER) const volume where { int py_material = PyObject_IsInstance($input, py_volume_object()); $1 = py_material; From c7818504ad1e6685876bfe64b870efa2beaa43c3 Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Wed, 11 Mar 2020 18:46:41 -0600 Subject: [PATCH 24/35] cleanup --- python/meep.i | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/python/meep.i b/python/meep.i index 762b75b85..89f9aaa9c 100644 --- a/python/meep.i +++ b/python/meep.i @@ -974,16 +974,6 @@ meep::volume_list *make_volume_list(const meep::volume &v, int c, $1 = (double *)array_data($input); } -/* -%typecheck(SWIG_TYPECHECK_POINTER, fragment="NumPy_Fragments") double* fs { - $1 = is_array($input); -} - -%typemap(in, fragment="NumPy_Macros") double* fs { - $1 = (double *)array_data($input); -} -*/ - %typecheck(SWIG_TYPECHECK_POINTER) const volume where { int py_material = PyObject_IsInstance($input, py_volume_object()); $1 = py_material; From 4dd0683e3301067f7f0857251d4e578e2766d369 Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Wed, 11 Mar 2020 19:58:01 -0600 Subject: [PATCH 25/35] this is certainly wrong --- python/meep.i | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/python/meep.i b/python/meep.i index 89f9aaa9c..4fbc70077 100644 --- a/python/meep.i +++ b/python/meep.i @@ -974,6 +974,17 @@ meep::volume_list *make_volume_list(const meep::volume &v, int c, $1 = (double *)array_data($input); } +%typecheck(SWIG_TYPECHECK_INTEGER) int Nfreqs { + $1 = PyInteger_Check($input); +} + +%typemap(in) int Nfreqs { + int py_nfreqs = PyInteger_Check($input); + $1 = py_nfreqs; +} + +%apply (double *IN_ARRAY1, int DIM1) {(double *freqs, int Nfreqs)}; + %typecheck(SWIG_TYPECHECK_POINTER) const volume where { int py_material = PyObject_IsInstance($input, py_volume_object()); $1 = py_material; From df25758dc08983e51e46ee74a1a691f3697484e8 Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Wed, 11 Mar 2020 21:22:36 -0600 Subject: [PATCH 26/35] Remove extra tolist() --- python/simulation.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/python/simulation.py b/python/simulation.py index 8ff31ef08..6022c8bb0 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -1637,7 +1637,7 @@ def add_dft_fields_uneven(self, components, freqs, nfreq, where=None, center=Non return dftf def add_dft_fields(self, components, freq_min, freq_max, nfreq, where=None, center=None, size=None, yee_grid=False): - freqs = np.linspace(start=freq_min, end=freq_max, num=nfreq).tolist() + freqs = np.linspace(start=freq_min, end=freq_max, num=nfreq) return self.add_dft_fields_uneven(components, freqs, nfreq, where, center, size, yee_grid) def _add_dft_fields(self, components, where, center, size, freqs, nfreq, use_centered_grid): @@ -1675,7 +1675,7 @@ def add_near2far_uneven(self, freqs, nfreq, *near2fars, **kwargs): return n2f def add_near2far(self, fcen, df, nfreq, *near2fars, **kwargs): - freqs = np.linspace(fcen - df / 2, fcen + df / 2, nfreq).tolist() + freqs = np.linspace(fcen - df / 2, fcen + df / 2, nfreq) return self.add_near2far_uneven(freqs, nfreq, *near2fars, **kwargs) def _add_near2far(self, freqs, nfreq, nperiods, near2fars): @@ -1691,7 +1691,7 @@ def add_energy_uneven(self, freqs, nfreq, *energys): return en def add_energy(self, fcen, df, nfreq, *energys): - freqs = np.linspace(fcen - df / 2, fcen + df / 2, nfreq).tolist() + freqs = np.linspace(fcen - df / 2, fcen + df / 2, nfreq) return self.add_energy_uneven(freqs, nfreq, *energys) def _add_energy(self, freqs, nfreq, energys): @@ -1793,7 +1793,7 @@ def add_force_uneven(self, freqs, nfreq, *forces): return force def add_force(self, fcen, df, nfreq, *forces): - freqs = np.linspace(fcen - df / 2, fcen + df / 2, nfreq).tolist() + freqs = np.linspace(fcen - df / 2, fcen + df / 2, nfreq) return self.add_force_uneven(freqs, nfreq, *forces) def _add_force(self, freqs, nfreq, forces): @@ -1842,7 +1842,7 @@ def add_flux_uneven(self, freqs, nfreq, *fluxes): return flux def add_flux(self, fcen, df, nfreq, *fluxes): - freqs = np.linspace(fcen - df / 2, fcen + df / 2, nfreq).tolist() + freqs = np.linspace(fcen - df / 2, fcen + df / 2, nfreq) return self.add_flux_uneven(freqs, nfreq, *fluxes) def _add_flux(self, freqs, nfreq, fluxes): @@ -1858,7 +1858,7 @@ def add_mode_monitor_uneven(self, freqs, nfreq, *fluxes): return flux def add_mode_monitor(self, fcen, df, nfreq, *fluxes): - freqs = np.linspace(fcen - df / 2, fcen + df / 2, nfreq).tolist() + freqs = np.linspace(fcen - df / 2, fcen + df / 2, nfreq) return self.add_mode_monitor_uneven(freqs, nfreq, *fluxes) def _add_mode_monitor(self, freqs, nfreq, fluxes): @@ -3003,7 +3003,7 @@ def dft_ldos(fcen=None, df=None, nfreq=None, ldos=None): if ldos is None: if fcen is None or df is None or nfreq is None: raise ValueError("Either fcen, df, and nfreq, or an Ldos is required for dft_ldos") - freqs = np.linspace(fcen - df / 2, fcen + df / 2, nfreq).tolist() + freqs = np.linspace(fcen - df / 2, fcen + df / 2, nfreq) return dft_ldos_uneven(freqs, nfreq, ldos) From db61366c82d3133661f9e68d550f38643bdbb27a Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Wed, 11 Mar 2020 23:04:50 -0600 Subject: [PATCH 27/35] python/meep.i still doesn't seem to be working --- python/meep.i | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/meep.i b/python/meep.i index 4fbc70077..690838699 100644 --- a/python/meep.i +++ b/python/meep.i @@ -965,7 +965,7 @@ meep::volume_list *make_volume_list(const meep::volume &v, int c, //-------------------------------------------------- // typemaps needed for add_dft_fields //-------------------------------------------------- - +/* %typecheck(SWIG_TYPECHECK_POINTER, fragment="NumPy_Fragments") double* freqs { $1 = is_array($input); } @@ -982,7 +982,7 @@ meep::volume_list *make_volume_list(const meep::volume &v, int c, int py_nfreqs = PyInteger_Check($input); $1 = py_nfreqs; } - +*/ %apply (double *IN_ARRAY1, int DIM1) {(double *freqs, int Nfreqs)}; %typecheck(SWIG_TYPECHECK_POINTER) const volume where { From 27f453621bd8324808d09b568a4fe2ac67639b4e Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Sun, 22 Mar 2020 18:02:30 -0600 Subject: [PATCH 28/35] Revert python/meep.i and python/simulation.py to master state --- python/meep.i | 32 +++------- python/simulation.py | 143 ++++++++++++++----------------------------- 2 files changed, 55 insertions(+), 120 deletions(-) diff --git a/python/meep.i b/python/meep.i index 5b6c75d82..06e4ad345 100644 --- a/python/meep.i +++ b/python/meep.i @@ -965,25 +965,6 @@ meep::volume_list *make_volume_list(const meep::volume &v, int c, //-------------------------------------------------- // typemaps needed for add_dft_fields //-------------------------------------------------- -/* -%typecheck(SWIG_TYPECHECK_POINTER, fragment="NumPy_Fragments") double* freqs { - $1 = is_array($input); -} - -%typemap(in, fragment="NumPy_Macros") double* freqs { - $1 = (double *)array_data($input); -} - -%typecheck(SWIG_TYPECHECK_INTEGER) int Nfreqs { - $1 = PyInteger_Check($input); -} - -%typemap(in) int Nfreqs { - int py_nfreqs = PyInteger_Check($input); - $1 = py_nfreqs; -} -*/ -%apply (double *IN_ARRAY1, int DIM1) {(double *freqs, int Nfreqs)}; %typecheck(SWIG_TYPECHECK_POINTER) const volume where { int py_material = PyObject_IsInstance($input, py_volume_object()); @@ -1365,8 +1346,11 @@ void _get_eigenmode(meep::fields *f, double omega_src, meep::direction d, const // Make omega members of meep::dft_ldos available as 'freq' in python %extend meep::dft_ldos { - double* get_omegas() { - return $self->omegas; + double get_omega_min() { + return $self->omega_min; + } + double get_domega() { + return $self->domega; } int get_Nomega() { return $self->Nomega; @@ -1379,10 +1363,12 @@ void _get_eigenmode(meep::fields *f, double omega_src, meep::direction d, const stop = start + (self.domega / (2 * math.pi)) * self.Nomega return np.linspace(start, stop, num=self.Nomega, endpoint=False).tolist() - __swig_getmethods__["freqs"] = get_omegas + __swig_getmethods__["freq_min"] = get_omega_min __swig_getmethods__["nfreq"] = get_Nomega - if _newclass: freqs = property(get_omegas) + __swig_getmethods__["dfreq"] = get_domega + if _newclass: freq_min = property(get_omega_min) if _newclass: nfreq = property(get_Nomega) + if _newclass: dfreq = property(get_domega) %} } diff --git a/python/simulation.py b/python/simulation.py index 6022c8bb0..a58e42430 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -316,8 +316,12 @@ def remove(self): return self.swigobj_attr('remove') @property - def freqs(self): - return self.swigobj_attr('freqs') + def freq_min(self): + return self.swigobj_attr('freq_min') + + @property + def dfreq(self): + return self.swigobj_attr('dfreq') @property def Nfreq(self): @@ -332,8 +336,8 @@ class DftFlux(DftObj): def __init__(self, func, args): super(DftFlux, self).__init__(func, args) - self.nfreqs = args[1] - self.regions = args[2] + self.nfreqs = args[2] + self.regions = args[3] self.num_components = 4 @property @@ -717,8 +721,8 @@ def _check_material_frequencies(self): dft_freqs = [] for dftf in self.dft_objects: - dft_freqs.append(dftf.freqs[0]) - dft_freqs.append(dftf.freqs[-1]) + dft_freqs.append(dftf.freq_min) + dft_freqs.append(dftf.freq_min + dftf.Nfreq * dftf.dfreq) warn_src = ('Note: your sources include frequencies outside the range of validity of the ' + 'material models. This is fine as long as you eventually only look at outputs ' + @@ -1626,28 +1630,22 @@ def _evaluate_dft_objects(self): if dft.swigobj is None: dft.swigobj = dft.func(*dft.args) - # TODO: make add_dft_fields_uneven, wrap add_dft_fields around it, change parameters of _add_dft_fields - def add_dft_fields_uneven(self, components, freqs, nfreq, where=None, center=None, size=None, yee_grid=False): - freqs.sort() + def add_dft_fields(self, components, freq_min, freq_max, nfreq, where=None, center=None, size=None, yee_grid=False): center_v3 = Vector3(*center) if center is not None else None size_v3 = Vector3(*size) if size is not None else None use_centered_grid = not yee_grid - dftf = DftFields(self._add_dft_fields, [components, where, center_v3, size_v3, freqs, nfreq, use_centered_grid]) + dftf = DftFields(self._add_dft_fields, [components, where, center_v3, size_v3, freq_min, freq_max, nfreq, use_centered_grid]) self.dft_objects.append(dftf) return dftf - def add_dft_fields(self, components, freq_min, freq_max, nfreq, where=None, center=None, size=None, yee_grid=False): - freqs = np.linspace(start=freq_min, end=freq_max, num=nfreq) - return self.add_dft_fields_uneven(components, freqs, nfreq, where, center, size, yee_grid) - - def _add_dft_fields(self, components, where, center, size, freqs, nfreq, use_centered_grid): + def _add_dft_fields(self, components, where, center, size, freq_min, freq_max, nfreq, use_centered_grid): if self.fields is None: self.init_sim() try: where = self._volume_from_kwargs(where, center, size) except ValueError: where = self.fields.total_volume() - return self.fields.add_dft_fields(components, where, freqs, nfreq, use_centered_grid) + return self.fields.add_dft_fields(components, where, freq_min, freq_max, nfreq, use_centered_grid) def output_dft(self, dft_fields, fname): if self.fields is None: @@ -1666,40 +1664,27 @@ def get_dft_data(self, dft_chunk): mp._get_dft_data(dft_chunk, arr) return arr - # TODO: make add_near2far_uneven, wrap add_near2far around it, change parameters of _add_near2far - def add_near2far_uneven(self, freqs, nfreq, *near2fars, **kwargs): - freqs.sort() + def add_near2far(self, fcen, df, nfreq, *near2fars, **kwargs): nperiods = kwargs.get('nperiods', 1) - n2f = DftNear2Far(self._add_near2far, [freqs, nfreq, nperiods, near2fars]) + n2f = DftNear2Far(self._add_near2far, [fcen, df, nfreq, nperiods, near2fars]) self.dft_objects.append(n2f) return n2f - def add_near2far(self, fcen, df, nfreq, *near2fars, **kwargs): - freqs = np.linspace(fcen - df / 2, fcen + df / 2, nfreq) - return self.add_near2far_uneven(freqs, nfreq, *near2fars, **kwargs) - - def _add_near2far(self, freqs, nfreq, nperiods, near2fars): + def _add_near2far(self, fcen, df, nfreq, nperiods, near2fars): if self.fields is None: self.init_sim() - return self._add_fluxish_stuff(self.fields.add_dft_near2far, freqs, nfreq, near2fars, nperiods) + return self._add_fluxish_stuff(self.fields.add_dft_near2far, fcen, df, nfreq, near2fars, nperiods) - # TODO: make add_energy_uneven, wrap add_energy around it, change parameters of _add_energy - def add_energy_uneven(self, freqs, nfreq, *energys): - freqs.sort() - en = DftEnergy(self._add_energy, [freqs, nfreq, energys]) + def add_energy(self, fcen, df, nfreq, *energys): + en = DftEnergy(self._add_energy, [fcen, df, nfreq, energys]) self.dft_objects.append(en) return en - def add_energy(self, fcen, df, nfreq, *energys): - freqs = np.linspace(fcen - df / 2, fcen + df / 2, nfreq) - return self.add_energy_uneven(freqs, nfreq, *energys) - - def _add_energy(self, freqs, nfreq, energys): + def _add_energy(self, fcen, df, nfreq, energys): if self.fields is None: self.init_sim() - return self._add_fluxish_stuff(self.fields.add_dft_energy, freqs, nfreq, energys) + return self._add_fluxish_stuff(self.fields.add_dft_energy, fcen, df, nfreq, energys) - # TODO: check if this needs to be changed def _display_energy(self, name, func, energys): if energys: freqs = get_energy_freqs(energys[0]) @@ -1785,23 +1770,16 @@ def load_minus_near2far_data(self, n2f, n2fdata): self.load_near2far_data(n2f, n2fdata) n2f.scale_dfts(complex(-1.0)) - # TODO: make add_force_uneven, wrap add_force around it, change parameters of _add_force - def add_force_uneven(self, freqs, nfreq, *forces): - freqs.sort() - force = DftForce(self._add_force, [freqs, nfreq, forces]) + def add_force(self, fcen, df, nfreq, *forces): + force = DftForce(self._add_force, [fcen, df, nfreq, forces]) self.dft_objects.append(force) return force - def add_force(self, fcen, df, nfreq, *forces): - freqs = np.linspace(fcen - df / 2, fcen + df / 2, nfreq) - return self.add_force_uneven(freqs, nfreq, *forces) - - def _add_force(self, freqs, nfreq, forces): + def _add_force(self, fcen, df, nfreq, forces): if self.fields is None: self.init_sim() - return self._add_fluxish_stuff(self.fields.add_dft_force, freqs, nfreq, forces) + return self._add_fluxish_stuff(self.fields.add_dft_force, fcen, df, nfreq, forces) - # TODO: check if display_forces needs to be changed def display_forces(self, *forces): force_freqs = get_force_freqs(forces[0]) display_csv(self, 'force', zip(force_freqs, *[get_forces(f) for f in forces])) @@ -1834,34 +1812,22 @@ def load_minus_force_data(self, force, fdata): self.load_force_data(force, fdata) force.scale_dfts(complex(-1.0)) - # TODO: make add_flux_uneven, wrap add_flux around it, change parameters of _add_flux - def add_flux_uneven(self, freqs, nfreq, *fluxes): - freqs.sort() - flux = DftFlux(self._add_flux, [freqs, nfreq, fluxes]) + def add_flux(self, fcen, df, nfreq, *fluxes): + flux = DftFlux(self._add_flux, [fcen, df, nfreq, fluxes]) self.dft_objects.append(flux) return flux - def add_flux(self, fcen, df, nfreq, *fluxes): - freqs = np.linspace(fcen - df / 2, fcen + df / 2, nfreq) - return self.add_flux_uneven(freqs, nfreq, *fluxes) - - def _add_flux(self, freqs, nfreq, fluxes): + def _add_flux(self, fcen, df, nfreq, fluxes): if self.fields is None: self.init_sim() - return self._add_fluxish_stuff(self.fields.add_dft_flux, freqs, nfreq, fluxes) + return self._add_fluxish_stuff(self.fields.add_dft_flux, fcen, df, nfreq, fluxes) - # TODO: make add_mode_monitor_uneven, wrap add_mode_monitor around it, change parameters of _add_mode_monitor - def add_mode_monitor_uneven(self, freqs, nfreq, *fluxes): - freqs.sort() - flux = DftFlux(self._add_mode_monitor, [freqs, nfreq, fluxes]) + def add_mode_monitor(self, fcen, df, nfreq, *fluxes): + flux = DftFlux(self._add_mode_monitor, [fcen, df, nfreq, fluxes]) self.dft_objects.append(flux) return flux - def add_mode_monitor(self, fcen, df, nfreq, *fluxes): - freqs = np.linspace(fcen - df / 2, fcen + df / 2, nfreq) - return self.add_mode_monitor_uneven(freqs, nfreq, *fluxes) - - def _add_mode_monitor(self, freqs, nfreq, fluxes): + def _add_mode_monitor(self, fcen, df, nfreq, fluxes): if self.fields is None: self.init_sim() @@ -1873,7 +1839,7 @@ def _add_mode_monitor(self, freqs, nfreq, fluxes): d0 = region.direction d = self.fields.normal_direction(v.swigobj) if d0 < 0 else d0 - return self.fields.add_mode_monitor(d, v.swigobj, freqs, nfreq) + return self.fields.add_mode_monitor(d, v.swigobj, fcen - df / 2, fcen + df / 2, nfreq) def add_eigenmode(self, fcen, df, nfreq, *fluxes): warnings.warn('add_eigenmode is deprecated. Please use add_mode_monitor instead.', DeprecationWarning) @@ -1970,8 +1936,7 @@ def solve_cw(self, tol=1e-8, maxiters=10000, L=2): self._evaluate_dft_objects() return self.fields.solve_cw(tol, maxiters, L) - # TODO: change parameters of _add_fluxish_stuff - def _add_fluxish_stuff(self, add_dft_stuff, freqs, nfreq, stufflist, *args): + def _add_fluxish_stuff(self, add_dft_stuff, fcen, df, nfreq, stufflist, *args): vol_list = None for s in stufflist: @@ -1984,7 +1949,7 @@ def _add_fluxish_stuff(self, add_dft_stuff, freqs, nfreq, stufflist, *args): is_cylindrical=self.is_cylindrical).swigobj vol_list = mp.make_volume_list(v2, c, s.weight, vol_list) - stuff = add_dft_stuff(vol_list, freqs, nfreq, *args) + stuff = add_dft_stuff(vol_list, fcen - df / 2, fcen + df / 2, nfreq, *args) vol_list.__swig_destroy__(vol_list) return stuff @@ -2980,13 +2945,11 @@ def Ldos(fcen, df, nfreq): return mp._dft_ldos(fcen - df / 2, fcen + df / 2, nfreq) -# TODO: make dft_ldos_uneven, wrap dft_ldos wrap around it -def dft_ldos_uneven(freqs=None, nfreq=None, ldos=None): +def dft_ldos(fcen=None, df=None, nfreq=None, ldos=None): if ldos is None: - if freqs is None or nfreq is None: - raise ValueError("Either freqs and nfreq or an Ldos is required for dft_ldos_uneven") - freqs.sort() - ldos = mp._dft_ldos(freqs, nfreq) + if fcen is None or df is None or nfreq is None: + raise ValueError("Either fcen, df, and nfreq, or an Ldos is required for dft_ldos") + ldos = mp._dft_ldos(fcen - df / 2, fcen + df / 2, nfreq) def _ldos(sim, todo): if todo == 'step': @@ -2999,21 +2962,12 @@ def _ldos(sim, todo): return _ldos -def dft_ldos(fcen=None, df=None, nfreq=None, ldos=None): - if ldos is None: - if fcen is None or df is None or nfreq is None: - raise ValueError("Either fcen, df, and nfreq, or an Ldos is required for dft_ldos") - freqs = np.linspace(fcen - df / 2, fcen + df / 2, nfreq) - return dft_ldos_uneven(freqs, nfreq, ldos) - - def scale_flux_fields(s, flux): flux.scale_dfts(s) -# TODO: change this def get_flux_freqs(f): - return f.freqs + return np.linspace(f.freq_min, f.freq_min + f.dfreq * f.Nfreq, num=f.Nfreq, endpoint=False).tolist() def get_fluxes(f): @@ -3024,14 +2978,12 @@ def scale_force_fields(s, force): force.scale_dfts(s) -# TODO: change this def get_eigenmode_freqs(f): - return f.freqs + return np.linspace(f.freq_min, f.freq_min + f.dfreq * f.Nfreq, num=f.Nfreq, endpoint=False).tolist() -# TODO: change this def get_force_freqs(f): - return f.freqs + return np.linspace(f.freq_min, f.freq_min + f.dfreq * f.Nfreq, num=f.Nfreq, endpoint=False).tolist() def get_forces(f): @@ -3042,19 +2994,16 @@ def scale_near2far_fields(s, n2f): n2f.scale_dfts(s) -# TODO: change this def get_near2far_freqs(f): - return f.freqs + return np.linspace(f.freq_min, f.freq_min + f.dfreq * f.Nfreq, num=f.Nfreq, endpoint=False).tolist() def scale_energy_fields(s, ef): - # TODO: check if this is a typo - ef.scale_dfts(s) + df.scale_dfts(s) -# TODO: change this def get_energy_freqs(f): - return f.freqs + return np.linspace(f.freq_min, f.freq_min + f.dfreq * f.Nfreq, num=f.Nfreq, endpoint=False).tolist() def get_electric_energy(f): From 772a90625e766619b6fa9fdc05454bb007334fad Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Sun, 22 Mar 2020 18:32:01 -0600 Subject: [PATCH 29/35] Replace freqs arrays with std::vectors --- src/meep.hpp | 62 +++++++++++++++++++++++++++------------------------- 1 file changed, 32 insertions(+), 30 deletions(-) diff --git a/src/meep.hpp b/src/meep.hpp index 7de9b112c..2e03073df 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1011,7 +1011,7 @@ class dft_chunk { // the frequencies to loop_in_chunks // TODO: change dft_chunk properties int Nomega; - double *omegas; + std::vector omegas; component c; // component to DFT (possibly transformed by symmetry) @@ -1083,8 +1083,9 @@ void load_dft_hdf5(dft_chunk *dft_chunks, const char *name, h5file *file, const class dft_flux { public: // TODO: add dft_flux constructor - dft_flux(const component cE_, const component cH_, dft_chunk *E_, dft_chunk *H_, double *fs, - int Nf, const volume &where_, direction normal_direction_, bool use_symmetry_); + dft_flux(const component cE_, const component cH_, dft_chunk *E_, dft_chunk *H_, + std::vector fs, int Nf, const volume &where_, direction normal_direction_, + bool use_symmetry_); dft_flux(const component cE_, const component cH_, dft_chunk *E_, dft_chunk *H_, double fmin, double fmax, int Nf, const volume &where_, direction normal_direction_, bool use_symmetry_); @@ -1109,7 +1110,7 @@ class dft_flux { // TODO: change dft_flux properties int Nfreq; - double *freqs; + std::vector freqs; dft_chunk *E, *H; component cE, cH; volume where; @@ -1121,8 +1122,8 @@ class dft_flux { class dft_energy { public: // TODO: add dft_energy constructor - dft_energy(dft_chunk *E_, dft_chunk *H_, dft_chunk *D_, dft_chunk *B_, double *fs, int Nf, - const volume &where_); + dft_energy(dft_chunk *E_, dft_chunk *H_, dft_chunk *D_, dft_chunk *B_, std::vector fs, + int Nf, const volume &where_); dft_energy(dft_chunk *E_, dft_chunk *H_, dft_chunk *D_, dft_chunk *B_, double fmin, double fmax, int Nf, const volume &where_); dft_energy(const dft_energy &f); @@ -1150,7 +1151,7 @@ class dft_energy { // TODO: change dft_energy properties int Nfreq; - double *freqs; + std::vector freqs; dft_chunk *E, *H, *D, *B; volume where; }; @@ -1159,8 +1160,8 @@ class dft_energy { class dft_force { public: // TODO: add dft_force constructor - dft_force(dft_chunk *offdiag1_, dft_chunk *offdiag2_, dft_chunk *diag_, double *fs, int Nf, - const volume &where_); + dft_force(dft_chunk *offdiag1_, dft_chunk *offdiag2_, dft_chunk *diag_, std::vector fs, + int Nf, const volume &where_); dft_force(dft_chunk *offdiag1_, dft_chunk *offdiag2_, dft_chunk *diag_, double fmin, double fmax, int Nf, const volume &where_); dft_force(const dft_force &f); @@ -1181,7 +1182,7 @@ class dft_force { // TODO: change dft_force properties int Nfreq; - double *freqs; + std::vector freqs; dft_chunk *offdiag1, *offdiag2, *diag; volume where; }; @@ -1192,8 +1193,8 @@ class dft_near2far { /* fourier tranforms of tangential E and H field components in a medium with the given scalar eps and mu */ // TODO: add dft_near2far constructor - dft_near2far(dft_chunk *F, double *fs, int Nf, double eps, double mu, const volume &where_, - const direction periodic_d_[2], const int periodic_n_[2], + dft_near2far(dft_chunk *F, std::vector fs, int Nf, double eps, double mu, + const volume &where_, const direction periodic_d_[2], const int periodic_n_[2], const double periodic_k_[2], const double period_[2]); dft_near2far(dft_chunk *F, double fmin, double fmax, int Nf, double eps, double mu, const volume &where_, const direction periodic_d_[2], const int periodic_n_[2], @@ -1234,7 +1235,7 @@ class dft_near2far { // TODO: change dft_near2far properties int Nfreq; - double *freqs; + std::vector freqs; dft_chunk *F; double eps, mu; volume where; @@ -1253,12 +1254,12 @@ class dft_near2far { class dft_ldos { public: // TODO: add dft_lods constructor - dft_ldos(double *freqs, int Nfreq); + dft_ldos(std::vector freqs, int Nfreq); dft_ldos(double freq_min, double freq_max, int Nfreq); ~dft_ldos() { delete[] Fdft; delete[] Jdft; - delete[] omegas; + // delete[] omegas; } void update(fields &f); // to be called after each timestep @@ -1273,14 +1274,14 @@ class dft_ldos { public: // TODO: change dft_ldos properties int Nomega; - double *omegas; + std::vector omegas; }; // dft.cpp (normally created with fields::add_dft_fields) class dft_fields { public: // TODO: add dft_fields constructor - dft_fields(dft_chunk *chunks, double *freqs, int Nfreq, const volume &where); + dft_fields(dft_chunk *chunks, std::vector freqs, int Nfreq, const volume &where); dft_fields(dft_chunk *chunks, double freq_min, double freq_max, int Nfreq, const volume &where); void scale_dfts(std::complex scale); @@ -1289,7 +1290,7 @@ class dft_fields { // TODO: change dft_fields properties int Nfreq; - double *freqs; + std::vector freqs; dft_chunk *chunks; volume where; }; @@ -1748,7 +1749,7 @@ class fields { // dft.cpp // TODO: overload add_dft - dft_chunk *add_dft(component c, const volume &where, double *freqs, int Nfreq, + dft_chunk *add_dft(component c, const volume &where, std::vector freqs, int Nfreq, bool include_dV_and_interp_weights = true, std::complex stored_weight = 1.0, dft_chunk *chunk_next = 0, bool sqrt_dV_and_interp_weights = false, @@ -1761,40 +1762,41 @@ class fields { std::complex extra_weight = 1.0, bool use_centered_grid = true, int vc = 0); // TODO: overload add_dft_pt - dft_chunk *add_dft_pt(component c, const vec &where, double *freqs, int Nfreq); + dft_chunk *add_dft_pt(component c, const vec &where, std::vector freqs, int Nfreq); dft_chunk *add_dft_pt(component c, const vec &where, double freq_min, double freq_max, int Nfreq); // TODO: overload add_dft - dft_chunk *add_dft(const volume_list *where, double *freqs, int Nfreq, bool include_dV = true); + dft_chunk *add_dft(const volume_list *where, std::vector freqs, int Nfreq, + bool include_dV = true); dft_chunk *add_dft(const volume_list *where, double freq_min, double freq_max, int Nfreq, bool include_dV = true); void update_dfts(); // TODO: overload add_dft_flux - dft_flux add_dft_flux(const volume_list *where, double *freqs, int Nfreq, + dft_flux add_dft_flux(const volume_list *where, std::vector freqs, int Nfreq, bool use_symmetry = true); dft_flux add_dft_flux(const volume_list *where, double freq_min, double freq_max, int Nfreq, bool use_symmetry = true); // TODO: overload add_dft_flux - dft_flux add_dft_flux(direction d, const volume &where, double *freqs, int Nfreq, + dft_flux add_dft_flux(direction d, const volume &where, std::vector freqs, int Nfreq, bool use_symmetry = true); dft_flux add_dft_flux(direction d, const volume &where, double freq_min, double freq_max, int Nfreq, bool use_symmetry = true); // TODO: overload add_dft_flux_box - dft_flux add_dft_flux_box(const volume &where, double *freqs, int Nfreq); + dft_flux add_dft_flux_box(const volume &where, std::vector freqs, int Nfreq); dft_flux add_dft_flux_box(const volume &where, double freq_min, double freq_max, int Nfreq); // TODO: overload add_dft_flux_plane - dft_flux add_dft_flux_plane(const volume &where, double *freqs, int Nfreq); + dft_flux add_dft_flux_plane(const volume &where, std::vector freqs, int Nfreq); dft_flux add_dft_flux_plane(const volume &where, double freq_min, double freq_max, int Nfreq); // a "mode monitor" is just a dft_flux with symmetry reduction turned off. // TODO: overload add_mode_monitor - dft_flux add_mode_monitor(direction d, const volume &where, double *freqs, int Nfreq); + dft_flux add_mode_monitor(direction d, const volume &where, std::vector freqs, int Nfreq); dft_flux add_mode_monitor(direction d, const volume &where, double freq_min, double freq_max, int Nfreq); // TODO: overload add_dft_fields dft_fields add_dft_fields(component *components, int num_components, const volume where, - double *freqs, int Nfreq, bool use_centered_grid=true); + std::vector freqs, int Nfreq, bool use_centered_grid=true); dft_fields add_dft_fields(component *components, int num_components, const volume where, double freq_min, double freq_max, int Nfreq, bool use_centered_grid=true); @@ -1844,17 +1846,17 @@ class fields { std::complex overlaps[2]); // TODO: overload add_dft_energy - dft_energy add_dft_energy(const volume_list *where, double *freqs, int Nfreq); + dft_energy add_dft_energy(const volume_list *where, std::vector freqs, int Nfreq); dft_energy add_dft_energy(const volume_list *where, double freq_min, double freq_max, int Nfreq); // stress.cpp // TODO: overload add_dft_force - dft_force add_dft_force(const volume_list *where, double *freqs, int Nfreq); + dft_force add_dft_force(const volume_list *where, std::vector freqs, int Nfreq); dft_force add_dft_force(const volume_list *where, double freq_min, double freq_max, int Nfreq); // near2far.cpp // TODO: overload add_dft_near2far - dft_near2far add_dft_near2far(const volume_list *where, double *freqs, int Nfreq, + dft_near2far add_dft_near2far(const volume_list *where, std::vector freqs, int Nfreq, int Nperiods = 1); dft_near2far add_dft_near2far(const volume_list *where, double freq_min, double freq_max, int Nfreq, int Nperiods = 1); From 0e549747cf0fa84895b2524e51879d746b770d95 Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Sun, 22 Mar 2020 18:40:13 -0600 Subject: [PATCH 30/35] Replace freqs arrays with std::vectors --- src/meep.hpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/meep.hpp b/src/meep.hpp index 2e03073df..8bf883b1b 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1012,6 +1012,7 @@ class dft_chunk { // TODO: change dft_chunk properties int Nomega; std::vector omegas; + omegas.reserve(Nomega); component c; // component to DFT (possibly transformed by symmetry) @@ -1111,6 +1112,7 @@ class dft_flux { // TODO: change dft_flux properties int Nfreq; std::vector freqs; + freqs.reserve(Nfreq); dft_chunk *E, *H; component cE, cH; volume where; @@ -1152,6 +1154,7 @@ class dft_energy { // TODO: change dft_energy properties int Nfreq; std::vector freqs; + freqs.reserve(Nfreq); dft_chunk *E, *H, *D, *B; volume where; }; @@ -1183,6 +1186,7 @@ class dft_force { // TODO: change dft_force properties int Nfreq; std::vector freqs; + freqs.reserve(Nfreq); dft_chunk *offdiag1, *offdiag2, *diag; volume where; }; @@ -1236,6 +1240,7 @@ class dft_near2far { // TODO: change dft_near2far properties int Nfreq; std::vector freqs; + freqs.reserve(Nfreq); dft_chunk *F; double eps, mu; volume where; @@ -1275,6 +1280,7 @@ class dft_ldos { // TODO: change dft_ldos properties int Nomega; std::vector omegas; + omegas.reserve(Nomega); }; // dft.cpp (normally created with fields::add_dft_fields) @@ -1291,6 +1297,7 @@ class dft_fields { // TODO: change dft_fields properties int Nfreq; std::vector freqs; + freqs.reserve(Nfreq); dft_chunk *chunks; volume where; }; From b6c2317bb3c1c7a5ef1de89909d35bf8094f55ab Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Sun, 22 Mar 2020 19:37:50 -0600 Subject: [PATCH 31/35] Replace freqs arrays with std::vectors --- src/dft.cpp | 87 +++++++++++++++++++++++++++-------------------------- 1 file changed, 45 insertions(+), 42 deletions(-) diff --git a/src/dft.cpp b/src/dft.cpp index 6986858b3..1c5560cf0 100644 --- a/src/dft.cpp +++ b/src/dft.cpp @@ -35,7 +35,8 @@ struct dft_chunk_data { // for passing to field::loop_in_chunks as void* int vc; // TODO: change dft_chunk_data properties int Nomega; - double *omegas = new double[Nomega]; + std::vector omegas; + omegas.reserve(Nomega); complex stored_weight, extra_weight; double dt_factor; bool include_dV_and_interp_weights; @@ -87,7 +88,7 @@ dft_chunk::dft_chunk(fields_chunk *fc_, ivec is_, ivec ie_, vec s0_, vec s1_, ve // TODO: change dft_chunk properties Nomega = data->Nomega; - omegas = data->omegas = new double [Nomega]; + omegas.assign(data->omegas.begin(), data->omegas.end()); dft_phase = new complex[Nomega]; N = 1; @@ -107,7 +108,7 @@ dft_chunk::~dft_chunk() { delete[] dft; delete[] dft_phase; // TODO: add omegas to destructor? - delete[] omegas; + //delete[] omegas; // delete from fields_chunk list dft_chunk *cur = fc->dft_chunks; @@ -122,7 +123,7 @@ dft_chunk::~dft_chunk() { void dft_flux::remove() { // TODO: add freqs to remove? - delete[] freqs; + // freqs.clear(); while (E) { dft_chunk *nxt = E->next_in_dft; delete E; @@ -151,7 +152,7 @@ static void add_dft_chunkloop(fields_chunk *fc, int ichunk, component cgrid, ive } // TODO: overload fields::add_dft -dft_chunk *fields::add_dft(component c, const volume &where, double *freqs, int Nfreq, +dft_chunk *fields::add_dft(component c, const volume &where, std::vector freqs, int Nfreq, bool include_dV_and_interp_weights, complex stored_weight, dft_chunk *chunk_next, bool sqrt_dV_and_interp_weights, complex extra_weight, bool use_centered_grid, int vc) { @@ -172,9 +173,9 @@ dft_chunk *fields::add_dft(component c, const volume &where, double *freqs, int // TODO: change how freqs becomes omegas if (Nfreq < 1) abort("Nfreq must be at least 1"); data.Nomega = Nfreq; - double omegas[Nfreq]; + std::vector omegas; for (int i = 0; i < Nfreq; i++) { omegas[i] = freqs[i] * 2 * pi; } - data.omegas = omegas; + data.omegas.assign(omegas.begin(), omegas.end()); data.stored_weight = stored_weight; data.extra_weight = extra_weight; data.dt_factor = dt / sqrt(2.0 * pi); @@ -196,7 +197,7 @@ dft_chunk *fields::add_dft(component c, const volume &where, double freq_min, do bool use_centered_grid, int vc) { if (Nfreq < 1) abort("Nfreq must be at least 1"); double dfreq = (freq_max - freq_min) / (Nfreq - 1); - double freqs[Nfreq]; + std::vector freqs; for (int i = 0; i < Nfreq; i++) { freqs[i] = freq_min + dfreq * i; } @@ -205,7 +206,7 @@ dft_chunk *fields::add_dft(component c, const volume &where, double freq_min, do } // TODO: overload add_dft -dft_chunk *fields::add_dft(const volume_list *where, double *freqs, int Nfreq, +dft_chunk *fields::add_dft(const volume_list *where, std::vector freqs, int Nfreq, bool include_dV_and_interp_weights) { dft_chunk *chunks = 0; while (where) { @@ -233,7 +234,7 @@ dft_chunk *fields::add_dft(const volume_list *where, double freq_min, double fre } // TODO: overload add_dft_pt -dft_chunk *fields::add_dft_pt(component c, const vec &where, double *freqs, int Nfreq) { +dft_chunk *fields::add_dft_pt(component c, const vec &where, std::vector freqs, int Nfreq) { // TODO: change here return add_dft(c, where, freqs, Nfreq, false); } @@ -383,13 +384,13 @@ void load_dft_hdf5(dft_chunk *dft_chunks, component c, h5file *file, const char // TODO: overload dft_flux constructor dft_flux::dft_flux(const component cE_, const component cH_, dft_chunk *E_, dft_chunk *H_, - double *fs, int Nf, const volume &where_, direction normal_direction_, + std::vector fs, int Nf, const volume &where_, direction normal_direction_, bool use_symmetry_) : Nfreq(Nf), E(E_), H(H_), cE(cE_), cH(cH_), where(where_), normal_direction(normal_direction_), use_symmetry(use_symmetry_) { // TODO: change here if (Nf < 1) abort("Nf must be at least 1"); - freqs = fs; + freqs.assign(fs.begin(), fs.end()); } dft_flux::dft_flux(const component cE_, const component cH_, dft_chunk *E_, dft_chunk *H_, @@ -399,17 +400,17 @@ dft_flux::dft_flux(const component cE_, const component cH_, dft_chunk *E_, dft_ use_symmetry(use_symmetry_) { if (Nf < 1) abort("Nf must be at least 1"); double dfreqs = (fmax - fmin) / (Nf - 1); - double fs[Nf]; + std::vector fs; for (int i = 0; i < Nf; i++) { fs[i] = fmin + dfreqs * i; } - freqs = fs; + freqs.assign(fs.begin(), fs.end()); } dft_flux::dft_flux(const dft_flux &f) : where(f.where) { // TODO: change here Nfreq = f.Nfreq; - freqs = f.freqs; + freqs.assign(f.freqs.begin(), f.freqs.end()); E = f.E; H = f.H; cE = f.cE; @@ -463,7 +464,7 @@ void dft_flux::scale_dfts(complex scale) { } // TODO: overload fields::add_dft_flux -dft_flux fields::add_dft_flux(const volume_list *where_, double *freqs, int Nfreq, +dft_flux fields::add_dft_flux(const volume_list *where_, std::vector freqs, int Nfreq, bool use_symmetry) { if (!where_) // handle empty list of volumes // TODO: change here @@ -521,7 +522,7 @@ dft_flux fields::add_dft_flux(const volume_list *where_, double freq_min, double int Nfreq, bool use_symmetry) { if (Nfreq < 1) abort("Nfreq must be at least 1"); double dfreq = (freq_max - freq_min) / (Nfreq - 1); - double freqs[Nfreq]; + std::vector freqs; for (int i = 0; i < Nfreq; i++) { freqs[i] = freq_min + dfreq * i; } @@ -529,12 +530,12 @@ dft_flux fields::add_dft_flux(const volume_list *where_, double freq_min, double } // TODO: overload dft_energy constructor -dft_energy::dft_energy(dft_chunk *E_, dft_chunk *H_, dft_chunk *D_, dft_chunk *B_, double *fs, - int Nf, const volume &where_) +dft_energy::dft_energy(dft_chunk *E_, dft_chunk *H_, dft_chunk *D_, dft_chunk *B_, + std::vector fs, int Nf, const volume &where_) : Nfreq(Nf), E(E_), H(H_), D(D_), B(B_), where(where_) { // TODO: change here if (Nf < 1) abort("Nf must be at least 1"); - freqs = fs; + freqs.assign(fs.begin(), fs.end()); } dft_energy::dft_energy(dft_chunk *E_, dft_chunk *H_, dft_chunk *D_, dft_chunk *B_, double fmin, @@ -542,18 +543,18 @@ dft_energy::dft_energy(dft_chunk *E_, dft_chunk *H_, dft_chunk *D_, dft_chunk *B : Nfreq(Nf), E(E_), H(H_), D(D_), B(B_), where(where_) { if (Nf < 1) abort("Nf must be at least 1"); double dfreq = (fmax - fmin) / (Nf - 1); - double fs[Nf]; + std::vector fs; for (int i = 0; i < Nf; i++) { fs[i] = fmin + dfreq * i; } - freqs = fs; + freqs.assign(fs.begin(), fs.end()); } // TODO: change here dft_energy::dft_energy(const dft_energy &f) : where(f.where) { // TODO: change here Nfreq = f.Nfreq; - freqs = f.freqs; + freqs.assign(f.freqs.begin(), f.freqs.end()); E = f.E; H = f.H; D = f.D; @@ -602,7 +603,7 @@ double *dft_energy::total() { } // TODO: overload fields::add_dft_energy -dft_energy fields::add_dft_energy(const volume_list *where_, double *freqs, int Nfreq) { +dft_energy fields::add_dft_energy(const volume_list *where_, std::vector freqs, int Nfreq) { if (!where_) // handle empty list of volumes // TODO: change here @@ -635,7 +636,7 @@ dft_energy fields::add_dft_energy(const volume_list *where_, double freq_min, do int Nfreq) { if (Nfreq < 1) abort("Nfreq must be at least 1"); double dfreq = (freq_max - freq_min) / (Nfreq - 1); - double freqs[Nfreq]; + std::vector freqs; for (int i = 0; i < Nfreq; i++) { freqs[i] = freq_min + dfreq * i; } @@ -683,7 +684,7 @@ void dft_energy::scale_dfts(complex scale) { void dft_energy::remove() { // TODO: add freqs to remove? - delete[] freqs; + //delete[] freqs; while (E) { dft_chunk *nxt = E->next_in_dft; delete E; @@ -723,8 +724,8 @@ direction fields::normal_direction(const volume &where) const { } // TODO: overload fields::add_dft_flux -dft_flux fields::add_dft_flux(direction d, const volume &where, double *freqs, int Nfreq, - bool use_symmetry) { +dft_flux fields::add_dft_flux(direction d, const volume &where, std::vector freqs, + int Nfreq, bool use_symmetry) { if (d == NO_DIRECTION) d = normal_direction(where); volume_list vl(where, direction_component(Sx, d)); // TODO: change here @@ -737,7 +738,7 @@ dft_flux fields::add_dft_flux(direction d, const volume &where, double freq_min, int Nfreq, bool use_symmetry) { if (Nfreq < 1) abort("Nfreq must be at least 1"); double dfreq = (freq_max - freq_min) / (Nfreq - 1); - double freqs[Nfreq]; + std::vector freqs]; for (int i = 0; i < Nfreq; i++) { freqs[i] = freq_min + dfreq * i; } @@ -745,7 +746,8 @@ dft_flux fields::add_dft_flux(direction d, const volume &where, double freq_min, } // TODO: overload fields::add_mode_monitor -dft_flux fields::add_mode_monitor(direction d, const volume &where, double *freqs, int Nfreq) { +dft_flux fields::add_mode_monitor(direction d, const volume &where, std::vector freqs, + int Nfreq) { // TODO: change here return add_dft_flux(d, where, freqs, Nfreq, /*use_symmetry=*/false); } @@ -754,7 +756,7 @@ dft_flux fields::add_mode_monitor(direction d, const volume &where, double freq_ double freq_max, int Nfreq) { if (Nfreq < 1) abort("Nfreq must be at least 1"); double dfreq = (freq_max - freq_min) / (Nfreq - 1); - double freqs[Nfreq]; + std::vector freqs; for (int i = 0; i < Nfreq; i++) { freqs[i] = freq_min + dfreq * i; } @@ -762,7 +764,7 @@ dft_flux fields::add_mode_monitor(direction d, const volume &where, double freq_ } // TODO: overload fields::add_dft_flux_box -dft_flux fields::add_dft_flux_box(const volume &where, double *freqs, int Nfreq) { +dft_flux fields::add_dft_flux_box(const volume &where, std::vector freqs, int Nfreq) { volume_list *faces = 0; LOOP_OVER_DIRECTIONS(where.dim, d) { if (where.in_direction(d) > 0) { @@ -786,7 +788,7 @@ dft_flux fields::add_dft_flux_box(const volume &where, double freq_min, double f int Nfreq) { if (Nfreq < 1) abort("Nfreq must be at least 1"); double dfreq = (freq_max - freq_min) / (Nfreq - 1); - double freqs[Nfreq]; + std::vector freqs; for (int i = 0; i < Nfreq; i++) { freqs[i] = freq_min + dfreq * i; } @@ -794,7 +796,7 @@ dft_flux fields::add_dft_flux_box(const volume &where, double freq_min, double f } // TODO: overload fields::add_dft_flux_plane -dft_flux fields::add_dft_flux_plane(const volume &where, double *freqs, int Nfreq) { +dft_flux fields::add_dft_flux_plane(const volume &where, std::vector freqs, int Nfreq) { // TODO: change here return add_dft_flux(NO_DIRECTION, where, freqs, Nfreq); } @@ -803,7 +805,7 @@ dft_flux fields::add_dft_flux_plane(const volume &where, double freq_min, double int Nfreq) { if (Nfreq < 1) abort("Nfreq must be at least 1"); double dfreq = (freq_max - freq_min) / (Nfreq - 1); - double freqs[Nfreq]; + std::vector freqs; for (int i = 0; i < Nfreq; i++) { freqs[i] = freq_min + dfreq * i; } @@ -811,12 +813,13 @@ dft_flux fields::add_dft_flux_plane(const volume &where, double freq_min, double } // TODO: overload dft_fields constructor -dft_fields::dft_fields(dft_chunk *chunks_, double *freqs_, int Nfreq_, const volume &where_) +dft_fields::dft_fields(dft_chunk *chunks_, std::vector freqs_, int Nfreq_, + const volume &where_) : where(where_) { chunks = chunks_; // TODO: change here Nfreq = Nfreq_; - freqs = freqs_; + freqs.assgin(freqs_.begin(), freqs_.end()); } dft_fields::dft_fields(dft_chunk *chunks_, double freq_min_, double freq_max_, int Nfreq_, @@ -825,18 +828,18 @@ dft_fields::dft_fields(dft_chunk *chunks_, double freq_min_, double freq_max_, i chunks = chunks_; Nfreq = Nfreq_; double dfreqs = (freq_max_ - freq_min_) / (Nfreq_ - 1); - double freqs_[Nfreq_]; + std::vector freqs_; for (int i = 0; i < Nfreq_; i++) { freqs_[i] = freq_min_ + dfreqs * i; } - freqs = freqs_; + freqs.assign(freqs_.begin(), freqs_.end()); } void dft_fields::scale_dfts(cdouble scale) { chunks->scale_dft(scale); } void dft_fields::remove() { // TODO: add freqs to remove? - delete[] freqs; + // delete[] freqs; while (chunks) { dft_chunk *nxt = chunks->next_in_dft; delete chunks; @@ -846,7 +849,7 @@ void dft_fields::remove() { // TODO: overload fields::add_dft_fields dft_fields fields::add_dft_fields(component *components, int num_components, const volume where, - double *freqs, int Nfreq, bool use_centered_grid) { + std::vector freqs, int Nfreq, bool use_centered_grid) { bool include_dV_and_interp_weights = false; bool sqrt_dV_and_interp_weights = false; // default option from meep.hpp (expose to user?) std::complex extra_weight = 1.0; // default option from meep.hpp (expose to user?) @@ -866,7 +869,7 @@ dft_fields fields::add_dft_fields(component *components, int num_components, con double freq_min, double freq_max, int Nfreq, bool use_centered_grid) { double dfreq = (freq_max - freq_min) / (Nfreq - 1); - double freqs[Nfreq]; + std::vector freqs; for (int i = 0; i < Nfreq; i++) { freqs[i] = freq_min + dfreq * i; } From ab83da1373da7f24de9035d4d839179a36f68a9f Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Sun, 22 Mar 2020 20:06:44 -0600 Subject: [PATCH 32/35] Replace freqs arrays with std::vectors --- src/dft_ldos.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/dft_ldos.cpp b/src/dft_ldos.cpp index ec780f231..03827b6fd 100644 --- a/src/dft_ldos.cpp +++ b/src/dft_ldos.cpp @@ -23,13 +23,13 @@ using namespace std; namespace meep { // TODO: overload dft_ldos constructor -dft_ldos::dft_ldos(double *freqs, int Nfreq) { +dft_ldos::dft_ldos(std::vector freqs, int Nfreq) { if (Nfreq < 1) { abort("Nfreq must be at least 1"); } else { Nomega = Nfreq; - double om[Nomega]; + std::vector om; for (int i = 0; i < Nfreq; i++) { om[i] = freqs[i] * 2 * pi; } - omegas = om; + omegas.assign(om.begin(), om.end()); } Fdft = new complex[Nomega]; Jdft = new complex[Nomega]; @@ -44,9 +44,9 @@ dft_ldos::dft_ldos(double freq_min, double freq_max, int Nfreq) { else { double dfreq = (freq_max - freq_min) / (Nfreq - 1); Nomega = Nfreq; - double om[Nomega]; + std::vector om; for (int i = 0; i < Nfreq; i++) { om[i] = (freq_min + dfreq * i) * 2 * pi; } - omegas = om; + omegas.assign(om.begin(), om.end()); } Fdft = new complex[Nomega]; Jdft = new complex[Nomega]; From 12c40aac659e8f11ad9b9c271e60528ba6d207ba Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Sun, 22 Mar 2020 20:09:11 -0600 Subject: [PATCH 33/35] Replace freqs arrays with std::vectors --- src/mpb.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mpb.cpp b/src/mpb.cpp index 234298d6b..cdc73270d 100644 --- a/src/mpb.cpp +++ b/src/mpb.cpp @@ -751,7 +751,7 @@ void fields::add_eigenmode_source(component c0, const src_time &src, direction d /* similarly, if kpoints is non-null it should point to a */ /* caller-allocated array of size num_bands*num_freqs, which on*/ /* return will be populated by the k-vectors for the modes. */ -/***************************************************************/ +/****************omegas.assign(om.begin(), om.end());***********************************************/ void fields::get_eigenmode_coefficients(dft_flux flux, const volume &eig_vol, int *bands, int num_bands, int parity, double eig_resolution, double eigensolver_tol, std::complex *coeffs, @@ -759,7 +759,7 @@ void fields::get_eigenmode_coefficients(dft_flux flux, const volume &eig_vol, in void *user_kpoint_data, vec *kpoints, vec *kdom_list, double *cscale, direction d) { int num_freqs = flux.Nfreq; - double *freqs = flux.freqs; + std::vector freqs(flux.freqs); bool match_frequency = true; if (flux.use_symmetry && S.multiplicity() > 1 && parity == 0) From fe3c65004ac81e491be5f083f15b47f342587c0f Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Sun, 22 Mar 2020 20:17:15 -0600 Subject: [PATCH 34/35] Replace freqs arrays with std::vectors --- src/near2far.cpp | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/near2far.cpp b/src/near2far.cpp index a0e2583be..751a57255 100644 --- a/src/near2far.cpp +++ b/src/near2far.cpp @@ -30,14 +30,14 @@ using namespace std; namespace meep { // TODO: overload dft_near2far constructor -dft_near2far::dft_near2far(dft_chunk *F_, double *fs, int Nf, double eps_, double mu_, +dft_near2far::dft_near2far(dft_chunk *F_, std::vector fs, int Nf, double eps_, double mu_, const volume &where_, const direction periodic_d_[2], const int periodic_n_[2], const double periodic_k_[2], const double period_[2]) : Nfreq(Nf), F(F_), eps(eps_), mu(mu_), where(where_) { // TODO: initialize freqs with fs if Nf >= 1 if (Nf < 1) abort("Nf must be at least 1"); - freqs = fs; + freqs.assign(fs.begin(), fs.end()); for (int i = 0; i < 2; ++i) { periodic_d[i] = periodic_d_[i]; periodic_n[i] = periodic_n_[i]; @@ -54,11 +54,9 @@ dft_near2far::dft_near2far(dft_chunk *F_, double fmin, double fmax, int Nf, doub // TODO: create fs array with which to initalize freqs if (Nf < 1) abort("Nf must be at least 1"); double dfreq = (fmax - fmin) / (Nf - 1); - double fs[Nf]; for (int i = 0; i < Nf; i++) { - fs[i] = fmin + dfreq * i; + freqs[i] = fmin + dfreq * i; } - freqs = fs; for (int i = 0; i < 2; ++i) { periodic_d[i] = periodic_d_[i]; periodic_n[i] = periodic_n_[i]; @@ -69,7 +67,9 @@ dft_near2far::dft_near2far(dft_chunk *F_, double fmin, double fmax, int Nf, doub // TODO: change here for updated properties of dft_near2far dft_near2far::dft_near2far(const dft_near2far &f) - : freqs(f.freqs), Nfreq(f.Nfreq), F(f.F), eps(f.eps), mu(f.mu), where(f.where) { + : Nfreq(f.Nfreq), F(f.F), eps(f.eps), mu(f.mu), where(f.where) { + if (Nf < 1) abort("Nf must be at least 1"); + freqs.assign(f.freqs.begin(), f.freqs.end()); for (int i = 0; i < 2; ++i) { periodic_d[i] = f.periodic_d[i]; periodic_n[i] = f.periodic_n[i]; @@ -80,7 +80,7 @@ dft_near2far::dft_near2far(const dft_near2far &f) void dft_near2far::remove() { // TODO: add freqs to remove? - delete[] freqs; + // delete[] freqs; while (F) { dft_chunk *nxt = F->next_in_dft; delete F; @@ -549,8 +549,8 @@ double *dft_near2far::flux(direction df, const volume &where, double resolution) static double approxeq(double a, double b) { return fabs(a - b) < 0.5e-11 * (fabs(a) + fabs(b)); } // TODO: overload fields::add_dft_near2far -dft_near2far fields::add_dft_near2far(const volume_list *where, double *freqs, int Nfreq, - int Nperiods) { +dft_near2far fields::add_dft_near2far(const volume_list *where, std::vector freqs, + int Nfreq, int Nperiods) { dft_chunk *F = 0; /* E and H chunks*/ double eps = 0, mu = 0; volume everywhere = where->v; @@ -640,7 +640,7 @@ dft_near2far fields::add_dft_near2far(const volume_list *where, double freq_min, int Nfreq, int Nperiods) { if (Nfreq < 1) abort("Nfreq must be at least 1"); double dfreq = (freq_max - freq_min) / (Nfreq - 1); - double freqs[Nfreq]; + std::vector freqs; for (int i = 0; i < Nfreq; i++) { freqs[i] = freq_min + dfreq * i; } From 9fdc3c91e6e32d1fea580ee666463914e1467bdf Mon Sep 17 00:00:00 2001 From: danielwboyce Date: Sun, 22 Mar 2020 20:23:32 -0600 Subject: [PATCH 35/35] Replace freqs arrays with std::vectors --- src/stress.cpp | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/src/stress.cpp b/src/stress.cpp index 364232f3f..04562c4a8 100644 --- a/src/stress.cpp +++ b/src/stress.cpp @@ -24,11 +24,11 @@ using namespace std; namespace meep { // TODO: overload dft_force constructor -dft_force::dft_force(dft_chunk *offdiag1_, dft_chunk *offdiag2_, dft_chunk *diag_, double *fs, - int Nf, const volume &where_) +dft_force::dft_force(dft_chunk *offdiag1_, dft_chunk *offdiag2_, dft_chunk *diag_, + std::vector fs, int Nf, const volume &where_) : where(where_) { if (Nf < 1) abort("Nf must be at least 1"); - freqs = fs; + freqs.assign(fs.begin(), fs.end()); Nfreq = Nf; offdiag1 = offdiag1_; offdiag2 = offdiag2_; @@ -42,11 +42,9 @@ dft_force::dft_force(dft_chunk *offdiag1_, dft_chunk *offdiag2_, dft_chunk *diag if (Nf < 1) abort("Nf must be at least 1"); Nfreq = Nf; double dfreq = (fmax - fmin) / (Nf - 1); - double fs[Nf]; for (int i = 0; i < Nf; i++) { - fs[i] = fmin + dfreq * i; + freqs[i] = fmin + dfreq * i; } - freqs = fs; dfreq = Nf <= 1 ? 0.0 : (fmax - fmin) / (Nf - 1); offdiag1 = offdiag1_; offdiag2 = offdiag2_; @@ -57,7 +55,7 @@ dft_force::dft_force(dft_chunk *offdiag1_, dft_chunk *offdiag2_, dft_chunk *diag dft_force::dft_force(const dft_force &f) : where(f.where) { // TODO: change here Nfreq = f.Nfreq; - freqs = f.freqs; + freqs.assign(f.freqs.begin(), f.freqs.end()); offdiag1 = f.offdiag1; offdiag2 = f.offdiag2; diag = f.diag; @@ -66,7 +64,7 @@ dft_force::dft_force(const dft_force &f) : where(f.where) { void dft_force::remove() { // TODO: add freqs to remove? - delete[] freqs; + // delete[] freqs; while (offdiag1) { dft_chunk *nxt = offdiag1->next_in_dft; delete offdiag1; @@ -152,7 +150,7 @@ void dft_force::scale_dfts(complex scale) { force to be computed, so they should be vector components (such as Ex, Ey, ... or Sx, ...) rather than pseudovectors (like Hx, ...). */ // TODO: overload fields::add_dft_force -dft_force fields::add_dft_force(const volume_list *where_, double *freqs, int Nfreq) { +dft_force fields::add_dft_force(const volume_list *where_, std::vector freqs, int Nfreq) { dft_chunk *offdiag1 = 0, *offdiag2 = 0, *diag = 0; volume_list *where = S.reduce(where_); @@ -202,7 +200,7 @@ dft_force fields::add_dft_force(const volume_list *where_, double freq_min, doub int Nfreq) { if (Nfreq < 1) abort("Nfreq must be at least 1"); double dfreq = (freq_max - freq_min) / (Nfreq - 1); - double freqs[Nfreq]; + std::vector freqs; for (int i = 0; i < Nfreq; i++) { freqs[i] = freq_min + dfreq * i; }