vec.hpp

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

// -*- C++ -*-
/* Copyright (C) 2006 Massachusetts Institute of Technology
%
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2, or (at your option)
%  any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software Foundation,
%  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
*/

#ifndef MEEP_VEC_H
#define MEEP_VEC_H

#include <complex>

using namespace std;

namespace meep {

const int NUM_FIELD_COMPONENTS = 15;
const int NUM_FIELD_TYPES = 4;

enum component { Ex=0, Ey, Er, Ep, Ez, Hx, Hy, Hr, Hp, Hz,
                 Dx, Dy, Dr, Dp, Dz, Dielectric };
enum derived_component { Sx=100, Sy, Sr, Sp, Sz, EnergyDensity,
			 D_EnergyDensity, H_EnergyDensity };
enum ndim { D1=0, D2, D3, Dcyl };
enum field_type { E_stuff=0, H_stuff=1, D_stuff=2, P_stuff=3 };
enum boundary_side { High=0, Low };
enum direction { X=0,Y,Z,R,P, NO_DIRECTION };
struct signed_direction {
  signed_direction(direction dd=X,bool f=false, complex<double> ph=1.0) {
    d = dd; flipped = f; phase = ph;
  };
  signed_direction(const signed_direction &sd) {
    d = sd.d; flipped = sd.flipped; phase = sd.phase;
  }
  signed_direction operator*(complex<double> ph);
  bool operator==(const signed_direction &sd) const { return (d == sd.d &&
						       flipped == sd.flipped
						       && phase == sd.phase); }
  bool operator!=(const signed_direction &sd) const { return !(*this == sd); }
  direction d;
  bool flipped;
  complex<double> phase;
};

inline int number_of_directions(ndim dim) {
  return (int) (dim + 1 - 2 * (dim == Dcyl));
}

inline direction start_at_direction(ndim dim) {
  return (direction) (((dim == D1) || (dim == Dcyl)) ? 2 : 0);
}

inline direction stop_at_direction(ndim dim) {
  return (direction) (dim + 1 + 2 * (dim == D1));
}

#define FOR_FIELD_TYPES(ft) for (field_type ft = E_stuff; \
                                 ft <= P_stuff; ft = (field_type) (ft+1))
#define FOR_ELECTRIC_COMPONENTS(c) for (component c = Ex; \
                                        c < Hx; c = (component) (c+1))
#define FOR_MAGNETIC_COMPONENTS(c) for (component c = Hz; \
                                        c > Ez; c = (component) (c-1))
#define FOR_D_COMPONENTS(c) for (component c = Dz; \
                                 c > Hz; c = (component) (c-1))
#define FOR_E_AND_D(e,d) for (component e = Ex, d = Dx; \
                              e <= Ez; e = (component) (e+1), d = (component) (d+1))
#define FOR_COMPONENTS(c) for (component c = Ex,loop_stop_co=Ey; \
                               c != loop_stop_co; \
                               c = (component)((c+1)%NUM_FIELD_COMPONENTS), \
                               loop_stop_co = Ex)
#define FOR_DIRECTIONS(d) for (direction d = X,loop_stop_di=Y; \
                               d != loop_stop_di; \
                               d = (direction)((d+1)%5), \
                               loop_stop_di = X)
#define FOR_SIDES(s) for (boundary_side s = High, loop_stop_bi=Low; \
			  s != loop_stop_bi; \
			  s = (boundary_side) ((s+1) % 2), \
			  loop_stop_bi = High)

#define LOOP_OVER_DIRECTIONS(dim, d) for (direction d = start_at_direction(dim), \
                                     loop_stop_directi = stop_at_direction(dim); \
                                     d < loop_stop_directi; d = (direction) (d+1))

