vec.cpp

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

  if (which < num_low)
    return split_at_fraction(false, split_point).split_by_effort(num_low,which, Ngv,gv,effort);
  else
    return split_at_fraction(true, split_point).split_by_effort(n-num_low,which-num_low, Ngv,gv,effort);
}

volume volume::split_at_fraction(bool want_high, int numer) const {
  int bestd = -1, bestlen = 1;
  for (int i=0;i<3;i++)
    if (num[i] > bestlen) {
      bestd = i;
      bestlen = num[i];
    }
  if (bestd == -1) {
    for (int i=0;i<3;i++) master_printf("num[%d] = %d\n", i, num[i]);
    abort("Crazy weird splitting error.\n");
  }
  volume retval(dim, a, 1,1,1);
  for (int i=0;i<3;i++) retval.num[i] = num[i];
  if (numer >= num[bestd])
    abort("Aaack bad bug in split_at_fraction.\n");
  direction d = (direction) bestd;
  if (dim == Dcyl && d == X) d = R;
  retval.set_origin(io);
  if (want_high)
    retval.shift_origin(d,numer*2);

  if (want_high) retval.num[bestd] -= numer;
  else retval.num[bestd] = numer;
  retval.num_changed();
  return retval;
}

// Halve the volume for symmetry exploitation...must contain icenter!
volume volume::halve(direction d) const {
  volume retval(*this);
  // note that icenter-io is always even by construction of volume::icenter
  retval.set_num_direction(d, (icenter().in_direction(d) 
			       - io.in_direction(d)) / 2);
  return retval;
}

volume volume::pad(direction d) const {
  volume v(*this);
  v.pad_self(d);
  return v;
}

void volume::pad_self(direction d) {
  num[d%3]+=2; // Pad in both directions by one grid point.
  num_changed();
  shift_origin(d, -2);
}

ivec volume::icenter() const {
  /* Find the center of the user's cell.  This will be used as the
     symmetry point, and therefore icenter-io must be *even*
     in all components in order that rotations preserve the Yee lattice. */
  switch (dim) {
  case D1: return io + ivec(nz()).round_up_to_even();
  case D2: return io + ivec(nx(), ny()).round_up_to_even();
  case D3: return io + ivec(nx(), ny(), nz()).round_up_to_even();
  case Dcyl: return io + iveccyl(0, nz()).round_up_to_even();
  }
  abort("Can't do symmetry with these dimensions.\n");
  return ivec(0); // This is never reached.
}

vec volume::center() const {
  return operator[](icenter());
}

symmetry rotate4(direction axis, const volume &v) {
  symmetry s = identity();
  if (axis > 2) abort("Can only rotate4 in 2D or 3D.\n");
  s.g = 4;
  FOR_DIRECTIONS(d) {
    s.S[d].d = d;
    s.S[d].flipped = false;
  }
  s.S[(axis+1)%3].d = (direction)((axis+2)%3);
  s.S[(axis+1)%3].flipped = true;
  s.S[(axis+2)%3].d = (direction)((axis+1)%3);
  s.symmetry_point = v.center();
  s.i_symmetry_point = v.icenter();
  return s;
}

symmetry rotate2(direction axis, const volume &v) {
  symmetry s = identity();
  if (axis > 2) abort("Can only rotate2 in 2D or 3D.\n");
  s.g = 2;
  s.S[(axis+1)%3].flipped = true;
  s.S[(axis+2)%3].flipped = true;
  s.symmetry_point = v.center();
  s.i_symmetry_point = v.icenter();
  return s;
}

symmetry mirror(direction axis, const volume &v) {
  symmetry s = identity();
  s.g = 2;
  s.S[axis].flipped = true;
  s.symmetry_point = v.center();
  s.i_symmetry_point = v.icenter();
  return s;
}

symmetry r_to_minus_r_symmetry(double m) {
  symmetry s = identity();
  s.g = 2;
  s.S[R].flipped = true;
  s.S[P].flipped = true;
  s.symmetry_point = zero_vec(Dcyl);
  s.i_symmetry_point = zero_ivec(Dcyl);
  if (m == int(m)) // phase is purely real (+/- 1) when m an integer
    s.ph = (int(m) & 1) ? -1.0 : 1.0;
  else
    s.ph = polar(1.0, m * pi); // general case
  return s;
}

