"Fossies" - the Fresh Open Source Software Archive

Member "meep-1.10.0/src/array_slice.cpp" (5 Jun 2019, 29210 Bytes) of package /linux/privat/meep-1.10.0.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "array_slice.cpp" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 1.9.0_vs_1.10.0.

    1 /* Copyright (C) 2005-2019 Massachusetts Institute of Technology
    2 %
    3 %  This program is free software; you can redistribute it and/or modify
    4 %  it under the terms of the GNU General Public License as published by
    5 %  the Free Software Foundation; either version 2, or (at your option)
    6 %  any later version.
    7 %
    8 %  This program is distributed in the hope that it will be useful,
    9 %  but WITHOUT ANY WARRANTY; without even the implied warranty of
   10 %  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   11 %  GNU General Public License for more details.
   12 %
   13 %  You should have received a copy of the GNU General Public License
   14 %  along with this program; if not, write to the Free Software Foundation,
   15 %  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
   16 */
   17 
   18 /* create and return arrays of field components on user-specified
   19    spatial slices. Uses fields::loop_in_chunks analogous to
   20    h5fields.cpp
   21 */
   22 
   23 #include <stdio.h>
   24 #include <string.h>
   25 #include <math.h>
   26 
   27 #include "meep_internals.hpp"
   28 
   29 using namespace std;
   30 
   31 typedef complex<double> cdouble;
   32 
   33 namespace meep {
   34 
   35 /***************************************************************************/
   36 
   37 typedef struct {
   38 
   39   // information related to the volume covered by the
   40   // array slice (its size, etcetera)
   41   // these fields are filled in by get_array_slice_dimensions
   42   // if the data parameter is non-null
   43   ivec min_corner, max_corner;
   44   int num_chunks;
   45   int rank;
   46   direction ds[3];
   47   size_t slice_size;
   48 
   49   // if non-null, min_max_loc[0,1] are filled in by get_array_slice_dimensions_chunkloop
   50   // with the (coordinate-wise) minimum and maximum grid points encountered
   51   // in looping over the slice region.
   52   vec *min_max_loc;
   53 
   54   // the function to output and related info (offsets for averaging, etc.)
   55   // note: either fun *or* rfun should be non-NULL (not both)
   56   field_function fun;
   57   field_rfunction rfun;
   58   void *fun_data;
   59   std::vector<component> components;
   60 
   61   void *vslice;
   62 
   63   // temporary internal storage buffers
   64   component *cS;
   65   cdouble *ph;
   66   cdouble *fields;
   67   ptrdiff_t *offsets;
   68 
   69   int ninveps;
   70   component inveps_cs[3];
   71   direction inveps_ds[3];
   72 
   73   int ninvmu;
   74   component invmu_cs[3];
   75   direction invmu_ds[3];
   76 
   77 } array_slice_data;
   78 
   79 #define UNUSED(x) (void)x // silence compiler warnings
   80 
   81 /* passthrough field function equivalent to component_fun in h5fields.cpp */
   82 static cdouble default_field_func(const cdouble *fields, const vec &loc, void *data_) {
   83   (void)loc;   // unused
   84   (void)data_; // unused
   85   return fields[0];
   86 }
   87 
   88 static double default_field_rfunc(const cdouble *fields, const vec &loc, void *data_) {
   89   (void)loc;   // unused
   90   (void)data_; // unused
   91   return real(fields[0]);
   92 }
   93 
   94 /***************************************************************/
   95 /* callback function passed to loop_in_chunks to compute       */
   96 /* dimensions of array slice                                   */
   97 /***************************************************************/
   98 static void get_array_slice_dimensions_chunkloop(fields_chunk *fc, int ichnk, component cgrid,
   99                                                  ivec is, ivec ie, vec s0, vec s1, vec e0, vec e1,
  100                                                  double dV0, double dV1, ivec shift,
  101                                                  complex<double> shift_phase, const symmetry &S,
  102                                                  int sn, void *data_) {
  103   UNUSED(ichnk);
  104   UNUSED(cgrid);
  105   UNUSED(s0);
  106   UNUSED(s1);
  107   UNUSED(e0);
  108   UNUSED(e1);
  109   UNUSED(dV0);
  110   UNUSED(dV1);
  111   UNUSED(shift_phase);
  112 
  113   array_slice_data *data = (array_slice_data *)data_;
  114   ivec isS = S.transform(is, sn) + shift;
  115   ivec ieS = S.transform(ie, sn) + shift;
  116   data->min_corner = min(data->min_corner, min(isS, ieS));
  117   data->max_corner = max(data->max_corner, max(isS, ieS));
  118   data->num_chunks++;
  119 
  120   if (data->min_max_loc) {
  121     vec rshift(shift * (0.5 * fc->gv.inva));
  122     LOOP_OVER_IVECS(fc->gv, is, ie, idx) {
  123       IVEC_LOOP_LOC(fc->gv, loc);
  124       loc = S.transform(loc, sn) + rshift;
  125       data->min_max_loc[0] = min(loc, data->min_max_loc[0]);
  126       data->min_max_loc[1] = max(loc, data->min_max_loc[1]);
  127     }
  128   }
  129 }
  130 
  131 /***************************************************************/
  132 /* chunkloop for get_source_slice ******************************/
  133 /***************************************************************/
  134 typedef struct {
  135  component source_component;
  136  ivec slice_imin, slice_imax;
  137  cdouble *slice;
  138 } source_slice_data;
  139 
  140 bool in_range(int imin, int i, int imax)
  141 { return (imin<=i && i<=imax); }
  142 
  143 bool in_subgrid(ivec ivmin, ivec iv, ivec ivmax)
  144 { LOOP_OVER_DIRECTIONS(iv.dim,d)
  145    if ( !in_range(ivmin.in_direction(d),iv.in_direction(d),ivmax.in_direction(d)) )
  146     return false;
  147   return true;
  148 }
  149 
  150 static void get_source_slice_chunkloop(fields_chunk *fc, int ichnk, component cgrid,
  151                                        ivec is, ivec ie, vec s0, vec s1, vec e0, vec e1,
  152                                        double dV0, double dV1, ivec shift,
  153                                        complex<double> shift_phase, const symmetry &S,
  154                                        int sn, void *data_) {
  155 
  156   UNUSED(ichnk); UNUSED(cgrid); UNUSED(is); UNUSED(ie); UNUSED(s0); UNUSED(s1);
  157   UNUSED(e0); UNUSED(e1); UNUSED(dV0); UNUSED(dV1); UNUSED(shift_phase);
  158 
  159   source_slice_data *data = (source_slice_data *)data_;
  160   ivec slice_imin = data->slice_imin, slice_imax=data->slice_imax;
  161   ndim dim=fc->gv.dim;
  162 
  163   // the following works in all cases except cylindrical coordinates
  164   ptrdiff_t NY = 1, NZ = 1;
  165   if (has_direction(dim, Z)) NZ = ((slice_imax - slice_imin).in_direction(Z) / 2) + 1;
  166   if (has_direction(dim, Y)) NY = ((slice_imax - slice_imin).in_direction(Y) / 2) + 1;
  167 
  168   for (int ft = 0; ft < NUM_FIELD_TYPES; ft++)
  169     for (src_vol *s = fc->sources[ft]; s; s = s->next) {
  170       component cS = S.transform(data->source_component, -sn);
  171       if (s->c != cS) continue;
  172 
  173       // loop over point sources in this src_vol. for each point source,
  174       // the src_vol stores the amplitude and the global index of the
  175       // symmetry-parent grid point, from which we need to compute the
  176       // local index of the symmetry-child grid point within this
  177       // slice (that is, if it even lies within the slice)
  178       for (size_t npt = 0; npt < s->npts; npt++) {
  179         cdouble amp = s->A[npt];
  180         ptrdiff_t chunk_index = s->index[npt];
  181         ivec iloc_parent = fc->gv.iloc(Dielectric, chunk_index);
  182         ivec iloc_child  = S.transform(iloc_parent,sn) + shift;
  183         if (!in_subgrid(slice_imin, iloc_child, slice_imax)) continue; // source point outside slice
  184         ivec slice_offset = iloc_child - slice_imin;
  185         ptrdiff_t slice_index = 0;
  186         // the following works to set the slice_index in all cases except cylindrical coordinates
  187         if (has_direction(dim, Z)) slice_index += slice_offset.in_direction(Z) / 2;
  188         if (has_direction(dim, Y)) slice_index += NZ * slice_offset.in_direction(Y) / 2;
  189         if (has_direction(dim, X)) slice_index += NY * NZ * slice_offset.in_direction(X) / 2;
  190         data->slice[slice_index] = amp;
  191       }
  192     }
  193 }
  194 
  195 /***************************************************************/
  196 /* callback function passed to loop_in_chunks to fill array slice */
  197 /***************************************************************/
  198 static void get_array_slice_chunkloop(fields_chunk *fc, int ichnk, component cgrid, ivec is,
  199                                       ivec ie, vec s0, vec s1, vec e0, vec e1, double dV0,
  200                                       double dV1, ivec shift, complex<double> shift_phase,
  201                                       const symmetry &S, int sn, void *data_) {
  202   UNUSED(ichnk);
  203   UNUSED(cgrid);
  204   UNUSED(s0);
  205   UNUSED(s1);
  206   UNUSED(e0);
  207   UNUSED(e1);
  208   UNUSED(dV0);
  209   UNUSED(dV1);
  210   array_slice_data *data = (array_slice_data *)data_;
  211 
  212   //-----------------------------------------------------------------------//
  213   // Find output chunk dimensions and strides, etc.
  214   //-----------------------------------------------------------------------//
  215 
  216   int start[3] = {0, 0, 0}, count[3] = {1, 1, 1};
  217   ptrdiff_t offset[3] = {0, 0, 0};
  218 
  219   ivec isS = S.transform(is, sn) + shift;
  220   ivec ieS = S.transform(ie, sn) + shift;
  221 
  222   // figure out what yucky_directions (in LOOP_OVER_IVECS)
  223   // correspond to what directions in the transformed vectors (in output).
  224   ivec permute(zero_ivec(fc->gv.dim));
  225   for (int i = 0; i < 3; ++i)
  226     permute.set_direction(fc->gv.yucky_direction(i), i);
  227   permute = S.transform_unshifted(permute, sn);
  228   LOOP_OVER_DIRECTIONS(permute.dim, d) { permute.set_direction(d, abs(permute.in_direction(d))); }
  229 
  230   // compute the size of the chunk to output, and its strides etc.
  231   size_t slice_size = 1;
  232   for (int i = 0; i < data->rank; ++i) {
  233     direction d = data->ds[i];
  234     int isd = isS.in_direction(d), ied = ieS.in_direction(d);
  235     start[i] = (min(isd, ied) - data->min_corner.in_direction(d)) / 2;
  236     count[i] = abs(ied - isd) / 2 + 1;
  237     slice_size *= count[i];
  238     if (ied < isd) offset[permute.in_direction(d)] = count[i] - 1;
  239   }
  240 
  241   // slightly confusing: for array_slice, in contrast to
  242   // h5fields, strides are computed using the dimensions of
  243   // the full array slice, not the dimensions of the chunk.
  244   size_t dims[3] = {1, 1, 1};
  245   for (int i = 0; i < data->rank; i++) {
  246     direction d = data->ds[i];
  247     dims[i] = (data->max_corner.in_direction(d) - data->min_corner.in_direction(d)) / 2 + 1;
  248   }
  249 
  250   ptrdiff_t stride[3] = {1, 1, 1};
  251   for (int i = 0; i < data->rank; ++i) {
  252     direction d = data->ds[i];
  253     int j = permute.in_direction(d);
  254     for (int k = i + 1; k < data->rank; ++k)
  255       stride[j] *= dims[k];
  256     offset[j] *= stride[j];
  257     if (offset[j]) stride[j] *= -1;
  258   }
  259 
  260   // sco="slice chunk offset"
  261   ptrdiff_t sco = start[0] * dims[1] * dims[2] + start[1] * dims[2] + start[2];
  262 
  263   //-----------------------------------------------------------------------//
  264   // Otherwise proceed to compute the function of field components to be   //
  265   // tabulated on the slice, exactly as in fields::integrate.              //
  266   //-----------------------------------------------------------------------//
  267   double *slice = 0;
  268   cdouble *zslice = 0;
  269   bool complex_data = (data->rfun == 0);
  270   if (complex_data)
  271     zslice = (cdouble *)data->vslice;
  272   else
  273     slice = (double *)data->vslice;
  274 
  275   ptrdiff_t *off = data->offsets;
  276   component *cS = data->cS;
  277   complex<double> *fields = data->fields, *ph = data->ph;
  278   const component *iecs = data->inveps_cs;
  279   const direction *ieds = data->inveps_ds;
  280   ptrdiff_t ieos[6];
  281   const component *imcs = data->invmu_cs;
  282   const direction *imds = data->invmu_ds;
  283   ptrdiff_t imos[6];
  284   int num_components = data->components.size();
  285 
  286   for (int i = 0; i < num_components; ++i) {
  287     cS[i] = S.transform(data->components[i], -sn);
  288     if (cS[i] == Dielectric || cS[i] == Permeability)
  289       ph[i] = 1.0;
  290     else {
  291       fc->gv.yee2cent_offsets(cS[i], off[2 * i], off[2 * i + 1]);
  292       ph[i] = shift_phase * S.phase_shift(cS[i], sn);
  293     }
  294   }
  295   for (int k = 0; k < data->ninveps; ++k)
  296     fc->gv.yee2cent_offsets(iecs[k], ieos[2 * k], ieos[2 * k + 1]);
  297   for (int k = 0; k < data->ninvmu; ++k)
  298     fc->gv.yee2cent_offsets(imcs[k], imos[2 * k], imos[2 * k + 1]);
  299 
  300   vec rshift(shift * (0.5 * fc->gv.inva));
  301   // main loop over all grid points owned by this field chunk.
  302   LOOP_OVER_IVECS(fc->gv, is, ie, idx) {
  303 
  304     // get real-space coordinates of grid point, taking into
  305     // account the complications of symmetries.
  306     IVEC_LOOP_LOC(fc->gv, loc);
  307     loc = S.transform(loc, sn) + rshift;
  308 
  309     // interpolate fields at the four nearest grid points
  310     // to get the value of the field component for this point
  311     for (int i = 0; i < num_components; ++i) {
  312       // special case for fetching grid point coordinates and weights
  313       if (cS[i] == NO_COMPONENT) {
  314         fields[i] = IVEC_LOOP_WEIGHT(s0, s1, e0, e1, dV0 + dV1 * loop_i2);
  315       } else if (cS[i] == Dielectric) {
  316         double tr = 0.0;
  317         for (int k = 0; k < data->ninveps; ++k) {
  318           const realnum *ie = fc->s->chi1inv[iecs[k]][ieds[k]];
  319           if (ie)
  320             tr += (ie[idx] + ie[idx + ieos[2 * k]] + ie[idx + ieos[1 + 2 * k]] +
  321                    ie[idx + ieos[2 * k] + ieos[1 + 2 * k]]);
  322           else
  323             tr += 4; // default inveps == 1
  324         }
  325         fields[i] = (4 * data->ninveps) / tr;
  326       } else if (cS[i] == Permeability) {
  327         double tr = 0.0;
  328         for (int k = 0; k < data->ninvmu; ++k) {
  329           const realnum *im = fc->s->chi1inv[imcs[k]][imds[k]];
  330           if (im)
  331             tr += (im[idx] + im[idx + imos[2 * k]] + im[idx + imos[1 + 2 * k]] +
  332                    im[idx + imos[2 * k] + imos[1 + 2 * k]]);
  333           else
  334             tr += 4; // default invmu == 1
  335         }
  336         fields[i] = (4 * data->ninvmu) / tr;
  337       } else {
  338         double f[2];
  339         for (int k = 0; k < 2; ++k)
  340           if (fc->f[cS[i]][k])
  341             f[k] = 0.25 * (fc->f[cS[i]][k][idx] + fc->f[cS[i]][k][idx + off[2 * i]] +
  342                            fc->f[cS[i]][k][idx + off[2 * i + 1]] +
  343                            fc->f[cS[i]][k][idx + off[2 * i] + off[2 * i + 1]]);
  344           else
  345             f[k] = 0;
  346         fields[i] = complex<double>(f[0], f[1]) * ph[i];
  347       }
  348     }
  349 
  350     // compute the index into the array for this grid point and store the result of the computation
  351     ptrdiff_t idx2 =
  352         sco + ((((offset[0] + offset[1] + offset[2]) + loop_i1 * stride[0]) + loop_i2 * stride[1]) +
  353                loop_i3 * stride[2]);
  354 
  355     if (complex_data)
  356       zslice[idx2] = data->fun(fields, loc, data->fun_data);
  357     else
  358       slice[idx2] = data->rfun(fields, loc, data->fun_data);
  359 
  360   } // LOOP_OVER_IVECS
  361 }
  362 
  363 /***************************************************************/
  364 /* repeatedly call sum_to_all to consolidate all entries of    */
  365 /* an array on all cores.                                      */
  366 /***************************************************************/
  367 #define BUFSIZE 1 << 16 // use 64k buffer
  368 double *array_to_all(double *array, size_t array_size) {
  369   double *buffer = new double[BUFSIZE];
  370   ptrdiff_t offset = 0;
  371   size_t remaining = array_size;
  372   while (remaining != 0) {
  373    size_t xfer_size = (remaining > BUFSIZE ? BUFSIZE : remaining);
  374    sum_to_all(array + offset, buffer, xfer_size);
  375    memcpy(array + offset, buffer, xfer_size*sizeof(double));
  376    remaining -= xfer_size;
  377    offset += xfer_size;
  378   }
  379   delete[] buffer;
  380   return array;
  381 }
  382 
  383 cdouble *array_to_all(cdouble *array, size_t array_size) {
  384  return (cdouble *)array_to_all( (double *)array, 2*array_size);
  385 }
  386 
  387 
  388 /***************************************************************/
  389 /* given a volume, fill in the dims[] and dirs[] arrays        */
  390 /* describing the array slice needed to store field data for   */
  391 /* all grid points in the volume.                              */
  392 /*                                                             */
  393 /* return value is rank of array slice.                        */
  394 /*                                                             */
  395 /* if caller_data is non-NULL, it should point to a            */
  396 /* caller-allocated array_slice_data structure which will be   */
  397 /* initialized appopriately for subsequent use in              */
  398 /* get_array_slice.                                            */
  399 /***************************************************************/
  400 int fields::get_array_slice_dimensions(const volume &where, size_t dims[3], direction dirs[3],
  401                                        bool collapse_empty_dimensions, bool snap_empty_dimensions,
  402                                        vec *min_max_loc, void *caller_data) {
  403   am_now_working_on(FieldOutput);
  404 
  405   // use a local data structure if the caller didn't provide one
  406   array_slice_data local_data;
  407   array_slice_data *data = (array_slice_data *)caller_data;
  408   if (data == 0) data = &local_data;
  409 
  410   data->min_corner = gv.round_vec(where.get_max_corner()) + one_ivec(gv.dim);
  411   data->max_corner = gv.round_vec(where.get_min_corner()) - one_ivec(gv.dim);
  412   data->num_chunks = 0;
  413 
  414   data->min_max_loc = min_max_loc;
  415   vec *min_loc = 0, *max_loc = 0;
  416   if (min_max_loc) {
  417     min_loc = min_max_loc + 0;
  418     max_loc = min_max_loc + 1;
  419     LOOP_OVER_DIRECTIONS(gv.dim, d) {
  420       min_loc->set_direction(d, +infinity);
  421       max_loc->set_direction(d, -infinity);
  422     }
  423   }
  424 
  425   bool use_symmetry = true;
  426   loop_in_chunks(get_array_slice_dimensions_chunkloop, (void *)data, where, Centered, use_symmetry,
  427                  snap_empty_dimensions);
  428 
  429   data->min_corner = -max_to_all(-data->min_corner); // i.e., min_to_all
  430   data->max_corner = max_to_all(data->max_corner);
  431   if (min_max_loc) LOOP_OVER_DIRECTIONS(gv.dim, d) {
  432       min_loc->set_direction(d, -1.0 * max_to_all(-1.0 * min_loc->in_direction(d)));
  433       max_loc->set_direction(d, max_to_all(max_loc->in_direction(d)));
  434     }
  435   data->num_chunks = sum_to_all(data->num_chunks);
  436   if (data->num_chunks == 0 || !(data->min_corner <= data->max_corner))
  437     return 0; // no data to write;
  438 
  439   int rank = 0;
  440   size_t slice_size = 1;
  441   LOOP_OVER_DIRECTIONS(gv.dim, d) {
  442     if (rank >= 3) abort("too many dimensions in array_slice");
  443     size_t n = (data->max_corner.in_direction(d) - data->min_corner.in_direction(d)) / 2 + 1;
  444     if (where.in_direction(d) == 0.0 && collapse_empty_dimensions) n = 1;
  445     if (n > 1) {
  446       data->ds[rank] = d;
  447       dims[rank++] = n;
  448       slice_size *= n;
  449     }
  450   }
  451   for (int r = 0; r < rank; r++)
  452     dirs[r] = data->ds[r];
  453   data->rank = rank;
  454   data->slice_size = slice_size;
  455   finished_working();
  456 
  457   return rank;
  458 }
  459 
  460 /**********************************************************************/
  461 /* precisely one of fun, rfun, should be non-NULL                     */
  462 /**********************************************************************/
  463 void *fields::do_get_array_slice(const volume &where, std::vector<component> components,
  464                                  field_function fun, field_rfunction rfun, void *fun_data,
  465                                  void *vslice) {
  466   am_now_working_on(FieldOutput);
  467 
  468   /***************************************************************/
  469   /* call get_array_slice_dimensions to get slice dimensions and */
  470   /* partially initialze an array_slice_data struct              */
  471   /***************************************************************/
  472   // by tradition, empty dimensions in time-domain field arrays are *not* collapsed
  473   // TODO make this a caller-specifiable parameter to get_array_slice()?
  474   bool collapse = false, snap = true;
  475   size_t dims[3];
  476   direction dirs[3];
  477   array_slice_data data;
  478   int rank = get_array_slice_dimensions(where, dims, dirs, collapse, snap, 0, &data);
  479   size_t slice_size = data.slice_size;
  480 
  481   bool complex_data = (rfun == 0);
  482   cdouble *zslice;
  483   double *slice;
  484   if (vslice == 0) {
  485     if (complex_data) {
  486       zslice = new cdouble[slice_size];
  487       memset(zslice, 0, slice_size * sizeof(cdouble));
  488       vslice = (void *)zslice;
  489     } else {
  490       slice = new double[slice_size];
  491       memset(slice, 0, slice_size * sizeof(double));
  492       vslice = (void *)slice;
  493     }
  494   }
  495 
  496   data.vslice = vslice;
  497   data.fun = fun;
  498   data.rfun = rfun;
  499   data.fun_data = fun_data;
  500   data.components = components;
  501   int num_components = components.size();
  502   data.cS = new component[num_components];
  503   data.ph = new cdouble[num_components];
  504   data.fields = new cdouble[num_components];
  505   data.offsets = new ptrdiff_t[2 * num_components];
  506   memset(data.offsets, 0, 2 * num_components * sizeof(ptrdiff_t));
  507 
  508   /* compute inverse-epsilon directions for computing Dielectric fields */
  509   data.ninveps = 0;
  510   bool needs_dielectric = false;
  511   for (int i = 0; i < num_components; ++i)
  512     if (components[i] == Dielectric) {
  513       needs_dielectric = true;
  514       break;
  515     }
  516   if (needs_dielectric) FOR_ELECTRIC_COMPONENTS(c) if (gv.has_field(c)) {
  517       if (data.ninveps == 3) abort("more than 3 field components??");
  518       data.inveps_cs[data.ninveps] = c;
  519       data.inveps_ds[data.ninveps] = component_direction(c);
  520       ++data.ninveps;
  521     }
  522 
  523   /* compute inverse-mu directions for computing Permeability fields */
  524   data.ninvmu = 0;
  525   bool needs_permeability = false;
  526   for (int i = 0; i < num_components; ++i)
  527     if (components[i] == Permeability) {
  528       needs_permeability = true;
  529       break;
  530     }
  531   if (needs_permeability) FOR_MAGNETIC_COMPONENTS(c) if (gv.has_field(c)) {
  532       if (data.ninvmu == 3) abort("more than 3 field components??");
  533       data.invmu_cs[data.ninvmu] = c;
  534       data.invmu_ds[data.ninvmu] = component_direction(c);
  535       ++data.ninvmu;
  536     }
  537 
  538   loop_in_chunks(get_array_slice_chunkloop, (void *)&data, where, Centered, true, true);
  539 
  540   /***************************************************************/
  541   /*consolidate full array on all cores                          */
  542   /***************************************************************/
  543   vslice=(void *)array_to_all( (double *)vslice, (complex_data ? 2:1)*slice_size);
  544 
  545   delete[] data.offsets;
  546   delete[] data.fields;
  547   delete[] data.ph;
  548   delete[] data.cS;
  549   finished_working();
  550 
  551   return vslice;
  552 }
  553 
  554 /***************************************************************/
  555 /* entry points to get_array_slice                             */
  556 /***************************************************************/
  557 double *fields::get_array_slice(const volume &where, std::vector<component> components,
  558                                 field_rfunction rfun, void *fun_data, double *slice) {
  559   return (double *)do_get_array_slice(where, components, 0, rfun, fun_data, (void *)slice);
  560 }
  561 
  562 cdouble *fields::get_complex_array_slice(const volume &where, std::vector<component> components,
  563                                          field_function fun, void *fun_data, cdouble *slice) {
  564   return (cdouble *)do_get_array_slice(where, components, fun, 0, fun_data, (void *)slice);
  565 }
  566 
  567 double *fields::get_array_slice(const volume &where, component c, double *slice) {
  568   std::vector<component> components(1);
  569   components[0] = c;
  570   return (double *)do_get_array_slice(where, components, 0, default_field_rfunc, 0, (void *)slice);
  571 }
  572 
  573 double *fields::get_array_slice(const volume &where, derived_component c, double *slice) {
  574   int nfields;
  575   component carray[12];
  576   field_rfunction rfun = derived_component_func(c, gv, nfields, carray);
  577   std::vector<component> cs(carray, carray + nfields);
  578   return (double *)do_get_array_slice(where, cs, 0, rfun, &nfields, (void *)slice);
  579 }
  580 
  581 cdouble *fields::get_complex_array_slice(const volume &where, component c, cdouble *slice) {
  582   std::vector<component> components(1);
  583   components[0] = c;
  584   return (cdouble *)do_get_array_slice(where, components, default_field_func, 0, 0, (void *)slice);
  585 }
  586 
  587 cdouble *fields::get_source_slice(const volume &where, component source_slice_component,
  588                                   cdouble *slice) {
  589 
  590   size_t dims[3];
  591   direction dirs[3];
  592   vec min_max_loc[2];
  593   bool collapse=false, snap=false;
  594   int rank=get_array_slice_dimensions(where, dims, dirs, collapse, snap, min_max_loc);
  595   size_t slice_size = dims[0] * (rank>=2 ? dims[1] : 1) * (rank==3 ? dims[2] : 1);
  596 
  597   source_slice_data data;
  598   data.source_component=source_slice_component;
  599   data.slice_imin=gv.round_vec(min_max_loc[0]);
  600   data.slice_imax=gv.round_vec(min_max_loc[1]);
  601   data.slice = slice ? slice : new cdouble[slice_size];
  602   if (!data.slice) abort("%s:%i: out of memory (%zu)", __FILE__, __LINE__, slice_size);
  603 
  604   loop_in_chunks(get_source_slice_chunkloop, (void *)&data, where, Centered, true, true);
  605   return array_to_all(data.slice,slice_size);
  606 }
  607 
  608 /**********************************************************************/
  609 /* increment a multi-index <n_1, n_2, ..., n_r> in which n_i runs over*/
  610 /* the range 0 <= n_i < nMax_i. return true if this is the increment  */
  611 /* that brings the multiindex to all zeros.                           */
  612 /**********************************************************************/
  613 bool increment(size_t *n, size_t *nMax, int rank) {
  614   for (rank--; rank >= 0; rank--)
  615     if (++n[rank] < nMax[rank])
  616       return false;
  617     else
  618       n[rank] = 0;
  619   return true;
  620 }
  621 
  622 // data_size = 1,2 for real,complex-valued array
  623 double *collapse_array(double *array, int *rank, size_t dims[3], direction dirs[3], volume where,
  624                        int data_size = 1) {
  625 
  626   /*--------------------------------------------------------------*/
  627   /*- detect empty dimensions and compute rank and strides for    */
  628   /*- collapsed array                                             */
  629   /*--------------------------------------------------------------*/
  630   int full_rank = *rank;
  631   if (full_rank == 0) return array;
  632 
  633   int reduced_rank = 0;
  634   size_t reduced_dims[3], reduced_stride[3] = {1, 1, 1};
  635   direction reduced_dirs[3];
  636   for (int r = 0; r < full_rank; r++) {
  637     if (where.in_direction(dirs[r]) == 0.0)
  638       reduced_stride[r] = 0; // degenerate dimension, to be collapsed
  639     else {
  640       reduced_dirs[reduced_rank] = dirs[r];
  641       reduced_dims[reduced_rank++] = dims[r];
  642     }
  643   }
  644   if (reduced_rank == full_rank) return array; // nothing to collapse
  645 
  646   /*--------------------------------------------------------------*/
  647   /*- set up strides into full and reduced arrays                 */
  648   /*--------------------------------------------------------------*/
  649   size_t stride[3] = {1, 1, 1}; // non-reduced array strides
  650   if (full_rank == 2)
  651     stride[0] = dims[1]; // rstride is already all set in this case
  652   else if (full_rank == 3) {
  653     stride[0] = dims[1] * dims[2];
  654     stride[1] = dims[2];
  655     if (reduced_stride[0] != 0)
  656       reduced_stride[0] = reduced_dims[1];
  657     else if (reduced_stride[1] != 0)
  658       reduced_stride[1] = reduced_dims[1];
  659     // else: two degenerate dimensions->reduced array is 1-diml, no strides needed
  660   }
  661 
  662   /*--------------------------------------------------------------*/
  663   /*- allocate reduced array and compress full array into it     -*/
  664   /*--------------------------------------------------------------*/
  665   size_t reduced_grid_size = reduced_dims[0] * (reduced_rank == 2 ? reduced_dims[1] : 1);
  666   size_t reduced_array_size = data_size * reduced_grid_size;
  667   double *reduced_array = new double[reduced_array_size];
  668   if (!reduced_array) abort("%s:%i: out of memory (%zu)", __FILE__, __LINE__, reduced_array_size);
  669   memset(reduced_array, 0, reduced_array_size * sizeof(double));
  670 
  671   size_t n[3] = {0, 0, 0};
  672   do {
  673     size_t index = n[0] * stride[0] + n[1] * stride[1] + n[2] * stride[2];
  674     size_t rindex = n[0] * reduced_stride[0] + n[1] * reduced_stride[1] + n[2] * reduced_stride[2];
  675     for (int i = 0; i < data_size; i++)
  676       reduced_array[data_size * rindex + i] += array[data_size * index + i];
  677   } while (!increment(n, dims, full_rank));
  678 
  679   *rank = reduced_rank;
  680   for (int r = 0; r < reduced_rank; r++) {
  681     dims[r] = reduced_dims[r];
  682     dirs[r] = reduced_dirs[r];
  683   }
  684   delete[] array;
  685   return reduced_array;
  686 }
  687 
  688 cdouble *collapse_array(cdouble *array, int *rank, size_t dims[3], direction dirs[3],
  689                         volume where) {
  690   return (cdouble *)collapse_array((double *)array, rank, dims, dirs, where, 2);
  691 }
  692 
  693 /***************************************************************/
  694 /***************************************************************/
  695 /***************************************************************/
  696 std::vector<double> fields::get_array_metadata(const volume &where, bool collapse_empty_dimensions,
  697                                                bool snap_empty_dimensions) {
  698 
  699   /* get extremal corners of subgrid and array of weights, collapsed if necessary */
  700   size_t dims[3];
  701   direction dirs[3];
  702   vec min_max_loc[2]; // extremal points in subgrid
  703   int rank = get_array_slice_dimensions(where, dims, dirs, false /*collapse_empty_dimensions*/,
  704                                         snap_empty_dimensions, min_max_loc);
  705   int full_rank = rank;
  706   direction full_dirs[3];
  707   for (int fr = 0; fr < rank; fr++)
  708     full_dirs[fr] = dirs[fr];
  709 
  710   double *weights = get_array_slice(where, NO_COMPONENT);
  711   if (collapse_empty_dimensions) weights = collapse_array(weights, &rank, dims, dirs, where);
  712 
  713   /* get length and endpoints of x,y,z tics arrays */
  714   size_t nxyz[3] = {1, 1, 1};
  715   double xyzmin[3] = {0.0, 0.0, 0.0}, xyzmax[3] = {0.0, 0.0, 0.0};
  716   for (int fr = 0, rr = 0; fr < full_rank; fr++) {
  717     direction d = full_dirs[fr];
  718     int nd = d - X;
  719     if (where.in_direction(d) == 0.0 && collapse_empty_dimensions) {
  720       xyzmin[nd] = xyzmax[nd] = where.in_direction_min(d);
  721       nxyz[nd] = 1;
  722     } else {
  723       nxyz[nd] = dims[rr++];
  724       xyzmin[nd] = min_max_loc[0].in_direction(d);
  725       xyzmax[nd] = min_max_loc[1].in_direction(d);
  726     }
  727   }
  728 
  729   /* pack all data into a single vector with each tics array preceded by its */
  730   /* length: [ NX, xtics[:], NY, ytics[:], NZ, ztics[:], weights[:] ]        */
  731   std::vector<double> xyzw;
  732   for (int nd = 0; nd < 3; nd++) {
  733     xyzw.push_back((double)nxyz[nd]);
  734     for (size_t n = 0; n < nxyz[nd]; n++)
  735       xyzw.push_back(xyzmin[nd] + n * gv.inva);
  736   }
  737   for (unsigned nw = 0; nw < (nxyz[0] * nxyz[1] * nxyz[2]); nw++)
  738     xyzw.push_back(weights[nw]);
  739 
  740   return xyzw;
  741 
  742 } // get_array_metadata
  743 
  744 } // namespace meep