// loop over indices idx from is to ie (inclusive) in v
#define LOOP_OVER_IVECS(v, is, ie, idx) \
  for (int loop_is1 = (is).yucky_val(0), \
           loop_is2 = (is).yucky_val(1), \
           loop_is3 = (is).yucky_val(2), \
           loop_n1 = ((ie).yucky_val(0) - loop_is1) / 2 + 1, \
           loop_n2 = ((ie).yucky_val(1) - loop_is2) / 2 + 1, \
           loop_n3 = ((ie).yucky_val(2) - loop_is3) / 2 + 1, \
           loop_d1 = (v).yucky_direction(0), \
           loop_d2 = (v).yucky_direction(1), \
           loop_d3 = (v).yucky_direction(2), \
           loop_s1 = (v).stride0((direction) loop_d1), \
           loop_s2 = (v).stride0((direction) loop_d2), \
           loop_s3 = (v).stride0((direction) loop_d3), \
           idx0 = (is - (v).little_corner()).yucky_val(0) / 2 * loop_s1 \
                + (is - (v).little_corner()).yucky_val(1) / 2 * loop_s2 \
                + (is - (v).little_corner()).yucky_val(2) / 2 * loop_s3,\
           loop_i1 = 0; loop_i1 < loop_n1; loop_i1++) \
    for (int loop_i2 = 0; loop_i2 < loop_n2; loop_i2++) \
      for (int idx = idx0 + loop_i1*loop_s1 + loop_i2*loop_s2, \
           loop_i3 = 0; loop_i3 < loop_n3; loop_i3++, idx+=loop_s3)

#define LOOP_OVER_VOL(v, c, idx) \
  LOOP_OVER_IVECS(v, (v).little_corner() + (v).iyee_shift(c), (v).big_corner() + (v).iyee_shift(c), idx)

#define LOOP_OVER_VOL_OWNED(v, c, idx) \
  LOOP_OVER_IVECS(v, (v).little_owned_corner(c), (v).big_corner(), idx)

#define LOOP_OVER_VOL_NOTOWNED(v, c, idx) \
 for (ivec loop_notowned_is((v).dim,0), loop_notowned_ie((v).dim,0); \
      loop_notowned_is == zero_ivec((v).dim);) \
   for (int loop_ibound = 0; (v).get_boundary_icorners(c, loop_ibound,     \
		  				       &loop_notowned_is,  \
						       &loop_notowned_ie); \
	loop_ibound++) \
     LOOP_OVER_IVECS(v, loop_notowned_is, loop_notowned_ie, idx)

#define IVEC_LOOP_AT_BOUNDARY 					\
 ((loop_s1 != 0 && (loop_i1 == 0 || loop_i1 == loop_n1-1)) ||	\
  (loop_s2 != 0 && (loop_i2 == 0 || loop_i2 == loop_n2-1)) ||	\
  (loop_s3 != 0 && (loop_i3 == 0 || loop_i3 == loop_n3-1)))

#define IVEC_LOOP_ILOC(v, iloc) \
  ivec iloc((v).dim); \
  iloc.set_direction(direction(loop_d1), loop_is1 + 2*loop_i1); \
  iloc.set_direction(direction(loop_d2), loop_is2 + 2*loop_i2); \
  iloc.set_direction(direction(loop_d3), loop_is3 + 2*loop_i3)

#define IVEC_LOOP_LOC(v, loc) \
  vec loc((v).dim); \
  loc.set_direction(direction(loop_d1), (0.5*loop_is1 + loop_i1) * (v).inva); \
  loc.set_direction(direction(loop_d2), (0.5*loop_is2 + loop_i2) * (v).inva); \
  loc.set_direction(direction(loop_d3), (0.5*loop_is3 + loop_i3) * (v).inva)

// integration weight for using LOOP_OVER_IVECS with field::integrate
#define IVEC_LOOP_WEIGHT1x(s0, s1, e0, e1, i, n, dir) ((i > 1 && i < n - 2) ? 1.0 : (i == 0 ? (s0).in_direction(direction(dir)) : (i == 1 ? (s1).in_direction(direction(dir)) : i == n - 1 ? (e0).in_direction(direction(dir)) : (i == n - 2 ? (e1).in_direction(direction(dir)) : 1.0))))
#define IVEC_LOOP_WEIGHT1(s0, s1, e0, e1, k) IVEC_LOOP_WEIGHT1x(s0, s1, e0, e1, loop_i##k,loop_n##k,loop_d##k)
#define IVEC_LOOP_WEIGHT(s0, s1, e0, e1, dV) (IVEC_LOOP_WEIGHT1(s0, s1, e0, e1, 3) * (IVEC_LOOP_WEIGHT1(s0, s1, e0, e1, 2) * ((dV) * IVEC_LOOP_WEIGHT1(s0, s1, e0, e1, 1))))