symmetry identity() {
  return symmetry();
}

symmetry::symmetry() {
  g = 1;
  ph = 1.0;
  FOR_DIRECTIONS(d) {
    S[d].d = d;
    S[d].flipped = false;
  }
  next = NULL;
}

symmetry::symmetry(const symmetry &s) {
  g = s.g;
  FOR_DIRECTIONS(d) {
    S[d].d = s.S[d].d;
    S[d].flipped = s.S[d].flipped;
  }
  ph = s.ph;
  symmetry_point = s.symmetry_point;
  i_symmetry_point = s.i_symmetry_point;
  if (s.next) next = new symmetry(*s.next);
  else next = NULL;
}

void symmetry::operator=(const symmetry &s) {
  g = s.g;
  FOR_DIRECTIONS(d) {
    S[d].d = s.S[d].d;
    S[d].flipped = s.S[d].flipped;
  }
  ph = s.ph;
  symmetry_point = s.symmetry_point;
  i_symmetry_point = s.i_symmetry_point;
  if (s.next) next = new symmetry(*s.next);
  else next = NULL;
}

bool symmetry::operator==(const symmetry &sym) const {
  int gtot = multiplicity();
  if (gtot != sym.multiplicity())
    return false;
  for (int sn = 1; sn < gtot; ++sn)
    FOR_DIRECTIONS(d)
      if (transform(d, sn) != sym.transform(d, sn))
	return false;
  return true;
}

symmetry::~symmetry() {
  delete next;
}

int symmetry::multiplicity() const {
  if (next) return g*next->multiplicity();
  else return g;
}

symmetry symmetry::operator+(const symmetry &b) const {
  // The following optimization ignores identity when adding symmetries
  // together.  This is important because identity has an undefined
  // symmetry point.
  if (multiplicity() == 1) return b;
  else if (b.multiplicity() == 1) return *this;
  symmetry s = *this;
  symmetry *sn = &s;
  for (; sn->next; sn = sn->next) ;
  sn->next = new symmetry(b);
  return s;
}

symmetry symmetry::operator*(complex<double> p) const {
  symmetry s = *this;
  s.ph *= p;
  return s;
}

signed_direction signed_direction::operator*(complex<double> p) {
  signed_direction sd = *this;
  sd.phase *= p;
  return sd;
}

signed_direction symmetry::transform(direction d, int n) const {
  // Returns transformed direction + phase/flip; -n indicates inverse transform
  if (n == 0 || d == NO_DIRECTION) return signed_direction(d);
  int nme, nrest;
  if (n < 0) {
       nme = (g - (-n) % g) % g;
       nrest = -((-n) / g);
  } else {
       nme = n % g;
       nrest = n / g;
  }
  if (nme == 0) {
    if (nrest == 0) return signed_direction(d);
    else return next->transform(d,nrest);
  } else {
    signed_direction sd;
    if (nme == 1) sd = S[d];
    if (S[d].flipped) sd = flip(transform(S[d].d, nme-1));
    else sd = transform(S[d].d, nme-1);

    if (next && nrest) {
      if (sd.flipped) return flip(next->transform(sd.d, nrest))*ph;
      else return next->transform(sd.d, nrest)*ph;
    } else {
      return sd*ph;
    }
  }
}

ivec symmetry::transform(const ivec &ov, int n) const {
  if (n == 0) return ov;
  ivec out = ov;
  LOOP_OVER_DIRECTIONS(ov.dim, d) {
    const signed_direction s = transform(d,n);
    const int sp_d  = i_symmetry_point.in_direction(d);
    const int sp_sd = i_symmetry_point.in_direction(s.d);
    const int delta = ov.in_direction(d) - sp_d;
    if (s.flipped) out.set_direction(s.d, sp_sd - delta);
    else out.set_direction(s.d, sp_sd + delta);
  }
  return out;
}

ivec symmetry::transform_unshifted(const ivec &ov, int n) const {
  if (n == 0) return ov;
  ivec out(ov.dim);
  LOOP_OVER_DIRECTIONS(ov.dim, d) {
    const signed_direction s = transform(d,n);
    if (s.flipped) out.set_direction(s.d, -ov.in_direction(d));
    else out.set_direction(s.d, ov.in_direction(d));
  }
  return out;
}

