diff --git a/doc/docs/Python_User_Interface.md b/doc/docs/Python_User_Interface.md
index 6edc44d03..d9413ecd2 100644
--- a/doc/docs/Python_User_Interface.md
+++ b/doc/docs/Python_User_Interface.md
@@ -1380,6 +1380,7 @@ The most common step function is an output function, which outputs some field co
Note that although the various field components are stored at different places in the [Yee lattice](Yee_Lattice.md), when they are outputted they are all linearly interpolated to the same grid: to the points at the *centers* of the Yee cells, i.e. $(i+0.5,j+0.5,k+0.5)\cdotΔ$ in 3d.
+
**`output_epsilon()`**
—
Output the dielectric function (relative permittivity) ε. Note that this only outputs the frequency-independent part of ε (the $\omega\to\infty$ limit).
@@ -1462,6 +1463,26 @@ Returns the Fourier-transformed fields as a NumPy array.
+ `num_freq`: The index of the frequency: (an integer in the range `0...nfreq-1`, where `nfreq` is the number of frequencies stored in `dft_obj,` as set by the `nfreq` parameter to `add_dft_fields`, `add_dft_flux`, etc.)
+#### Source slices
+
+**`get_source_slice(component, vol=None, center=None, size=None)`**
+
+This routine returns an array of the same dimensions as that
+returned by `get_array` for the given `vol` or `center/size`,
+but containing information on the [*sources*](#source) at grid points,
+not the fields there.
+(Because sources, unlike fields, are *inputs* rather
+than outputs of meep calculations, this routine
+is not a means of extracting results from meep calculations,
+but is rather a tool to sanity-check visualization
+and confirm that the source distribution you specified is the one you
+wanted; philosophically it is closer in
+spirit to [`output_epsilon()`](#output_epsilon) than to
+ `get_array/get_dft_array`.
+
+The array returned by `get_source_slice` is always complex,
+and corresponds to the given `component`.
+
#### Harminv
The following step function collects field data from a given point and runs [Harminv](https://github.com/stevengj/harminv) on that data to extract the frequencies, decay rates, and other information.
diff --git a/libpympb/pympb.cpp b/libpympb/pympb.cpp
index 4f83726bf..01af180da 100644
--- a/libpympb/pympb.cpp
+++ b/libpympb/pympb.cpp
@@ -2376,9 +2376,7 @@ double mode_solver::compute_energy_in_objects(geometric_object_list objects) {
return 0.0;
}
- for (int i = 0; i < objects.num_items; ++i) {
- geom_fix_object(objects.items[i]);
- }
+ geom_fix_objects0(objects);
int n1 = mdata->nx;
int n2 = mdata->ny;
diff --git a/python/meep.i b/python/meep.i
index 7005f092b..234f98c69 100644
--- a/python/meep.i
+++ b/python/meep.i
@@ -69,6 +69,7 @@ using namespace meep_geom;
extern boolean point_in_objectp(vector3 p, GEOMETRIC_OBJECT o);
extern boolean point_in_periodic_objectp(vector3 p, GEOMETRIC_OBJECT o);
void display_geometric_object_info(int indentby, GEOMETRIC_OBJECT o);
+
%}
%include "numpy.i"
diff --git a/python/simulation.py b/python/simulation.py
index 476adbbc2..88e474cbc 100644
--- a/python/simulation.py
+++ b/python/simulation.py
@@ -1233,11 +1233,11 @@ def add_source(self, src):
if not fname.endswith('.h5'):
fname += '.h5'
- add_vol_src(fname, dset, src.amplitude * 1.0,)
+ add_vol_src(fname, dset, src.amplitude * 1.0)
elif src.amp_func:
add_vol_src(src.amp_func, src.amplitude * 1.0)
elif src.amp_data is not None:
- add_vol_src(src.amp_data, src.amplitude * 1.0,)
+ add_vol_src(src.amp_data, src.amplitude * 1.0)
else:
add_vol_src(src.amplitude * 1.0)
@@ -1594,6 +1594,18 @@ def get_dft_array(self, dft_obj, component, num_freq):
else:
raise ValueError("Invalid type of dft object: {}".format(dft_swigobj))
+ def get_source_slice(self, component, vol=None, center=None, size=None):
+ if vol is None and center is None and size is None:
+ v = self.fields.total_volume()
+ else:
+ v = self._volume_from_kwargs(vol, center, size)
+ dim_sizes = np.zeros(3, dtype=np.uintp)
+ self.fields.get_array_slice_dimensions(v, dim_sizes)
+ dims = [s for s in dim_sizes if s != 0]
+ arr = np.zeros(dims, dtype=np.complex128)
+ self.fields.get_source_slice(v, component ,arr)
+ return arr
+
def get_eigenmode_coefficients(self, flux, bands, eig_parity=mp.NO_PARITY, eig_vol=None,
eig_resolution=0, eig_tolerance=1e-12, kpoint_func=None, verbose=False):
if self.fields is None:
diff --git a/python/source.py b/python/source.py
index 94b3fb031..bf5d5fe27 100644
--- a/python/source.py
+++ b/python/source.py
@@ -10,7 +10,6 @@ def check_positive(prop, val):
else:
raise ValueError("{} must be positive. Got {}".format(prop, val))
-
class Source(object):
def __init__(self, src, component, center, size=Vector3(), amplitude=1.0, amp_func=None,
diff --git a/python/tests/simulation.py b/python/tests/simulation.py
index cecd6382f..e79f1b1a1 100644
--- a/python/tests/simulation.py
+++ b/python/tests/simulation.py
@@ -574,6 +574,11 @@ def print_field(sim):
self.assertAlmostEqual(result[0], -0.0599602798684155)
+ def test_source_slice(self):
+ sim = self.init_simple_simulation()
+ sim.run(until=5)
+ slice = sim.get_source_slice(mp.Ez)
+ print(slice)
if __name__ == '__main__':
unittest.main()
diff --git a/python/vec.i b/python/vec.i
index 09ad405d4..c43e1f875 100644
--- a/python/vec.i
+++ b/python/vec.i
@@ -57,6 +57,8 @@
%rename(vec_abs) meep::abs;
%rename(vec_max) meep::max;
%rename(vec_min) meep::min;
+%rename(vec_str) meep::vec::str;
+%rename(ivec_str) meep::ivec::str;
%rename(print_grid_volume) meep::grid_volume::print;
%rename(symmetry_reduce) meep::symmetry::reduce;
diff --git a/src/array_slice.cpp b/src/array_slice.cpp
index ea689f08a..6e508a413 100644
--- a/src/array_slice.cpp
+++ b/src/array_slice.cpp
@@ -53,6 +53,8 @@ typedef struct {
field_rfunction rfun;
void *fun_data;
std::vector components;
+ component source_slice_component;
+ bool get_source_slice;
void *vslice;
@@ -111,22 +113,80 @@ static void get_array_slice_dimensions_chunkloop(fields_chunk *fc, int ichnk, co
data->min_corner = min(data->min_corner, min(isS, ieS));
data->max_corner = max(data->max_corner, max(isS, ieS));
data->num_chunks++;
+}
+
+/*****************************************************************/
+/* populate the array slice with information about sources with */
+/* the component specified in source_slice_component. */
+/*****************************************************************/
+void fill_chunk_source_slice(fields_chunk *fc, array_slice_data *data)
+{
+ ndim dim=fc->gv.dim;
+ if (dim==Dcyl)
+ { fprintf(stderr,"warning: source slices not implemented for cylindrical coordinates; array will be all zeros\n");
+ return;
+ }
+ component slice_component = data->source_slice_component;
+ cdouble *slice = (cdouble *)data->vslice;
+ ivec slice_imin=data->min_corner, slice_imax=data->max_corner;
+
+ ptrdiff_t NY=1, NZ=1;
+ if (has_direction(fc->gv.dim,Z))
+ NZ = ((slice_imax-slice_imin).in_direction(Z)/2) + 1;
+ if (has_direction(fc->gv.dim,Y))
+ NY = ((slice_imax-slice_imin).in_direction(Y)/2) + 1;
+
+ for(int ft=0; ftsources[ft]; s; s=s->next)
+ {
+ if (slice_component!=s->c)
+ 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 grid point,
+ // from which we need to compute the local index of the grid point within the
+ // slice (that is, if it even lies within the slice)
+ for(size_t npt=0; nptnpts; npt++)
+ { cdouble amp = s->A[npt];
+ ptrdiff_t chunk_index = s->index[npt];
+ ivec iloc = fc->gv.iloc(Dielectric, chunk_index);
+ if (ilocslice_imax) continue; // source point outside slice
+ ivec slice_offset = iloc-slice_imin;
+ ptrdiff_t slice_index=0;
+ if (has_direction(dim,Z))
+ slice_index += slice_offset.in_direction(Z)/2;
+ if (has_direction(dim,Y))
+ slice_index += NZ*slice_offset.in_direction(Y)/2;
+ if (has_direction(dim,X))
+ slice_index += NY*NZ*slice_offset.in_direction(X)/2;
+
+ slice[slice_index] = amp;
+ }
+ }
}
/***************************************************************/
/* callback function passed to loop_in_chunks to fill array slice */
/***************************************************************/
static void get_array_slice_chunkloop(fields_chunk *fc, int ichnk, component cgrid,
- ivec is, ivec ie, vec s0, vec s1, vec e0, vec e1,
- double dV0, double dV1,
- ivec shift, complex shift_phase,
- const symmetry &S, int sn, void *data_)
+ ivec is, ivec ie, vec s0, vec s1, vec e0, vec e1,
+ double dV0, double dV1,
+ ivec shift, complex shift_phase,
+ const symmetry &S, int sn, void *data_)
{
UNUSED(ichnk);UNUSED(cgrid);UNUSED(s0);UNUSED(s1);UNUSED(e0);UNUSED(e1);
UNUSED(dV0);UNUSED(dV1);
array_slice_data *data = (array_slice_data *) data_;
+ //-----------------------------------------------------------------------//
+ //- If we're fetching a 'source slice,' branch off here to handle that. //
+ //-----------------------------------------------------------------------//
+ if (data->get_source_slice)
+ { fill_chunk_source_slice(fc,data);
+ return;
+ }
+
//-----------------------------------------------------------------------//
// Find output chunk dimensions and strides, etc.
@@ -146,11 +206,13 @@ static void get_array_slice_chunkloop(fields_chunk *fc, int ichnk, component cgr
permute.set_direction(d, abs(permute.in_direction(d)));
// compute the size of the chunk to output, and its strides etc.
+ size_t slice_size=1;
for (int i = 0; i < data->rank; ++i) {
direction d = data->ds[i];
int isd = isS.in_direction(d), ied = ieS.in_direction(d);
start[i] = (min(isd, ied) - data->min_corner.in_direction(d)) / 2;
count[i] = abs(ied - isd) / 2 + 1;
+ slice_size *= count[i];
if (ied < isd) offset[permute.in_direction(d)] = count[i] - 1;
}
@@ -162,7 +224,7 @@ static void get_array_slice_chunkloop(fields_chunk *fc, int ichnk, component cgr
direction d = data->ds[i];
dims[i]= (data->max_corner.in_direction(d)
- data->min_corner.in_direction(d)) / 2 + 1;
- };
+ }
ptrdiff_t stride[3]={1,1,1};
for (int i = 0; i < data->rank; ++i) {
@@ -171,13 +233,23 @@ static void get_array_slice_chunkloop(fields_chunk *fc, int ichnk, component cgr
for (int k = i + 1; k < data->rank; ++k) stride[j] *= dims[k];
offset[j] *= stride[j];
if (offset[j]) stride[j] *= -1;
- };
+ }
// sco="slice chunk offset"
ptrdiff_t sco=start[0]*dims[1]*dims[2] + start[1]*dims[2] + start[2];
//-----------------------------------------------------------------------//
- // Compute the function to output, exactly as in fields::integrate.
+ // Otherwise proceed to compute the function of field components to be //
+ // tabulated on the slice, exactly as in fields::integrate. //
+ //-----------------------------------------------------------------------//
+ double *slice=0;
+ cdouble *zslice=0;
+ bool complex_data = (data->rfun==0);
+ if (complex_data)
+ zslice = (cdouble *)data->vslice;
+ else
+ slice = (double *)data->vslice;
+
ptrdiff_t *off = data->offsets;
component *cS = data->cS;
complex *fields = data->fields, *ph = data->ph;
@@ -189,14 +261,6 @@ static void get_array_slice_chunkloop(fields_chunk *fc, int ichnk, component cgr
ptrdiff_t imos[6];
int num_components=data->components.size();
- double *slice=0;
- cdouble *zslice=0;
- bool complex_data = (data->rfun==0);
- if (complex_data)
- zslice = (cdouble *)data->vslice;
- else
- slice = (double *)data->vslice;
-
for (int i = 0; i < num_components; ++i) {
cS[i] = S.transform(data->components[i], -sn);
if (cS[i] == Dielectric || cS[i] == Permeability)
@@ -212,10 +276,16 @@ static void get_array_slice_chunkloop(fields_chunk *fc, int ichnk, component cgr
fc->gv.yee2cent_offsets(imcs[k], imos[2*k], imos[2*k+1]);
vec rshift(shift * (0.5*fc->gv.inva));
+ // main loop over all grid points owned by this field chunk.
LOOP_OVER_IVECS(fc->gv, is, ie, idx) {
+
+ // get real-space coordinates of grid point, taking into
+ // account the complications of symmetries.
IVEC_LOOP_LOC(fc->gv, loc);
loc = S.transform(loc, sn) + rshift;
+ // interpolate fields at the four nearest grid points
+ // to get the value of the field component for this point
for (int i = 0; i < num_components; ++i) {
if (cS[i] == Dielectric) {
double tr = 0.0;
@@ -250,22 +320,24 @@ static void get_array_slice_chunkloop(fields_chunk *fc, int ichnk, component cgr
fields[i] = complex(f[0], f[1]) * ph[i];
}
}
- int idx2 = sco + ((((offset[0] + offset[1] + offset[2])
- + loop_i1 * stride[0])
- + loop_i2 * stride[1])
- + loop_i3 * stride[2]);
+
+ // compute the index into the array for this grid point and store the result of the computation
+ ptrdiff_t idx2 = sco + ((((offset[0] + offset[1] + offset[2])
+ + loop_i1 * stride[0])
+ + loop_i2 * stride[1])
+ + loop_i3 * stride[2]);
if (complex_data)
zslice[idx2] = data->fun(fields, loc, data->fun_data);
else
slice[idx2] = data->rfun(fields, loc, data->fun_data);
- };
+ } // LOOP_OVER_IVECS
}
/***************************************************************/
-/* given a volume, fill in the dims[] and directions[] arrays */
+/* given a volume, fill in the dims[] and dirs[] arrays */
/* describing the array slice needed to store field data for */
/* all grid points in the volume. */
/* */
@@ -276,7 +348,8 @@ static void get_array_slice_chunkloop(fields_chunk *fc, int ichnk, component cgr
/* initialized appopriately for subsequent use in */
/* get_array_slice. */
/***************************************************************/
-int fields::get_array_slice_dimensions(const volume &where, size_t dims[3], void *caller_data)
+int fields::get_array_slice_dimensions(const volume &where, size_t dims[3],
+ void *caller_data)
{
am_now_working_on(FieldOutput);
@@ -305,6 +378,7 @@ int fields::get_array_slice_dimensions(const volume &where, size_t dims[3], void
if (rank >= 3) abort("too many dimensions in array_slice");
size_t n = (data->max_corner.in_direction(d)
- data->min_corner.in_direction(d)) / 2 + 1;
+
if (n > 1) {
data->ds[rank] = d;
dims[rank++] = n;
@@ -318,22 +392,25 @@ int fields::get_array_slice_dimensions(const volume &where, size_t dims[3], void
return rank;
}
-/***************************************************************/
-/* precisely one of fun, rfun should be non-NULL */
-/***************************************************************/
+/**********************************************************************/
+/* precisely one of fun, rfun, source_slice should be non-NULL / true */
+/**********************************************************************/
void *fields::do_get_array_slice(const volume &where,
std::vector components,
field_function fun,
field_rfunction rfun,
void *fun_data,
- void *vslice) {
-
+ void *vslice,
+ component source_slice_component,
+ bool get_source_slice) {
am_now_working_on(FieldOutput);
/***************************************************************/
/* call get_array_slice_dimensions to get slice dimensions and */
/* partially initialze an array_slice_data struct */
/***************************************************************/
+ // by tradition, empty dimensions in time-domain field arrays are *not* collapsed;
+ // TODO make this a caller-specifiable parameter to get_array_slice()?
size_t dims[3];
array_slice_data data;
int rank=get_array_slice_dimensions(where, dims, &data);
@@ -346,29 +423,29 @@ void *fields::do_get_array_slice(const volume &where,
if (vslice==0)
{ if (complex_data)
{ zslice = new cdouble[slice_size];
+ memset(zslice,0,slice_size*sizeof(cdouble));
vslice = (void *)zslice;
}
else
{ slice = new double[slice_size];
+ memset(slice,0,slice_size*sizeof(double));
vslice = (void *)slice;
- };
- };
-
- data.vslice = vslice;
- data.fun = fun;
- data.rfun = rfun;
- data.fun_data = fun_data;
- data.components = components;
-
- int num_components = components.size();
-
- data.cS = new component[num_components];
- data.ph = new cdouble[num_components];
- data.fields = new cdouble[num_components];
+ }
+ }
- data.offsets = new ptrdiff_t[2 * num_components];
- for (int i = 0; i < 2 * num_components; ++i)
- data.offsets[i] = 0;
+ data.vslice = vslice;
+ data.fun = fun;
+ data.rfun = rfun;
+ data.fun_data = fun_data;
+ data.source_slice_component = source_slice_component;
+ data.get_source_slice = get_source_slice;
+ data.components = components;
+ int num_components = components.size();
+ data.cS = new component[num_components];
+ data.ph = new cdouble[num_components];
+ data.fields = new cdouble[num_components];
+ data.offsets = new ptrdiff_t[2 * num_components];
+ memset(data.offsets, 0, 2*num_components*sizeof(ptrdiff_t));
/* compute inverse-epsilon directions for computing Dielectric fields */
data.ninveps = 0;
@@ -416,7 +493,7 @@ void *fields::do_get_array_slice(const volume &where,
memcpy(slice+offset, buffer, size*sizeof(cdouble));
remaining-=size;
offset+=size;
- };
+ }
delete[] buffer;
}
else
@@ -430,9 +507,9 @@ void *fields::do_get_array_slice(const volume &where,
memcpy(slice+offset, buffer, size*sizeof(double));
remaining-=size;
offset+=size;
- };
+ }
delete[] buffer;
- };
+ }
delete[] data.offsets;
delete[] data.fields;
@@ -499,4 +576,10 @@ cdouble *fields::get_complex_array_slice(const volume &where, component c,
(void *)slice);
}
+
+cdouble *fields::get_source_slice(const volume &where, component source_slice_component, cdouble *slice) {
+ vector cs; // empty
+ return (cdouble *)do_get_array_slice(where,cs,0,0,0,(void *)slice,source_slice_component, true);
+}
+
} // namespace meep
diff --git a/src/energy_and_flux.cpp b/src/energy_and_flux.cpp
index 2444c4a94..a4dcaa602 100644
--- a/src/energy_and_flux.cpp
+++ b/src/energy_and_flux.cpp
@@ -198,7 +198,7 @@ double fields::flux_in_box_wrongH(direction d, const volume &where) {
if (coordinate_mismatch(gv.dim, d))
return 0.0;
- component cE[2], cH[2];
+ component cE[2]={Ey, Ez}, cH[2]={Hz, Hy};
switch (d) {
case X: cE[0] = Ey, cE[1] = Ez, cH[0] = Hz, cH[1] = Hy; break;
case Y: cE[0] = Ez, cE[1] = Ex, cH[0] = Hx, cH[1] = Hz; break;
diff --git a/src/meep.hpp b/src/meep.hpp
index dd87986e4..5c3fdab6c 100644
--- a/src/meep.hpp
+++ b/src/meep.hpp
@@ -1304,6 +1304,7 @@ class chunkloop_field_components {
/***************************************************************/
typedef vec (*kpoint_func)(double freq, int mode, void *user_data);
+
class fields {
public:
int num_chunks;
@@ -1441,13 +1442,18 @@ class fields {
component c,
std::complex *slice=0);
+ // like get_array_slice, but for *sources* instead of fields
+ std::complex *get_source_slice(const volume &where, component source_slice_component, std::complex *slice=0);
+
// master routine for all above entry points
void *do_get_array_slice(const volume &where,
std::vector components,
field_function fun,
field_rfunction rfun,
void *fun_data,
- void *vslice);
+ void *vslice,
+ component source_slice_component=Ex,
+ bool get_source_slice=false);
// step.cpp methods:
double last_step_output_wall_time;
@@ -1471,6 +1477,7 @@ class fields {
int is_continuous = 0);
void add_point_source(component c, const src_time &src,
const vec &, std::complex amp = 1.0);
+
void add_volume_source(component c, const src_time &src, const volume &where_,
std::complex *arr, size_t dim1, size_t dim2, size_t dim3,
std::complex amp);
@@ -1482,8 +1489,7 @@ class fields {
std::complex A(const vec &),
std::complex amp = 1.0);
void add_volume_source(component c, const src_time &src,
- const volume &,
- std::complex amp = 1.0);
+ const volume &, std::complex amp = 1.0);
void require_component(component c);
// mpb.cpp
diff --git a/src/meep/vec.hpp b/src/meep/vec.hpp
index d82d0cf45..3c439658a 100644
--- a/src/meep/vec.hpp
+++ b/src/meep/vec.hpp
@@ -403,12 +403,13 @@ vec zero_vec(ndim);
class vec {
public:
- vec() {};
- vec(ndim di) { dim = di; };
+ vec() { init_t(); };
+ vec(ndim di) { init_t(); dim = di; };
vec(ndim di, double val) { dim = di; t[0]=t[1]=t[2]=t[3]=t[4]=val; };
- vec(double zz) { dim = D1; t[Z] = zz; };
- vec(double xx, double yy) { dim = D2; t[X] = xx; t[Y] = yy; };
+ vec(double zz) { init_t(); dim = D1; t[Z] = zz; };
+ vec(double xx, double yy) { init_t(); dim = D2; t[X] = xx; t[Y] = yy; };
vec(double xx, double yy, double zz) {
+ init_t();
dim = D3; t[X] = xx; t[Y] = yy; t[Z] = zz; };
friend vec veccyl(double rr, double zz);
~vec() {};
@@ -484,11 +485,20 @@ class vec {
double in_direction(direction d) const { return t[d]; };
void set_direction(direction d, double val) { t[d] = val; };
+ // pretty-print to a user-supplied buffer (if provided) or to a static internal buffer (in which case not thread-safe)
+ const char *str(char *buffer=0, size_t buflen=0);
+
double project_to_boundary(direction, double boundary_loc);
friend vec zero_vec(ndim);
friend vec one_vec(ndim);
private:
double t[5];
+ void init_t() {
+ for (int i = 0; i < 5; ++i) {
+ t[i] = 0;
+ }
+ }
+
};
inline double abs(const vec &pt) { return sqrt(pt & pt); }
@@ -526,12 +536,13 @@ ivec one_ivec(ndim);
class ivec {
public:
- ivec() { dim = D2; t[X]=t[Y]=t[Z]=0; };
- ivec(ndim di) { dim = di; t[X]=t[Y]=t[Z]=0; };
+ ivec() { init_t(); dim = D2; };
+ ivec(ndim di) { init_t(); dim = di; };
ivec(ndim di, int val) { dim = di; t[0]=t[1]=t[2]=t[3]=t[4]=val; };
- ivec(int zz) { dim = D1; t[X]=t[Y]=0; t[Z] = zz; };
- ivec(int xx, int yy) { dim = D2; t[X] = xx; t[Y] = yy; t[Z]=0; };
+ ivec(int zz) { init_t(); dim = D1; t[Z] = zz; };
+ ivec(int xx, int yy) { init_t(); dim = D2; t[X] = xx; t[Y] = yy; };
ivec(int xx, int yy, int zz) {
+ init_t();
dim = D3; t[X] = xx; t[Y] = yy; t[Z] = zz; };
friend ivec iveccyl(int xx, int yy);
~ivec() {};
@@ -618,6 +629,9 @@ class ivec {
int in_direction(direction d) const { return t[d]; };
void set_direction(direction d, int val) { t[d] = val; };
+ // pretty-print to a user-supplied buffer (if provided) or to a static internal buffer (in which case not thread-safe)
+ const char *str(char *buffer=0, size_t buflen=0);
+
ivec round_up_to_even(void) const {
ivec result(dim);
LOOP_OVER_DIRECTIONS(dim, d)
@@ -629,6 +643,11 @@ class ivec {
friend ivec one_ivec(ndim);
private:
int t[5];
+ void init_t() {
+ for (int i = 0; i < 5; ++i) {
+ t[i] = 0;
+ }
+ }
};
inline ivec zero_ivec(ndim di) {
@@ -728,6 +747,9 @@ class grid_volume {
public:
grid_volume() {};
+ grid_volume subvolume(ivec is, ivec ie);
+ void init_subvolume(ivec is, ivec ie);
+
ndim dim;
double a, inva /* = 1/a */;
diff --git a/src/meepgeom.cpp b/src/meepgeom.cpp
index d39f398fc..68c777bb0 100644
--- a/src/meepgeom.cpp
+++ b/src/meepgeom.cpp
@@ -656,7 +656,7 @@ static bool get_front_object(const meep::volume &v,
const geometric_object *o1 = 0, *o2 = 0;
vector3 shiftby1 = {0,0,0}, shiftby2 = {0,0,0};
geom_box pixel;
- material_type mat1, mat2;
+ material_type mat1 = vacuum, mat2 = vacuum;
int id1 = -1, id2 = -1;
const int num_neighbors[3] = { 3, 5, 9 };
const int neighbors[3][9][3] = {
@@ -2147,6 +2147,6 @@ dft_data::dft_data(int freqs, int components, std::vector volumes)
num_freqs(freqs),
num_components(components),
vols(volumes) {
-
}
+
} // namespace meep_geom
diff --git a/src/monitor.cpp b/src/monitor.cpp
index 80855f658..072f0e5f6 100644
--- a/src/monitor.cpp
+++ b/src/monitor.cpp
@@ -165,12 +165,9 @@ complex fields_chunk::get_field(component c, const ivec &iloc) const {
largest box in which you can interpolate the fields without communication.
It is *not* necessarily non-overlapping with other chunks. */
volume fields_chunk::get_field_gv(component c) const {
- switch (c) {
- case Dielectric: case Permeability:
- c = gv.eps_component();
- default:
- return volume(gv.loc(c, 0), gv.loc(c, gv.ntot() - 1));
- }
+ if (c==Dielectric || c==Permeability)
+ c = gv.eps_component();
+ return volume(gv.loc(c, 0), gv.loc(c, gv.ntot() - 1));
}
/* Non-collective, zero-communication get_field... loc *must*
diff --git a/src/multilevel-atom.cpp b/src/multilevel-atom.cpp
index 91bb33307..a05b4e13d 100644
--- a/src/multilevel-atom.cpp
+++ b/src/multilevel-atom.cpp
@@ -264,7 +264,7 @@ void multilevel_susceptibility::update_P
}
// compute E*8 at point i
- double E8[3][2];
+ double E8[3][2]={{0.0,0.0},{0.0,0.0},{0.0,0.0}};
for (idot = 0; idot < 3 && cdot[idot] != Dielectric; ++idot) {
realnum *w = W[cdot[idot]][0], *wp = W_prev[cdot[idot]][0];
E8[idot][0] = w[i]+w[i+o1[idot]]+w[i+o2[idot]]+w[i+o1[idot]+o2[idot]]
diff --git a/src/sources.cpp b/src/sources.cpp
index 6e130e306..7bc991f19 100644
--- a/src/sources.cpp
+++ b/src/sources.cpp
@@ -292,9 +292,10 @@ static void src_vol_chunkloop(fields_chunk *fc, int ichunk, component c,
if (!fc->gv.owns(iloc)) continue;
IVEC_LOOP_LOC(fc->gv, loc);
- loc += shift * (0.5*inva) - data->center;
+ loc += shift * (0.5*inva);
- amps_array[idx_vol] = IVEC_LOOP_WEIGHT(s0,s1,e0,e1,1) * amp * data->A(loc);
+ vec rel_loc = loc - data->center;
+ amps_array[idx_vol] = IVEC_LOOP_WEIGHT(s0,s1,e0,e1,1) * amp * data->A(rel_loc);
/* for "D" sources, multiply by epsilon. FIXME: this is not quite
right because it doesn't handle non-diagonal chi1inv!
diff --git a/src/vec.cpp b/src/vec.cpp
index 05dd0e179..3b22cfaa0 100644
--- a/src/vec.cpp
+++ b/src/vec.cpp
@@ -1471,5 +1471,62 @@ field_rfunction derived_component_func(derived_component c, const grid_volume &g
}
/***************************************************************************/
+/* utility methods for pretty-printing. may be called with no arguments, */
+/* in which case static internal buffers are used; NUMBUFS defines the */
+/* number of str() calls with no arguments that may be appear */
+/* simultaneously as e.g. arguments to a single invocation of printf(). */
+/***************************************************************************/
+#define BUFLEN 100
+#define NUMBUFS 10
+const char *ivec::str(char *buffer, size_t buflen)
+{ static char bufring[NUMBUFS][BUFLEN];
+ static int nbuf=0;
+ if (buffer==0)
+ { buffer=bufring[nbuf]; buflen=BUFLEN; nbuf=(nbuf+1)%NUMBUFS; }
+ if (dim==Dcyl)
+ snprintf(buffer,buflen,"{%i,%i}",t[R],t[Z]);
+ else
+ snprintf(buffer,buflen,"{%i,%i,%i}",t[X],t[Y],t[Z]);
+ return buffer;
+}
+
+const char *vec::str(char *buffer, size_t buflen)
+{ static char bufring[NUMBUFS][BUFLEN];
+ static int nbuf=0;
+ if (buffer==0)
+ { buffer=bufring[nbuf]; buflen=BUFLEN; nbuf=(nbuf+1)%NUMBUFS; }
+ if (dim==Dcyl)
+ snprintf(buffer,buflen,"{%f,%f}",t[R],t[Z]);
+ else
+ snprintf(buffer,buflen,"{%f,%f,%f}",t[X],t[Y],t[Z]);
+ return buffer;
+}
+
+/********************************************************************/
+/********************************************************************/
+/********************************************************************/
+grid_volume grid_volume::subvolume(ivec is, ivec ie)
+{
+ if ( !(contains(is) && contains(ie) ) )
+ abort("invalid extents in subvolume");
+ grid_volume sub;
+ sub.dim=dim;
+ sub.a=a;
+ sub.inva=inva;
+ sub.init_subvolume(is, ie);
+ return sub;
+}
+
+void grid_volume::init_subvolume(ivec is, ivec ie)
+{
+ ivec origin(dim,0);
+ LOOP_OVER_DIRECTIONS(dim,d)
+ { num[(int)d] = (ie-is).in_direction(d)/2;
+ origin.set_direction(d, is.in_direction(d));
+ }
+ num_changed();
+ center_origin();
+ shift_origin(origin);
+}
} // namespace meep