vec.hpp

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

  ~ivec() {};

  // Only an idiot (or a macro) would use a yucky function.  Don't be an
  // idiot.
  int yucky_val(int) const;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  ivec round_up_to_even(void) const { 
    ivec result(dim);
    LOOP_OVER_DIRECTIONS(dim, d) 
      result.t[d] = t[d] + (t[d] >= 0 ? t[d] : -t[d]) % 2; 
    return result;
  }

  friend ivec zero_ivec(ndim);
  friend ivec one_ivec(ndim);
 private:
  int t[5];
};

inline ivec zero_ivec(ndim di) {
  ivec v; v.dim = di; LOOP_OVER_DIRECTIONS(di, d) v.set_direction(d, 0);
  return v;
}

inline ivec one_ivec(ndim di) {
  ivec v; v.dim = di; LOOP_OVER_DIRECTIONS(di, d) v.set_direction(d, 1);
  return v;
}

inline ivec unit_ivec(ndim di, direction d) {
  ivec v(zero_ivec(di));
  v.set_direction(d, 1);
  return v;
}

inline ivec iveccyl(int rr, int zz) {
  ivec v(Dcyl); v.t[R] = rr; v.t[Z] = zz; return v;
}

vec max(const vec &vec1, const vec &vec2);
vec min(const vec &vec1, const vec &vec2);
ivec max(const ivec &ivec1, const ivec &ivec2);
ivec min(const ivec &ivec1, const ivec &ivec2);
ivec max_to_all(const ivec &v); // in mympi.cpp

class geometric_volume {
 public:
  ndim dim;
  geometric_volume(ndim di) { dim = di; min_corner.dim = di; max_corner.dim = di; };
  geometric_volume(const vec &vec1, const vec &vec2);
  geometric_volume(const vec &pt);
  void set_direction_min(direction d, double val) { min_corner.set_direction(d, val); };
  void set_direction_max(direction d, double val) { max_corner.set_direction(d, val); };
  double in_direction_min(direction d) const { return min_corner.in_direction(d); };
  double in_direction_max(direction d) const { return max_corner.in_direction(d); };
  double in_direction(direction d) const { return in_direction_max(d) - in_direction_min(d); }
  double computational_volume() const; 
  double integral_volume() const;
  double full_volume() const;
  vec center() const { return (min_corner + max_corner) * 0.5; }
  double diameter() const;
  bool contains(const vec &h) const;
  geometric_volume intersect_with(const geometric_volume &a) const;
  geometric_volume operator&(const geometric_volume &a) const {
    return intersect_with(a);
  };
  geometric_volume operator|(const geometric_volume &a) const {
    return geometric_volume(min(min_corner, a.min_corner),
			    max(max_corner, a.max_corner));
  };
  geometric_volume operator+(const vec &a) const {
    return geometric_volume(min_corner + a, max_corner + a);
  }
  geometric_volume operator+=(const vec &a) {
    min_corner += a; max_corner += a;
    return *this;
  }
  geometric_volume operator-(const vec &a) const {
    return geometric_volume(min_corner - a, max_corner - a);
  }
  geometric_volume operator-=(const vec &a) {
    min_corner -= a; max_corner -= a;
    return *this;
  }
  bool operator==(const geometric_volume &a) const {
    return (min_corner == a.min_corner && max_corner == a.max_corner);
  }
  bool operator!=(const geometric_volume &a) const { return !(*this == a); };
  geometric_volume round_float(void) const {
    return geometric_volume(min_corner.round_float(),max_corner.round_float());
  }
  bool intersects(const geometric_volume &a) const;
  bool operator&&(const geometric_volume &a) const {
    return intersects(a);
  };
  vec get_min_corner() const { return min_corner; };
  vec get_max_corner() const { return max_corner; };
  direction normal_direction() const;
 private:
  vec min_corner, max_corner;
};