vec symmetry::transform(const vec &ov, int n) const {
  if (n == 0) return ov;
  vec delta = ov;
  LOOP_OVER_DIRECTIONS(ov.dim, d) {
    const signed_direction s = transform(d,n);
    double deltad = ov.in_direction(d) - symmetry_point.in_direction(d);
    if (s.flipped) delta.set_direction(s.d, -deltad);
    else delta.set_direction(s.d, deltad);
  }
  return symmetry_point + delta;
}

geometric_volume symmetry::transform(const geometric_volume &gv, int n) const {
  return geometric_volume(transform(gv.get_min_corner(),n),
                          transform(gv.get_max_corner(),n));
}

component symmetry::transform(component c, int n) const {
  return direction_component(c,transform(component_direction(c),n).d);
}

derived_component symmetry::transform(derived_component c, int n) const {
  return direction_component(c,transform(component_direction(c),n).d);
}

int symmetry::transform(int c, int n) const {
  return (is_derived(c) ? int(transform(derived_component(c), n))
	  : int(transform(component(c), n)));
}

complex<double> symmetry::phase_shift(component c, int n) const {
  if (c == Dielectric) return 1.0;
  complex<double> phase = transform(component_direction(c),n).phase;
  // flip tells us if we need to flip the sign.  For vectors (E), it is
  // just this simple:
  bool flip = transform(component_direction(c),n).flipped;
  if (is_magnetic(c)) {
    // Because H is a pseudovector, here we have to figure out if the
    // transformation changes the handedness of the basis.
    bool have_one = false, have_two = false;
    FOR_DIRECTIONS(d) {
      if (transform(d,n).flipped) flip = !flip;
      int shift = (transform(d,n).d - d + 6) % 3;
      if (shift == 1) have_one = true;
      if (shift == 2) have_two = true;
    }
    if (have_one && have_two) flip = !flip;
  }
  if (flip) return -phase;
  else return phase;
}

complex<double> symmetry::phase_shift(derived_component c, int n) const {
  if (is_poynting(c)) {
    signed_direction ds = transform(component_direction(c),n);
    complex<double> ph = conj(ds.phase) * ds.phase; // E x H gets |phase|^2
    return (ds.flipped ? -ph : ph);
  }
  else /* energy density */
    return 1.0;
}

complex<double> symmetry::phase_shift(int c, int n) const {
  return (is_derived(c) ? phase_shift(derived_component(c), n)
	  : phase_shift(component(c), n));
}

bool symmetry::is_primitive(const ivec &p) const {
  // This is only correct if p is somewhere on the yee lattice.
  if (multiplicity() == 1) return true;
  for (int i=1;i<multiplicity();i++) {
    const ivec pp = transform(p,i);
    switch (p.dim) {
    case D2:
      if (pp.x()+pp.y() < p.x()+p.y()) return false;
      if (pp.x()+pp.y() == p.x()+p.y() &&
          p.y() > p.x() && pp.y() <= pp.x()) return false;
      break;
    case D3:
      if (pp.x()+pp.y()+pp.z() <  p.x()+p.y()+p.z()) return false;
      if (pp.x()+pp.y()+pp.z() == p.x()+p.y()+p.z() &&
          pp.x()+pp.y()-pp.z() <  p.x()+p.y()-p.z()) return false;
      if (pp.x()+pp.y()+pp.z() == p.x()+p.y()+p.z() &&
          pp.x()+pp.y()-pp.z() == p.x()+p.y()-p.z() &&
          pp.x()-pp.y()-pp.z() <  p.x()-p.y()-p.z()) return false;
      break;
    case D1: case Dcyl:
      if (pp.z() < p.z()) return false;
      break;
    }
  }
  return true;
}

/* given a list of geometric volumes, produce a new list with appropriate
   weights that is minimized according to the symmetry.  */
