From 80810a59a7f99cde1c480fd44298ec2d15546d55 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Thu, 8 Apr 2021 13:42:47 -0400 Subject: [PATCH] add needs_boundary_fix to add_srcdata --- src/fix_boundary_sources.cpp | 22 +++++++++++----------- src/meep.hpp | 2 +- src/meep_internals.hpp | 5 ++++- src/sources.cpp | 9 ++++----- 4 files changed, 20 insertions(+), 18 deletions(-) diff --git a/src/fix_boundary_sources.cpp b/src/fix_boundary_sources.cpp index 2046ba8a6..5de8d192b 100644 --- a/src/fix_boundary_sources.cpp +++ b/src/fix_boundary_sources.cpp @@ -45,28 +45,28 @@ void fields::fix_boundary_sources() { // they are actually supposed to be located, storing info in boundarysources. for (int i = 0; i < num_chunks; i++) { FOR_FIELD_TYPES(ft) { - for (src_vol *src = chunks[i]->sources[ft]; src; src = src->next) - if (src->needs_boundary_fix) { - for (size_t ipt = 0; ipt < src->npts; ++ipt) { - component c = src->c; - ivec here = chunks[i]->gv.iloc(c, src->index[ipt]); - if (!chunks[i]->gv.owns(here) && src->A[ipt] != 0.0) { - if (src->t->id == 0) abort("bug: fix_boundary_sources called for non-registered source"); + for (src_vol &src : chunks[i]->sources[ft]) + if (src.needs_boundary_fix) { + for (size_t ipt = 0; ipt < src.num_points(); ++ipt) { + component c = src.c; + ivec here = chunks[i]->gv.iloc(c, src.index_at(ipt)); + if (!chunks[i]->gv.owns(here) && src.amplitude(ipt) != 0.0) { + if (src.t()->id == 0) abort("bug: fix_boundary_sources called for non-registered source"); // find the chunk that owns this point, similar to logic in boundaries.cpp std::complex thephase; if (locate_component_point(&c, &here, &thephase) && !on_metal_boundary(here)) { for (int j = 0; j < num_chunks; j++) if (chunks[j]->gv.owns(here)) { - srcpt_info s = { src->A[ipt]*conj(thephase), chunks[j]->gv.index(c, here), src->t->id, chunks[j]->chunk_idx, c }; + srcpt_info s = { src.amplitude(ipt)*conj(thephase), chunks[j]->gv.index(c, here), src.t()->id, chunks[j]->chunk_idx, c }; boundarysources.push_back(s); break; } } - src->A[ipt] = 0.0; // will no longer be needed + src.set_amplitude(ipt, 0.0); // will no longer be needed } } - src->needs_boundary_fix = false; + src.needs_boundary_fix = false; } } } @@ -141,7 +141,7 @@ for (int psrc = 0; psrc < P; ++psrc) sourcedata srcdata = { (component) c, idx_arr, chunk_idx, amp_arr }; src_time *srctime = lookup_src_time(src_time_id); if (srctime == NULL) abort("bug: unknown src_time_id (missing registration?)"); - add_srcdata(srcdata, srctime); + add_srcdata(srcdata, srctime, size_t(0), NULL, false); if (idx < N) { chunk_idx = srcpts[idx].chunk_idx; src_time_id = srcpts[idx].src_time_id; diff --git a/src/meep.hpp b/src/meep.hpp index 7516ace10..fd2ccce3f 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1897,7 +1897,7 @@ class fields { void require_source_components(); void _require_component(component c, bool aniso2d); void require_component(component c) { _require_component(c, is_aniso2d()); sync_chunk_connections(); } - void add_srcdata(struct sourcedata cur_data, src_time *src, size_t n, std::complex* amp_arr); + void add_srcdata(struct sourcedata cur_data, src_time *src, size_t n, std::complex* amp_arr, bool needs_boundary_fix=false); void register_src_time(src_time *src); src_time *lookup_src_time(size_t id); diff --git a/src/meep_internals.hpp b/src/meep_internals.hpp index 566e217e8..c3461b3b1 100644 --- a/src/meep_internals.hpp +++ b/src/meep_internals.hpp @@ -55,13 +55,16 @@ class src_vol { // 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); + std::vector > &&s, bool fix_boundaries=false); // 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]; } + std::complex amplitude(size_t j) const { return amp[j]; }; + void set_amplitude(size_t j, std::complex a) { amp[j] = a; }; + void set_amplitude(size_t j, double a) { amp[j] = a; }; 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; }; diff --git a/src/sources.cpp b/src/sources.cpp index 0c7b9ac5a..cc1dc91d3 100644 --- a/src/sources.cpp +++ b/src/sources.cpp @@ -166,13 +166,13 @@ bool custom_src_time::is_equal(const src_time &t) const { /*********************************************************************/ src_vol::src_vol(component cc, src_time *st, std::vector &&ind, - std::vector > &&s) + std::vector > &&s, bool fix_boundaries) : 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)), needs_boundary_fix(false) { + src_t(st), index(std::move(ind)), amp(std::move(amps)), needs_boundary_fix(fix_boundaries) { assert(index.size() == amp.size()); } @@ -289,7 +289,7 @@ static void src_vol_chunkloop(fields_chunk *fc, int ichunk, component c, ivec is 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){ +void fields::add_srcdata(struct sourcedata cur_data, src_time *src, size_t n, std::complex* amp_arr, bool needs_boundary_fix){ if (n == 0) { n = cur_data.idx_arr.size(); assert(amp_arr == NULL); @@ -300,13 +300,12 @@ void fields::add_srcdata(struct sourcedata cur_data, src_time *src, size_t n, st std::vector index_arr(cur_data.idx_arr); std::vector> amplitudes(amp_arr, amp_arr+n); component c = cur_data.near_fd_comp; - 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->add_source(ft, src_vol(c, src, std::move(index_arr), std::move(amplitudes))); + fc->add_source(ft, src_vol(c, src, std::move(index_arr), std::move(amplitudes), needs_boundary_fix)); // 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()