class volume;
volume volcyl(double rsize, double zsize, double a);
volume volone(double zsize, double a);
volume vol1d(double zsize, double a);
volume voltwo(double xsize, double ysize, double a);
volume vol2d(double xsize, double ysize, double a);
volume vol3d(double xsize, double ysize, double zsize, double a);

class volume {
 public:
  volume() {};

  ndim dim;
  double a, inva /* = 1/a */;

  void print() const;
  int stride(direction d) const { return the_stride[d]; };
  int stride0(direction d) const { return the_stride0[d]; };
  int num_direction(direction d) const {
    return num[((int) d) % 3];
  };
  // Only an idiot (or a macro) would use a yucky function.  Don't be an
  // idiot.
  int yucky_num(int) const;
  direction yucky_direction(int) const;
  void set_num_direction(direction d, int value);
  int nr() const { return num_direction(R); }
  int nx() const { return num_direction(X); }
  int ny() const { return num_direction(Y); }
  int nz() const { return num_direction(Z); }

  bool has_field(component c) const {
    if (dim == D1) return c == Ex || c == Hy || c == Dx;
    return (dim == Dcyl)?component_direction(c)>Y:component_direction(c)<R;
  }
  int has_boundary(boundary_side,direction) const;

  vec dr() const;
  vec dx() const;
  vec dy() const;
  vec dz() const;

  int ntot() const { return the_ntot; }
  int nowned_min() const { int n = 1; LOOP_OVER_DIRECTIONS(dim,d) n *= num_direction(d); return n; }
  int nowned(component c) const;
  vec operator[](const ivec &p) const { return p*(0.5*inva); };
  int index(component, const ivec &) const;
  ivec round_vec(const vec &) const;
  void interpolate(component, const vec &, int indices[8], double weights[8]) const;
  void interpolate(component, const vec &, ivec locs[8], double weights[8]) const;

  geometric_volume dV(component c, int index) const;
  geometric_volume dV(const ivec &, double diameter = 1.0) const;
  bool intersect_with(const volume &vol_in, volume *intersection = NULL, volume *others = NULL, int *num_others = NULL) const;
  double rmin() const;
  double rmax() const;
  double xmin() const;
  double xmax() const;
  double ymin() const;
  double ymax() const;
  double zmin() const;
  double zmax() const;
  vec center() const;
  ivec icenter() const;
  vec loc(component, int index) const;
  vec loc_at_resolution(int index, double res) const;
  int ntot_at_resolution(double res) const;
  ivec iloc(component, int index) const;

  int yee_index(component c) const {
    int idx = 0;
    LOOP_OVER_DIRECTIONS(dim,d)
      idx += (1-iyee_shift(c).in_direction(d))*stride(d);
    return idx;
  }
  vec yee_shift(component) const;
  component eps_component() const;
  void yee2diel_offsets(component c, int &offset1, int &offset2);

  double boundary_location(boundary_side, direction) const;
  ivec big_corner() const;
  ivec little_corner() const { return io; };
  vec corner(boundary_side b) const;

  bool contains(const vec &) const;
  bool contains(const ivec &) const;
  ivec little_owned_corner(component c) const;
  bool owns(const ivec &) const;
  geometric_volume surroundings() const;
  geometric_volume interior() const;

  bool get_boundary_icorners(component c, int ib, ivec *cs, ivec *ce) const;

  friend volume volcyl(double rsize, double zsize, double a);
  friend volume volone(double zsize, double a);
  friend volume vol1d(double zsize, double a);
  friend volume voltwo(double xsize, double ysize, double a);
  friend volume vol2d(double xsize, double ysize, double a);
  friend volume vol3d(double xsize, double ysize, double zsize, double a);

  volume split(int num, int which) const;
  volume split_by_effort(int num, int which, int Ngv = 0, const volume *gv = NULL, double *effort = NULL) const;
  volume split_at_fraction(bool want_high, int numer) const;
  volume halve(direction d) const;
  void pad_self(direction d);
  volume pad(direction d) const;
  volume pad() const {
       volume v(*this);
       LOOP_OVER_DIRECTIONS(dim,d)
	    v.pad_self(d);
       return v;
  }
  ivec iyee_shift(component c) const {
    ivec out = zero_ivec(dim);
    LOOP_OVER_DIRECTIONS(dim,d)
      if (c == Dielectric ||
          ((is_electric(c) || is_D(c)) && d == component_direction(c)) ||
          (is_magnetic(c) && d != component_direction(c)))
        out.set_direction(d,1);
    return out;
  }