geometric_volume_list *symmetry::reduce(const geometric_volume_list *gl) const {
  geometric_volume_list *glnew = 0;
  for (const geometric_volume_list *g = gl; g; g = g->next) {
    int sn;
    for (sn = 0; sn < multiplicity(); ++sn) {
      geometric_volume gS(transform(g->gv, sn));
      int cS = transform(g->c, sn);
      geometric_volume_list *gn;
      for (gn = glnew; gn; gn = gn->next)
	if (gn->c == cS && gn->gv.round_float() == gS.round_float())
	  break;
      if (gn) { // found a match
	gn->weight += g->weight * phase_shift(g->c, sn);
	break;
      }
    }
    if (sn == multiplicity() && g->weight != 0.0) { // no match, add to glnew
      geometric_volume_list *gn = 
	new geometric_volume_list(g->gv, g->c, g->weight, glnew);
      glnew = gn;
    }
  }

  // reduce gv's redundant with themselves & delete elements with zero weight:
  geometric_volume_list *gprev = 0, *g = glnew;
  while (g) {
    // first, see if g->gv is redundant with itself
    bool halve[5] = {false,false,false,false,false};
    complex<double> weight = g->weight;
    for (int sn = 1; sn < multiplicity(); ++sn)
      if (g->c == transform(g->c, sn) && 
	  g->gv.round_float() == transform(g->gv, sn).round_float()) {
	LOOP_OVER_DIRECTIONS(g->gv.dim, d)
	  if (transform(d,sn).flipped) {
	    halve[d] = true;
	    break;
	  }
	g->weight += weight * phase_shift(g->c, sn);
      }
    LOOP_OVER_DIRECTIONS(g->gv.dim, d)
      if (halve[d])
	g->gv.set_direction_max(d, g->gv.in_direction_min(d) +
				0.5 * g->gv.in_direction(d));
    
      // now, delete it if it has zero weight
    if (g->weight == 0.0) {
      if (gprev)
	gprev->next = g->next;
      else // g == glnew
	glnew = g->next;
      g->next = 0; // necessary so that g->next is not deleted recursively
      delete g;
      g = gprev ? gprev->next : glnew;
    }
    else
      g = (gprev = g)->next;
  }

  return glnew;
}

/***************************************************************************/

static double poynting_fun(const complex<double> *fields,
		       const vec &loc, void *data_)
{
     (void) loc; // unused
     (void) data_; // unused
     return (real(conj(fields[0]) * fields[1])
	     - real(conj(fields[2])*fields[3]));
}

static double energy_fun(const complex<double> *fields,
		       const vec &loc, void *data_)
{
     (void) loc; // unused
     int nfields = *((int *) data_) / 2;
     double sum = 0;
     for (int k = 0; k < nfields; ++k)
       sum += real(conj(fields[2*k]) * fields[2*k+1]);
     return sum * 0.5;
}

field_rfunction derived_component_func(derived_component c, const volume &v,
				       int &nfields, component cs[6]) {
  switch (c) {
  case Sx: case Sy: case Sz: case Sr: case Sp:
    switch (c) {
    case Sx: cs[0] = Ey; cs[1] = Hz; break;
    case Sy: cs[0] = Ez; cs[1] = Hx; break;
    case Sz: cs[0] = Ex; cs[1] = Hy; break;
    case Sr: cs[0] = Ep; cs[1] = Hz; break;
    case Sp: cs[0] = Ez; cs[1] = Hr; break;
    default: break; // never reached
    }
    nfields = 4;
    cs[2] = direction_component(Ex, component_direction(cs[1]));
    cs[3] = direction_component(Hx, component_direction(cs[0]));
    return poynting_fun;

  case EnergyDensity: case D_EnergyDensity: case H_EnergyDensity:
    nfields = 0;
    if (c != H_EnergyDensity)
      FOR_ELECTRIC_COMPONENTS(c0) if (v.has_field(c0)) {
	cs[nfields++] = c0;
	cs[nfields++] = direction_component(Dx, component_direction(c0));
      }
    if (c != D_EnergyDensity)
      FOR_MAGNETIC_COMPONENTS(c0) if (v.has_field(c0)) {
	cs[nfields++] = c0;
	cs[nfields++] = c0;
      }
    if (nfields > 6) abort("too many field components");
    return energy_fun;

  default: 
    abort("unknown derived_component in derived_component_func");
  }
  return 0;
}

/***************************************************************************/

} // namespace meep