Function src_vol_chunkloop

src/sources.cpp:243–312

Adding source volumes can be treated as a kind of "integration" problem, since we need to loop over all the chunks that intersect the source grid_volume, with appropriate interpolation weights at the boundaries so that the integral of the current is fixed regardless of resolution. Unlike most uses of fields::loop_in_chunks, however, we set use_symmetry=false: we only find the inter

   with the untransformed chunks (since the transformed versions are
   implicit). */
static void src_vol_chunkloop(fields_chunk *fc, int ichunk, component c, ivec is, ivec ie, vec s0,
                              vec s1, vec e0, vec e1, double dV0, double dV1, ivec shift,
                              complex<double> shift_phase, const symmetry &S, int sn, void *data_) {
  src_vol_chunkloop_data *data = (src_vol_chunkloop_data *)data_;

  (void)S;
  (void)sn; // these should be the identity
  (void)dV0;
  (void)dV1; // grid_volume weighting is included in data->amp
  (void)ichunk;

  size_t npts = 1;
  LOOP_OVER_DIRECTIONS(is.dim, d) { npts *= (ie.in_direction(d) - is.in_direction(d)) / 2 + 1; }
  std::vector<ptrdiff_t> index_array(npts);
  std::vector<complex<double> > amps_array(npts);

  complex<double> amp = data->amp * conj(shift_phase);

  direction cd = component_direction(c);

  double inva = fc->gv.inva;
  size_t idx_vol = 0;
  LOOP_OVER_IVECS(fc->gv, is, ie, idx) {
    IVEC_LOOP_ILOC(fc->gv, iloc);
    if (!fc->gv.owns(iloc)) continue;

    IVEC_LOOP_LOC(fc->gv, loc);
    loc += shift * (0.5 * inva);

    vec rel_loc = loc - data->center;
    amps_array[idx_vol] = IVEC_LOOP_WEIGHT(s0, s1, e0, e1, 1) * amp * data->A(rel_loc);

    // check for invalid sources at r=0 in cylindrical coordinates
    if (fc->gv.dim == Dcyl && loc.r() == 0 && amps_array[idx_vol] != 0.0) {
      if (fc->m == 0) {
        if (component_direction(c) == R || component_direction(c) == P)
          meep::abort("Not possible to place a %s source at r=0 in "
                      "cylindrical coordinates for m = 0.",
                      component_name(c));
      }
      else if (fabs(fc->m) == 1.0) {
        if (component_direction(c) == Z)
          meep::abort("Not possible to place a %s source at r=0 in "
                      "cylindrical coordinates for |m| = 1.0.",
                      component_name(c));
      }
      else {
        meep::abort("Not possible to place a source at r=0 in "
                    "cylindrical coordinates for m = %g.",
                    fc->m);
      }
    }

    /* for "D" sources, multiply by epsilon.  FIXME: this is not quite
       right because it doesn't handle non-diagonal chi1inv!
       similarly, for "B" sources, multiply by mu. */
    if (is_D(c) && fc->s->chi1inv[c - Dx + Ex][cd])
      amps_array[idx_vol] /= fc->s->chi1inv[c - Dx + Ex][cd][idx];