  vec get_origin() const { return origin; }
  void set_origin(const vec &o);
  void set_origin(const ivec &o);
  void shift_origin(const vec &s) { set_origin(origin + s); }
  void shift_origin(const ivec &s) { set_origin(io + s); }
  void shift_origin(direction d, int s) {shift_origin(unit_ivec(dim, d) * s);}
  void set_origin(direction d, int o);
  void center_origin(void) { shift_origin(-icenter()); }
  double origin_in_direction(direction d) const{return origin.in_direction(d);}
  double origin_r() const { return origin.r(); }
  double origin_x() const { return origin.x(); }
  double origin_y() const { return origin.y(); }
  double origin_z() const { return origin.z(); }

 private:
  volume(ndim d, double ta, int na, int nb, int nc);
  ivec io; // integer origin ... always change via set_origin etc.!
  vec origin; // cache of operator[](io), for performance
  void update_ntot();
  void set_strides();
  void num_changed() { update_ntot(); set_strides(); }
  int num[3];
  int the_stride0[5], the_stride[5];
  int the_ntot;
};

class geometric_volume_list;

class symmetry;
symmetry identity();
symmetry rotate4(direction,const volume &);
symmetry rotate2(direction,const volume &);
symmetry mirror(direction,const volume &);
symmetry r_to_minus_r_symmetry(double m);

class symmetry {
 public:
  symmetry();
  symmetry(const symmetry &);
  ~symmetry();
  friend symmetry identity();
  friend symmetry rotate4(direction,const volume &);
  friend symmetry rotate2(direction,const volume &);
  friend symmetry mirror(direction,const volume &);

  signed_direction transform(direction d, int n) const;
  ivec transform(const ivec &, int n) const;
  vec transform(const vec &, int n) const;
  ivec transform_unshifted(const ivec &, int n) const;
  geometric_volume transform(const geometric_volume &, int n) const;
  component transform(component, int n) const;
  complex<double> phase_shift(component, int n) const;
  derived_component transform(derived_component, int n) const;
  complex<double> phase_shift(derived_component, int n) const;
  int transform(int, int n) const;
  complex<double> phase_shift(int, int n) const;
  int multiplicity() const;
  bool is_primitive(const ivec &) const;

  geometric_volume_list *reduce(const geometric_volume_list *gl) const;

  symmetry operator+(const symmetry &) const;
  symmetry operator*(complex<double>) const;
  symmetry operator-(const symmetry &b) const { return *this + b * (-1.0); }
  symmetry operator-(void) const { return *this * (-1.0); }
  void operator=(const symmetry &);
  bool operator==(const symmetry &) const;
  bool operator!=(const symmetry &S) const { return !(*this == S); };
 private:
  signed_direction S[5];
  complex<double> ph;
  vec symmetry_point;
  ivec i_symmetry_point;
  int g; // g is the multiplicity of the symmetry.
  symmetry *next;
  friend symmetry r_to_minus_r_symmetry(double m);
};

symmetry identity();
symmetry rotate4(direction,const volume &);
symmetry rotate2(direction,const volume &);
symmetry mirror(direction,const volume &);

class geometric_volume_list {
public:
  geometric_volume_list(const geometric_volume &gv, int c, complex<double> weight = 1.0, geometric_volume_list *next = 0) : gv(gv), c(c), weight(weight), next(next) {}
  ~geometric_volume_list() { delete next; }
  
  geometric_volume gv;
  int c; // component or derived component associated with gv (e.g. for flux)
  complex<double> weight;
  geometric_volume_list *next;
};

} /* namespace meep */

#endif /* MEEP_VEC_H */