inline signed_direction flip(signed_direction d) {
  signed_direction d2 = d;
  d2.flipped = !d.flipped;
  return d2;
}

inline bool has_direction(ndim dim, direction d) {
  LOOP_OVER_DIRECTIONS(dim, dd) if (dd == d) return true;
  return false;
}

// true if d is polar while dim is cartesian, or vice versa 
inline bool coordinate_mismatch(ndim dim, direction d) {
  return (d != NO_DIRECTION &&
	  ((dim >= D1 && dim <= D3 && d != X && d != Y && d != Z) ||
	   (dim == Dcyl && d != R && d != P && d != Z)));
}

extern void abort(const char *, ...); // mympi.cpp

inline bool is_electric(component c) { return c < Hx; }
inline bool is_magnetic(component c) { return c >= Hx && c < Dx; }
inline bool is_D(component c) { return c >= Dx && c < Dielectric; }
inline bool is_derived(int c) { return c >= Sx; }
inline bool is_poynting(derived_component c) { return c < EnergyDensity; }
inline bool is_energydensity(derived_component c) { return c>=EnergyDensity; }
inline field_type type(component c) {
  if (is_electric(c)) return E_stuff;
  else if (is_magnetic(c)) return H_stuff;
  else if (is_D(c)) return D_stuff;
  abort("Invalid field in type.\n");
  return E_stuff; // This is never reached.
}
const char *component_name(component c);
const char *component_name(derived_component c);
const char *component_name(int c);
const char *direction_name(direction);
const char *dimension_name(ndim);

direction component_direction(int c);
int direction_component(int c, direction d);
inline direction component_direction(component c) {
  switch (c) {
  case Ex: case Hx: case Dx: return X;
  case Ey: case Hy: case Dy: return Y;
  case Ez: case Hz: case Dz: return Z;
  case Er: case Hr: case Dr: return R;
  case Ep: case Hp: case Dp: return P;
  case Dielectric: return NO_DIRECTION;
  }
  return X; // This code is never reached...
}
inline direction component_direction(derived_component c) {
  switch (c) {
  case Sx: return X;
  case Sy: return Y;
  case Sz: return Z;
  case Sr: return R;
  case Sp: return P;
  case EnergyDensity: case D_EnergyDensity: case H_EnergyDensity:
    return NO_DIRECTION;
  }
  return X; // This code is never reached...
}
inline direction component_direction(int c) {
  if (is_derived(c))
    return component_direction(derived_component(c));
  else
    return component_direction(component(c));
}
inline component direction_component(component c, direction d) {
  component start_point;
  if (is_electric(c)) start_point = Ex;
  else if (is_magnetic(c)) start_point = Hx;
  else if (is_D(c)) start_point = Dx;
  else if (c == Dielectric && d == NO_DIRECTION) return Dielectric;
  else abort("unknown field component %d", c);
  switch (d) {
  case X: return start_point;
  case Y: return (component) (start_point + 1);
  case Z: return (component) (start_point + 4);
  case R: return (component) (start_point + 2);
  case P: return (component) (start_point + 3);
  case NO_DIRECTION: abort("vector %d component in NO_DIRECTION", c);
  }
  return Ex; // This is never reached.
}
inline derived_component direction_component(derived_component c, direction d) {
  derived_component start_point;
  if (is_poynting(c)) start_point = Sx;
  else if (is_energydensity(c) && d == NO_DIRECTION) return c;
  else abort("unknown field component %d", c);
  switch (d) {
  case X: return start_point;
  case Y: return (derived_component) (start_point + 1);
  case Z: return (derived_component) (start_point + 4);
  case R: return (derived_component) (start_point + 2);
  case P: return (derived_component) (start_point + 3);
  case NO_DIRECTION: abort("vector %d derived_component in NO_DIRECTION", c);
  }
  return Sx; // This is never reached.
}
inline int direction_component(int c, direction d) {
  if (is_derived(c))
    return int(direction_component(derived_component(c), d));
  else
    return int(direction_component(component(c), d));
}

