diff --git a/python/meep.i b/python/meep.i index 80ef75024..de0166b2a 100644 --- a/python/meep.i +++ b/python/meep.i @@ -1470,6 +1470,7 @@ void _get_gradient(PyObject *grad, PyObject *fields_a, PyObject *fields_f, PyObj %ignore is_medium; %ignore is_metal; %ignore meep::choose_chunkdivision; +%ignore meep::fields_chunk; %ignore meep::infinity; %ignore std::vector::vector(size_type); diff --git a/scheme/Makefile.am b/scheme/Makefile.am index dfdfc1464..9b8329b54 100644 --- a/scheme/Makefile.am +++ b/scheme/Makefile.am @@ -39,7 +39,7 @@ endif # workaround missing namespace prefix in swig meep_renames.i: $(LIBHDRS) - (echo "// AUTOMATICALLY GENERATED -- DO NOT EDIT"; perl -pe 's/^ *class +([A-Za-z_0-9:]*)( *| *:[^{]*){.*$$/%rename(meep_\1) meep::\1;/' $(LIBHDRS) | grep "%rename" | sort -u; echo; grep -hv typedef $(LIBHDRS) | perl -pe 's/(inline|const|extern|static) +//g' | perl -pe 's/^[A-Za-z_0-9:<>]+[* ]+([A-Za-z_0-9:]*) *\(.*$$/%rename(meep_\1) meep::\1;/' | grep "%rename" | sed '/choose_chunkdivision/d' | sort -u; ) > $@ + (echo "// AUTOMATICALLY GENERATED -- DO NOT EDIT"; perl -pe 's/^ *class +([A-Za-z_0-9:]*)( *| *:[^{]*){.*$$/%rename(meep_\1) meep::\1;/' $(LIBHDRS) | grep "%rename" | sed '/meep_fields_chunk/d' | sort -u; echo; grep -hv typedef $(LIBHDRS) | perl -pe 's/(inline|const|extern|static) +//g' | perl -pe 's/^[A-Za-z_0-9:<>]+[* ]+([A-Za-z_0-9:]*) *\(.*$$/%rename(meep_\1) meep::\1;/' | grep "%rename" | sed '/choose_chunkdivision/d' | sort -u; ) > $@ # work around bug in swig, where it doesn't prepend namespace to friend funcs meep_swig_bug_workaround.i: $(LIBHDRS) diff --git a/scheme/meep.i b/scheme/meep.i index 6526da1d6..a03fde7f6 100644 --- a/scheme/meep.i +++ b/scheme/meep.i @@ -281,6 +281,7 @@ static meep::vec my_kpoint_func(double freq, int mode, void *user_data) { %warnfilter(302,325,451,503,509); %ignore meep::choose_chunkdivision; +%ignore meep::fields_chunk; %ignore meep_geom::fragment_stats; %include "meep_renames.i" diff --git a/src/array_slice.cpp b/src/array_slice.cpp index 37a8487a9..e39484d37 100644 --- a/src/array_slice.cpp +++ b/src/array_slice.cpp @@ -209,18 +209,18 @@ static void get_source_slice_chunkloop(fields_chunk *fc, int ichnk, component cg if (has_direction(dim, Y)) NY = ((slice_imax - slice_imin).in_direction(Y) / 2) + 1; for (int ft = 0; ft < NUM_FIELD_TYPES; ft++) - for (src_vol *s = fc->sources[ft]; s; s = s->next) { + for (const src_vol &s : fc->get_sources(static_cast(ft))) { component cS = S.transform(data->source_component, -sn); - if (s->c != cS) continue; + if (s.c != cS) continue; // loop over point sources in this src_vol. for each point source, // the src_vol stores the amplitude and the global index of the // symmetry-parent grid point, from which we need to compute the // local index of the symmetry-child grid point within this // slice (that is, if it even lies within the slice) - for (size_t npt = 0; npt < s->npts; npt++) { - complex amp = s->A[npt]; - ptrdiff_t chunk_index = s->index[npt]; + for (size_t npt = 0; npt < s.num_points(); npt++) { + const complex & = s.amplitude_at(npt); + ptrdiff_t chunk_index = s.index_at(npt); ivec iloc_parent = fc->gv.iloc(Dielectric, chunk_index); ivec iloc_child = S.transform(iloc_parent, sn) + shift; if (!in_subgrid(slice_imin, iloc_child, slice_imax)) continue; // source point outside slice diff --git a/src/dft_ldos.cpp b/src/dft_ldos.cpp index a5348f4e7..b892efa36 100644 --- a/src/dft_ldos.cpp +++ b/src/dft_ldos.cpp @@ -102,41 +102,41 @@ void dft_ldos::update(fields &f) { for (int ic = 0; ic < f.num_chunks; ic++) if (f.chunks[ic]->is_mine()) { - for (src_vol *sv = f.chunks[ic]->sources[D_stuff]; sv; sv = sv->next) { - component c = direction_component(Ex, component_direction(sv->c)); + for (const src_vol &sv : f.chunks[ic]->get_sources(D_stuff)) { + component c = direction_component(Ex, component_direction(sv.c)); realnum *fr = f.chunks[ic]->f[c][0]; realnum *fi = f.chunks[ic]->f[c][1]; if (fr && fi) // complex E - for (size_t j = 0; j < sv->npts; j++) { - const ptrdiff_t idx = sv->index[j]; - const complex A = sv->A[j]; + for (size_t j = 0; j < sv.num_points(); j++) { + const ptrdiff_t idx = sv.index_at(j); + const complex& A = sv.amplitude_at(j); EJ += complex(fr[idx], fi[idx]) * conj(A); Jsum += abs(A); } else if (fr) { // E is purely real - for (size_t j = 0; j < sv->npts; j++) { - const ptrdiff_t idx = sv->index[j]; - const complex A = sv->A[j]; + for (size_t j = 0; j < sv.num_points(); j++) { + const ptrdiff_t idx = sv.index_at(j); + const complex& A = sv.amplitude_at(j); EJ += double(fr[idx]) * conj(A); Jsum += abs(A); } } } - for (src_vol *sv = f.chunks[ic]->sources[B_stuff]; sv; sv = sv->next) { - component c = direction_component(Hx, component_direction(sv->c)); + for (const src_vol& sv : f.chunks[ic]->get_sources(B_stuff)) { + component c = direction_component(Hx, component_direction(sv.c)); realnum *fr = f.chunks[ic]->f[c][0]; realnum *fi = f.chunks[ic]->f[c][1]; if (fr && fi) // complex H - for (size_t j = 0; j < sv->npts; j++) { - const ptrdiff_t idx = sv->index[j]; - const complex A = sv->A[j]; + for (size_t j = 0; j < sv.num_points(); j++) { + const ptrdiff_t idx = sv.index_at(j); + const complex& A = sv.amplitude_at(j); HJ += complex(fr[idx], fi[idx]) * conj(A); Jsum += abs(A); } else if (fr) { // H is purely real - for (size_t j = 0; j < sv->npts; j++) { - const ptrdiff_t idx = sv->index[j]; - const complex A = sv->A[j]; + for (size_t j = 0; j < sv.num_points(); j++) { + const ptrdiff_t idx = sv.index_at(j); + const complex& A = sv.amplitude_at(j); HJ += double(fr[idx]) * conj(A); Jsum += abs(A); } diff --git a/src/fields.cpp b/src/fields.cpp index ae6b63ec5..d48c2c9d0 100644 --- a/src/fields.cpp +++ b/src/fields.cpp @@ -15,6 +15,8 @@ % Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ +#include +#include #include #include #include @@ -197,7 +199,6 @@ fields_chunk::~fields_chunk() { dft_chunks = nxt; } FOR_FIELD_TYPES(ft) { - delete sources[ft]; delete[] zeroes[ft]; } FOR_FIELD_TYPES(ft) { @@ -246,7 +247,6 @@ fields_chunk::fields_chunk(structure_chunk *the_s, const char *od, double m, dou } doing_solve_cw = false; solve_cw_omega = 0.0; - FOR_FIELD_TYPES(ft) { sources[ft] = NULL; } FOR_COMPONENTS(c) DOCMP2 { f[c][cmp] = NULL; f_u[c][cmp] = NULL; @@ -308,7 +308,6 @@ fields_chunk::fields_chunk(const fields_chunk &thef, int chunkidx) : gv(thef.gv) } doing_solve_cw = thef.doing_solve_cw; solve_cw_omega = thef.solve_cw_omega; - FOR_FIELD_TYPES(ft) { sources[ft] = NULL; } FOR_COMPONENTS(c) DOCMP2 { f[c][cmp] = NULL; f_u[c][cmp] = NULL; @@ -484,8 +483,9 @@ void fields::require_source_components() { memset(needed, 0, sizeof(needed)); for (int i = 0; i < num_chunks; i++) { FOR_FIELD_TYPES(ft) { - for (src_vol *src = chunks[i]->sources[ft]; src; src = src->next) - needed[src->c] = 1; + for (const auto& src : chunks[i]->get_sources(ft)) { + needed[src.c] = 1; + } } } int allneeded[NUM_FIELD_COMPONENTS]; @@ -544,10 +544,22 @@ void fields::_require_component(component c, bool aniso2d) { } } +void fields_chunk::add_source(field_type ft, src_vol &&src) { + auto it = std::find_if(sources[ft].begin(), sources[ft].end(), [&src](const src_vol &other) { + return src_vol::combinable(src, other); + }); + + if (it != sources[ft].end()) { + it->add_amplitudes_from(src); + return; + } + + sources[ft].push_back(std::move(src)); +} + void fields_chunk::remove_sources() { FOR_FIELD_TYPES(ft) { - delete sources[ft]; - sources[ft] = NULL; + sources[ft].clear(); } } diff --git a/src/meep.hpp b/src/meep.hpp index 25bdeb575..93ae1c5b0 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1370,7 +1370,7 @@ class fields_chunk { bool zero_fields_near_cylorigin; // fields=0 m pixels near r=0 for stability double beta; int is_real; - src_vol *sources[NUM_FIELD_TYPES]; + std::vector sources[NUM_FIELD_TYPES]; structure_chunk *new_s; structure_chunk *s; const char *outdir; @@ -1398,6 +1398,10 @@ class fields_chunk { std::complex get_chi1inv(component, direction, const ivec &iloc, double frequency = 0) const; + // Returns the vector of sources volumes for field type `ft`. + const std::vector &get_sources(field_type ft) const { return sources[ft]; } + // Adds a source volume of field type `ft` and takes ownership of `src`. + void add_source(field_type ft, src_vol &&src); void backup_component(component c); void average_with_backup(component c); diff --git a/src/meep_internals.hpp b/src/meep_internals.hpp index 45941d45c..051b5e0d3 100644 --- a/src/meep_internals.hpp +++ b/src/meep_internals.hpp @@ -15,6 +15,7 @@ % Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ +#include #include #include "meep.hpp" @@ -53,35 +54,37 @@ inline int rmin_bulk(int m) { return r; } +// A source volume +// Moveable and copyable class src_vol { public: - src_vol(component cc, src_time *st, size_t n, ptrdiff_t *ind, std::complex *amps); - src_vol(const src_vol &sv); - ~src_vol() { - delete next; - delete[] index; - delete[] A; - } - - src_time *t; - ptrdiff_t *index; // list of locations of sources in grid (indices) - size_t npts; // number of points in list - component c; // field component the source applies to - std::complex *A; // list of amplitudes - - std::complex dipole(size_t j) { return A[j] * t->dipole(); } - std::complex current(size_t j) { return A[j] * t->current(); } - void update(double time, double dt) { t->update(time, dt); } - - bool operator==(const src_vol &sv) const { - // note: don't compare sv.A, since this is used to see if we can just - // add one source's amplitudes to another in src_vol::add_to - return sv.npts == npts && sv.c == c && sv.t == t && - memcmp(sv.index, index, npts * sizeof(ptrdiff_t)) == 0; - } - - src_vol *add_to(src_vol *others); - src_vol *next; + // Constructs a new source volume. Takes ownership of `ind` and `amps`. + // Requirement: ind.size() == amps.size() + src_vol(component cc, src_time *st, std::vector &&ind, + std::vector > &&s); + + // Checks whether `a` and `b` are combinable, i.e. have the same indices and point to the same + // `src_time` instance, but have potentially different amplitudes. + static bool combinable(const src_vol &a, const src_vol &b); + + ptrdiff_t index_at(size_t pos) const { return index[pos]; } + const std::complex &litude_at(size_t pos) const { return amp[pos]; } + size_t num_points() const { return index.size(); }; + const src_time *t() const { return src_t; }; + + std::complex dipole(size_t j) const { return amp[j] * src_t->dipole(); }; + std::complex current(size_t j) const { return amp[j] * src_t->current(); }; + void update(double time, double dt) { src_t->update(time, dt); }; + // Merges the amplitudes from volume source `other` into `this`. + // Requirement: other.num_points() == this->num_points(). + // It is recommended to use `combinable` before calling this method. + void add_amplitudes_from(const src_vol &other); + + const component c; // field component the source applies to +private: + src_time *src_t; // Not owned by us. + std::vector index; // locations of sources in grid (indices) + std::vector > amp; // amplitudes }; const int num_bandpts = 32; diff --git a/src/sources.cpp b/src/sources.cpp index b149b4192..4308377b0 100644 --- a/src/sources.cpp +++ b/src/sources.cpp @@ -15,6 +15,8 @@ % Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ +#include +#include #include #include #include @@ -152,52 +154,25 @@ bool custom_src_time::is_equal(const src_time &t) const { /*********************************************************************/ -src_vol::src_vol(component cc, src_time *st, size_t n, ptrdiff_t *ind, complex *amps) { - c = cc; - if (is_D(c)) c = direction_component(Ex, component_direction(c)); - if (is_B(c)) c = direction_component(Hx, component_direction(c)); - t = st; - next = NULL; - npts = n; - index = ind; - A = amps; +src_vol::src_vol(component cc, src_time *st, std::vector &&ind, + std::vector > &&s) + : c([](component c) -> component { + if (is_D(c)) c = direction_component(Ex, component_direction(c)); + if (is_B(c)) c = direction_component(Hx, component_direction(c)); + return c; + }(cc)), + src_t(st), index(std::move(ind)), amp(std::move(amps)) { + assert(index.size() == amp.size()); } -src_vol::src_vol(const src_vol &sv) { - c = sv.c; - t = sv.t; - npts = sv.npts; - index = new ptrdiff_t[npts]; - A = new complex[npts]; - for (size_t j = 0; j < npts; j++) { - index[j] = sv.index[j]; - A[j] = sv.A[j]; - } - if (sv.next) - next = new src_vol(*sv.next); - else - next = NULL; +bool src_vol::combinable(const src_vol &a, const src_vol &b) { + return (a.c == b.c) && (a.src_t == b.src_t) && (a.index == b.index); } -src_vol *src_vol::add_to(src_vol *others) { - if (others) { - if (*this == *others) { - if (npts != others->npts) - meep::abort("Cannot add grid_volume sources with different number of points\n"); - /* Compare all of the indices...if this ever becomes too slow, - we can just compare the first and last indices. */ - for (size_t j = 0; j < npts; j++) { - if (others->index[j] != index[j]) meep::abort("Different indices\n"); - others->A[j] += A[j]; - } - } - else - others->next = add_to(others->next); - return others; - } - else { - next = others; - return this; +void src_vol::add_amplitudes_from(const src_vol &other) { + assert(amp.size() == other.num_points()); + for (size_t i = 0; i < amp.size(); ++i) { + amp[i] += other.amp[i]; } } @@ -267,8 +242,8 @@ static void src_vol_chunkloop(fields_chunk *fc, int ichunk, component c, ivec is size_t npts = 1; LOOP_OVER_DIRECTIONS(is.dim, d) { npts *= (ie.in_direction(d) - is.in_direction(d)) / 2 + 1; } - ptrdiff_t *index_array = new ptrdiff_t[npts]; - complex *amps_array = new complex[npts]; + std::vector index_array(npts); + std::vector> amps_array(npts); complex amp = data->amp * conj(shift_phase); @@ -299,26 +274,22 @@ static void src_vol_chunkloop(fields_chunk *fc, int ichunk, component c, ivec is if (idx_vol > npts) meep::abort("add_volume_source: computed wrong npts (%zd vs. %zd)", npts, idx_vol); - src_vol *tmp = new src_vol(c, data->src, idx_vol, index_array, amps_array); field_type ft = is_magnetic(c) ? B_stuff : D_stuff; - fc->sources[ft] = tmp->add_to(fc->sources[ft]); + fc->add_source(ft, src_vol(c, data->src, std::move(index_array), std::move(amps_array))); } void fields::add_srcdata(struct sourcedata cur_data, src_time *src, size_t n, std::complex* amp_arr){ sources = src->add_to(sources, &src); - ptrdiff_t* index_arr = new ptrdiff_t[n]; - std::complex* amp_arr_copy = new std::complex[n]; - for (size_t i = 0; i < n; i++){ - index_arr[i] = cur_data.idx_arr[i]; - amp_arr_copy[i] = amp_arr[i]; - } + std::vector index_arr(cur_data.idx_arr); + std::vector> amplitudes(amp_arr, amp_arr+n); component c = cur_data.near_fd_comp; - src_vol *tmp = new src_vol(c, src, n, index_arr, amp_arr_copy); + field_type ft = is_magnetic(c) ? B_stuff : D_stuff; if (0 > cur_data.fc_idx or cur_data.fc_idx >= num_chunks) meep::abort("fields chunk index out of range"); fields_chunk *fc = chunks[cur_data.fc_idx]; if (!fc->is_mine()) meep::abort("wrong fields chunk"); - fc->sources[ft] = tmp->add_to(fc->sources[ft]); + + fc->add_source(ft, src_vol(c, src, std::move(index_arr), std::move(amplitudes))); // We can't do require_component(c) since that only works if all processes are adding // srcdata for the same components in the same order, which may not be true. // ... instead, the caller should call fields::require_source_components() diff --git a/src/step.cpp b/src/step.cpp index 55a4980a2..9d452d20b 100644 --- a/src/step.cpp +++ b/src/step.cpp @@ -201,22 +201,22 @@ void fields::step_source(field_type ft, bool including_integrated) { } void fields_chunk::step_source(field_type ft, bool including_integrated) { if (doing_solve_cw && !including_integrated) return; - for (src_vol *sv = sources[ft]; sv; sv = sv->next) { - component c = direction_component(first_field_component(ft), component_direction(sv->c)); - const realnum *cndinv = s->condinv[c][component_direction(sv->c)]; - if ((including_integrated || !sv->t->is_integrated) && f[c][0] && - ((ft == D_stuff && is_electric(sv->c)) || (ft == B_stuff && is_magnetic(sv->c)))) { + for (const src_vol &sv : sources[ft]) { + component c = direction_component(first_field_component(ft), component_direction(sv.c)); + const realnum *cndinv = s->condinv[c][component_direction(sv.c)]; + if ((including_integrated || !sv.t()->is_integrated) && f[c][0] && + ((ft == D_stuff && is_electric(sv.c)) || (ft == B_stuff && is_magnetic(sv.c)))) { if (cndinv) - for (size_t j = 0; j < sv->npts; j++) { - const ptrdiff_t i = sv->index[j]; - const complex A = sv->current(j) * dt * double(cndinv[i]); + for (size_t j = 0; j < sv.num_points(); j++) { + const ptrdiff_t i = sv.index_at(j); + const complex A = sv.current(j) * dt * double(cndinv[i]); f[c][0][i] -= real(A); if (!is_real) f[c][1][i] -= imag(A); } else - for (size_t j = 0; j < sv->npts; j++) { - const complex A = sv->current(j) * dt; - const ptrdiff_t i = sv->index[j]; + for (size_t j = 0; j < sv.num_points(); j++) { + const complex A = sv.current(j) * dt; + const ptrdiff_t i = sv.index_at(j); f[c][0][i] -= real(A); if (!is_real) f[c][1][i] -= imag(A); } diff --git a/src/update_eh.cpp b/src/update_eh.cpp index 3f619eb0a..858e5f267 100644 --- a/src/update_eh.cpp +++ b/src/update_eh.cpp @@ -47,11 +47,12 @@ bool fields_chunk::update_eh(field_type ft, bool skip_w_components) { bool have_int_sources = false; if (!doing_solve_cw) { - for (src_vol *sv = sources[ft2]; sv; sv = sv->next) - if (sv->t->is_integrated) { + for (const src_vol &sv : sources[ft2]) { + if (sv.t()->is_integrated) { have_int_sources = true; break; } + } } FOR_FT_COMPONENTS(ft, ec) { @@ -103,12 +104,12 @@ bool fields_chunk::update_eh(field_type ft, bool skip_w_components) { // Next, subtract time-integrated sources (i.e. polarizations, not currents) if (have_f_minus_p && !doing_solve_cw) { - for (src_vol *sv = sources[ft2]; sv; sv = sv->next) { - if (sv->t->is_integrated && f[sv->c][0] && ft == type(sv->c)) { - component c = field_type_component(ft2, sv->c); - for (size_t j = 0; j < sv->npts; ++j) { - const complex A = sv->dipole(j); - DOCMP { f_minus_p[c][cmp][sv->index[j]] -= (cmp) ? imag(A) : real(A); } + for (const src_vol &sv : sources[ft2]) { + if (sv.t()->is_integrated && f[sv.c][0] && ft == type(sv.c)) { + component c = field_type_component(ft2, sv.c); + for (size_t j = 0; j < sv.num_points(); ++j) { + const complex A = sv.dipole(j); + DOCMP { f_minus_p[c][cmp][sv.index_at(j)] -= (cmp) ? imag(A) : real(A); } } } }