loop_in_chunks.cpp

麻省理工的计算光子晶体的程序

   tests/integrate.cpp file.

   The parameters USE_SYMMETRY (default = true) and SNAP_EMPTY_DIMS
   (default = false) are for use with not-quite-integration-like
   operations.  If use_symmetry is false, then we do *not* loop over
   all possible symmetry transformations of the chunks to see if they
   intersect WHERE; we only use chunks that, untransformed, already
   intersect the volume.  If SNAP_EMPTY_DIMS is true, then for empty
   (min = max) dimensions of WHERE, instead of interpolating, we
   "snap" them to the nearest grid point.  */

void fields::loop_in_chunks(field_chunkloop chunkloop, void *chunkloop_data,
			    const geometric_volume &where, 
			    component cgrid,
			    bool use_symmetry, bool snap_empty_dims)
{
  if (coordinate_mismatch(v.dim, cgrid))
    abort("Invalid fields::loop_in_chunks grid type %s for dimensions %s\n",
	  component_name(cgrid), dimension_name(v.dim));
  if (where.dim != v.dim)
    abort("Invalid dimensions %d for WHERE in fields::loop_in_chunks", where.dim);
  
  /*
    We handle looping on an arbitrary component grid by shifting
    to the dielectric grid and then shifting back.  The looping
    coordinates are internally calculated on the odd-indexed
    "dielectric grid", which has the virtue that it is disjoint for
    each chunk and each chunk has enough information to interpolate all
    of its field components onto this grid without communication.
    Another virtue of this grid is that it is invariant under all of
    our symmetry transformations, so we can uniquely decide which
    transformed chunk gets to loop_in_chunks which grid point. 
  */
  vec yee_c(v.yee_shift(Dielectric) - v.yee_shift(cgrid));
  ivec iyee_c(v.iyee_shift(Dielectric) - v.iyee_shift(cgrid));
  geometric_volume wherec(where + yee_c);

  /* Find the corners (is and ie) of the smallest bounding box for
     wherec, on the grid of odd-coordinate ivecs (i.e. the
     "epsilon grid"). */
  ivec is(vec2diel_floor(wherec.get_min_corner(), v.a, zero_ivec(v.dim)));
  ivec ie(vec2diel_ceil(wherec.get_max_corner(), v.a, zero_ivec(v.dim)));
  
  /* Integration weights at boundaries (c.f. long comment at top). */
  vec s0(v.dim), e0(v.dim), s1(v.dim), e1(v.dim);
  LOOP_OVER_DIRECTIONS(v.dim, d) {
    double w0, w1;
    w0 = 1. - wherec.in_direction_min(d)*v.a + 0.5*is.in_direction(d);
    w1 = 1. + wherec.in_direction_max(d)*v.a - 0.5*ie.in_direction(d);
    if (ie.in_direction(d) >= is.in_direction(d) + 3*2) {
      s0.set_direction(d, w0*w0 / 2);
      s1.set_direction(d, 1 - (1-w0)*(1-w0) / 2);
      e0.set_direction(d, w1*w1 / 2);
      e1.set_direction(d, 1 - (1-w1)*(1-w1) / 2);
    }
    else if (ie.in_direction(d) == is.in_direction(d) + 2*2) {
      s0.set_direction(d, w0*w0 / 2);
      s1.set_direction(d, 1 - (1-w0)*(1-w0) / 2 - (1-w1)*(1-w1) / 2);
      e0.set_direction(d, w1*w1 / 2);
      e1.set_direction(d, s1.in_direction(d));
    }
    else if (wherec.in_direction_min(d) == wherec.in_direction_max(d)) {
      if (snap_empty_dims) {
	if (w0 > w1) ie.set_direction(d, is.in_direction(d));
	else is.set_direction(d, ie.in_direction(d));
	wherec.set_direction_min(d, is.in_direction(d) * (0.5*v.inva));
	wherec.set_direction_max(d, is.in_direction(d) * (0.5*v.inva));
	w0 = w1 = 1.0;
      }
      s0.set_direction(d, w0);
      s1.set_direction(d, w1);
      e0.set_direction(d, w1);
      e1.set_direction(d, w0);
    }
    else if (ie.in_direction(d) == is.in_direction(d) + 1*2) {
      s0.set_direction(d, w0*w0 / 2 - (1-w1)*(1-w1) / 2);
      e0.set_direction(d, w1*w1 / 2 - (1-w0)*(1-w0) / 2);
      s1.set_direction(d, e0.in_direction(d));
      e1.set_direction(d, s0.in_direction(d));
    }
    else
      abort("bug: impossible(?) looping boundaries");
  }

  // loop over symmetry transformations of the chunks:
  for (int sn = 0; sn < (use_symmetry ? S.multiplicity() : 1); ++sn) {
    component cS = S.transform(cgrid, -sn);
    ivec iyee_cS(S.transform_unshifted(iyee_c, -sn));

    geometric_volume vS = S.transform(v.surroundings(), sn);
    vec L(v.dim);
    ivec iL(v.dim);

    // n.b. we can't just S.transform(lattice_vector,sn), 'cause of signs
    LOOP_OVER_DIRECTIONS(v.dim, d) {
      direction dS = S.transform(d, -sn).d;
      L.set_direction(d, fabs(lattice_vector(dS).in_direction(dS)));
      iL.set_direction(d, iabs(ilattice_vector(dS).in_direction(dS)));
    }

    // figure out range of lattice shifts for which vS intersects wherec:
    ivec min_ishift(v.dim), max_ishift(v.dim);
    LOOP_OVER_DIRECTIONS(v.dim, d) {
      if (boundaries[High][S.transform(d, -sn).d] == Periodic) {
	min_ishift.set_direction(d, 
	 int(floor((wherec.in_direction_min(d) - vS.in_direction_max(d))
		   / L.in_direction(d))));
	max_ishift.set_direction(d,
	 int(ceil((wherec.in_direction_max(d) - vS.in_direction_min(d))
		  / L.in_direction(d))));
      }
      else {
	min_ishift.set_direction(d, 0);
	max_ishift.set_direction(d, 0);
      }
    }
    
    // loop over lattice shifts
    ivec ishift(min_ishift);
    do {
      complex<double> ph = 1.0;
      vec shift(v.dim, 0.0);
      ivec shifti(v.dim, 0);
      LOOP_OVER_DIRECTIONS(v.dim, d) {
	shift.set_direction(d, L.in_direction(d) * ishift.in_direction(d));
	shifti.set_direction(d, iL.in_direction(d) * ishift.in_direction(d));
	ph *= pow(eikna[d], ishift.in_direction(d));
      }

      for (int i = 0; i < num_chunks; ++i) {
	if (!chunks[i]->is_mine()) continue;
	// Chunk looping boundaries:
	geometric_volume gvS(v.dim);

	if (use_symmetry)
	  gvS = S.transform(chunks[i]->gv, sn);
	else {
	  /* If we're not using symmetry, it's because (as in src_vol)
	     we don't care about correctly counting the points in the
	     volume.  Rather, we just want to make sure to get *all*
	     of the chunk points that intersect where.  Hence, add a little
	     padding to make sure we don't miss any points due to rounding. */
	  vec pad(one_ivec(v.dim) * v.inva * 1e-3);
	  gvS = geometric_volume(chunks[i]->v.loc(Dielectric,0) - pad,
				 chunks[i]->v.loc(Dielectric,
						  chunks[i]->v.ntot()-1) +pad);
	}

	ivec iscS(max(is-shifti, vec2diel_ceil(gvS.get_min_corner(),
					       v.a, one_ivec(v.dim) * 2)));
	ivec iecS(min(ie-shifti, vec2diel_floor(gvS.get_max_corner(),
						v.a, zero_ivec(v.dim))));
	if (iscS <= iecS) {
	  // Determine weights at chunk looping boundaries:
	  ivec isc(S.transform(iscS, -sn)), iec(S.transform(iecS, -sn));
	  vec s0c(v.dim,1.0), s1c(v.dim,1.0), e0c(v.dim,1.0), e1c(v.dim,1.0);
	  iscS += shifti;
	  iecS += shifti;
	  LOOP_OVER_DIRECTIONS(v.dim, d) {
	    direction dS = S.transform(d, sn).d;
	    if (iscS.in_direction(dS) == is.in_direction(dS)) {
	      s0c.set_direction(d, s0.in_direction(dS));
	      s1c.set_direction(d, s1.in_direction(dS));
	    }
	    else if (iscS.in_direction(dS) == is.in_direction(dS) + 2) {
	      s0c.set_direction(d, s1.in_direction(dS));
	    }
	    if (iecS.in_direction(dS) == ie.in_direction(dS)) {
	      e0c.set_direction(d, e0.in_direction(dS));
	      e1c.set_direction(d, e1.in_direction(dS));
	    }
	    else if (iecS.in_direction(dS) == ie.in_direction(dS) - 2) {
	      e0c.set_direction(d, e1.in_direction(dS));
	    }
	    if (iecS.in_direction(dS) == iscS.in_direction(dS)) {
	      double w = min(s0c.in_direction(d), e0c.in_direction(d));
	      s0c.set_direction(d, w); e0c.set_direction(d, w);
	      s1c.set_direction(d, w); e1c.set_direction(d, w);
	    }
	    else if (iecS.in_direction(dS) == iscS.in_direction(dS) + 1*2) {
	      double w = min(s0c.in_direction(d), e1c.in_direction(d));
	      s0c.set_direction(d, w); e1c.set_direction(d, w);
	      w = min(s1c.in_direction(d), e0c.in_direction(d));
	      s1c.set_direction(d, w); e0c.set_direction(d, w);
	    }
	    else if (iecS.in_direction(dS) == iscS.in_direction(dS) + 2*2) {
	      double w = min(s1c.in_direction(d), e1c.in_direction(d));
	      s1c.set_direction(d, w); e1c.set_direction(d, w);
	    }

	    // swap endpoints/weights if in wrong order due to S.transform
	    if (isc.in_direction(d) > iec.in_direction(d)) {
	      int iswap = isc.in_direction(d);
	      isc.set_direction(d, iec.in_direction(d));
	      iec.set_direction(d, iswap);
	      double swap = s0c.in_direction(d);
	      s0c.set_direction(d, e0c.in_direction(d));
	      e0c.set_direction(d, swap);
	      swap = s1c.in_direction(d);
	      s1c.set_direction(d, e1c.in_direction(d));
	      e1c.set_direction(d, swap);
	    }
	  }
	  

	  // Determine integration "volumes" dV0 and dV1;
	  double dV0 = 1.0, dV1 = 0.0;
	  LOOP_OVER_DIRECTIONS(v.dim, d)
	    if (wherec.in_direction(d) > 0.0)
	      dV0 *= v.inva;
	  if (v.dim == Dcyl) {
	    dV1 = dV0 * 2*pi * v.inva;
	    dV0 *= 2*pi * fabs((S.transform(chunks[i]->v[isc], sn) + shift
				- yee_c).in_direction(R));
	  }
	 
	  chunkloop(chunks[i], i, cS,
		    isc - iyee_cS, iec - iyee_cS,
		    s0c, s1c, e0c, e1c,
		    dV0, dV1,
		    shifti, ph,
		    S, sn,
		    chunkloop_data);
	}
      }
      
      
      LOOP_OVER_DIRECTIONS(v.dim, d) {
	if (ishift.in_direction(d) + 1 <= max_ishift.in_direction(d)) {
	  ishift.set_direction(d, ishift.in_direction(d) + 1);
	  break;
	}
	ishift.set_direction(d, min_ishift.in_direction(d));
      }
    } while (ishift != min_ishift);
  }
}

} // namespace meep