inline bool coordinate_mismatch(ndim dim, component c) {
  return coordinate_mismatch(dim, component_direction(c));
}
inline bool coordinate_mismatch(ndim dim, derived_component c) {
  return coordinate_mismatch(dim, component_direction(c));
}

class vec;
vec veccyl(double rr, double zz);
vec zero_vec(ndim);

class vec {
 public:
  vec() {};
  vec(ndim di) { 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 xx, double yy, double zz) {
    dim = D3; t[X] = xx; t[Y] = yy; t[Z] = zz; };
  friend vec veccyl(double rr, double zz);
  ~vec() {};

  vec operator+(const vec &a) const {
    vec result = a;
    LOOP_OVER_DIRECTIONS(dim, d) result.t[d] += t[d];
    return result;
  };

  vec operator+=(const vec &a) {
    LOOP_OVER_DIRECTIONS(dim, d) t[d] += a.t[d];
    return *this;
  };

  vec operator-(const vec &a) const {
    vec result = a;
    LOOP_OVER_DIRECTIONS(dim, d) result.t[d] = t[d] - result.t[d];
    return result;
  };

  vec operator-(void) const {
    vec result(dim);
    LOOP_OVER_DIRECTIONS(dim, d) result.t[d] = -t[d];
    return result;
  };

  vec operator-=(const vec &a) {
    LOOP_OVER_DIRECTIONS(dim, d) t[d] -= a.t[d];
    return *this;
  };

  bool operator!=(const vec &a) const {
    LOOP_OVER_DIRECTIONS(dim, d) if (t[d] != a.t[d]) return true;
    return false;
  };

  bool operator==(const vec &a) const {
    LOOP_OVER_DIRECTIONS(dim, d) if (t[d] != a.t[d]) return false;
    return true;
  };

  vec round_float(void) const {
    vec result = *this;
    LOOP_OVER_DIRECTIONS(dim, d) result.t[d] = float(result.t[d]);
    return result;
  }

  vec operator*(double s) const {
    vec result = *this;
    LOOP_OVER_DIRECTIONS(dim, d) result.t[d] *= s;
    return result;
  };

  vec operator/(double s) const {
    vec result = *this;
    LOOP_OVER_DIRECTIONS(dim, d) result.t[d] *= (1.0/s);
    return result;
  };

  // I use & as a dot product.
  double operator&(const vec &a) const {
    double result = 0.0;
    LOOP_OVER_DIRECTIONS(dim, d) result += t[d] * a.t[d];
    return result;
  };
  ndim dim;

  double r() const { return t[R]; };
  double x() const { return t[X]; };
  double y() const { return t[Y]; };
  double z() const { return t[Z]; };
  double in_direction(direction d) const { return t[d]; };
  void set_direction(direction d, double val) { t[d] = val; };

  double project_to_boundary(direction, double boundary_loc);
  friend vec zero_vec(ndim);
 private:
  double t[5];
};

inline double abs(const vec &v) { return sqrt(v & v); }

inline vec zero_vec(ndim di) {
  vec v(di); LOOP_OVER_DIRECTIONS(di, d) v.set_direction(d, 0.0);
  return v;
}

inline vec unit_vec(ndim di, direction d) {
  vec v(zero_vec(di));
  v.set_direction(d, 1.0);
  return v;
}

inline vec clean_vec(const vec &v, double val_unused = 0.0) {
  vec vc(v.dim, val_unused);
  LOOP_OVER_DIRECTIONS(v.dim, d) vc.set_direction(d, v.in_direction(d));
  return vc;
}

inline vec veccyl(double rr, double zz) {
  vec v(Dcyl); v.t[R] = rr; v.t[Z] = zz; return v;
}

class ivec;
ivec iveccyl(int xx, int yy);
ivec zero_ivec(ndim);
ivec one_ivec(ndim);

class ivec {
 public:
  ivec() { dim = D2; t[X] = t[Y] = 0; };
  ivec(ndim di) { 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[Z] = zz; };
  ivec(int xx, int yy) { dim = D2; t[X] = xx; t[Y] = yy; };
  ivec(int xx, int yy, int zz) {
    dim = D3; t[X] = xx; t[Y] = yy; t[Z] = zz; };
  friend ivec iveccyl(int xx, int yy);