vec.cpp

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

/* 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.
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex>

#include "meep_internals.hpp"

namespace meep {

ivec volume::round_vec(const vec &p) const {
  ivec result(dim);
  LOOP_OVER_DIRECTIONS(dim, d)
    result.set_direction(d, my_round(p.in_direction(d) * 2 * a));
  return result;
}

void volume::set_origin(const ivec &o) {
  io = o;
  origin = operator[](io); // adjust origin to match io
}

void volume::set_origin(direction d, int o) {
  io.set_direction(d, o);
  origin = operator[](io); // adjust origin to match io
}

void volume::set_origin(const vec &o) {
  set_origin(round_vec(o));
}

const char *dimension_name(ndim dim) {
  switch (dim) {
  case D1: return "1D";
  case D2: return "2D";
  case D3: return "3D";
  case Dcyl: return "Cylindrical";
  }
  return "Error in dimension_name";
}

const char *direction_name(direction d) {
  switch (d) {
  case X: return "x";
  case Y: return "y";
  case Z: return "z";
  case R: return "r";
  case P: return "phi";
  case NO_DIRECTION: return "no_direction";
  }
  return "Error in direction_name";
}

const char *component_name(component c) {
  if (is_derived(int(c))) return component_name(derived_component(c));
  switch (c) {
  case Er: return "er";
  case Ep: return "ep";
  case Ez: return "ez";
  case Hr: return "hr";
  case Hp: return "hp";
  case Hz: return "hz";
  case Ex: return "ex";
  case Ey: return "ey";
  case Hx: return "hx";
  case Hy: return "hy";
  case Dx: return "dx";
  case Dy: return "dy";
  case Dz: return "dz";
  case Dr: return "dr";
  case Dp: return "dp";
  case Dielectric: return "eps";
  }
  return "Error in component_name";
}

const char *component_name(derived_component c) {
  if (!is_derived(int(c))) return component_name(component(c));
  switch (c) {
  case Sr: return "sr";
  case Sp: return "sp";
  case Sz: return "sz";
  case Sx: return "sx";
  case Sy: return "sy";
  case EnergyDensity: return "energy";
  case D_EnergyDensity: return "denergy";
  case H_EnergyDensity: return "henergy";
  }
  return "Error in component_name";
}

const char *component_name(int c) {
  return (is_derived(c) ? component_name(derived_component(c))
	  : component_name(component(c)));
}

vec min(const vec &vec1, const vec &vec2) {
  vec m(vec1.dim);
  LOOP_OVER_DIRECTIONS(vec1.dim, d)
    m.set_direction(d, min(vec1.in_direction(d), vec2.in_direction(d)));
  return m;
}

vec max(const vec &vec1, const vec &vec2) {
  vec m(vec1.dim);
  LOOP_OVER_DIRECTIONS(vec1.dim, d)
    m.set_direction(d, max(vec1.in_direction(d), vec2.in_direction(d)));
  return m;
}

ivec min(const ivec &ivec1, const ivec &ivec2) {
  ivec m(ivec1.dim);
  LOOP_OVER_DIRECTIONS(ivec1.dim, d)
    m.set_direction(d, min(ivec1.in_direction(d), ivec2.in_direction(d)));
  return m;
}

ivec max(const ivec &ivec1, const ivec &ivec2) {
  ivec m(ivec1.dim);
  LOOP_OVER_DIRECTIONS(ivec1.dim, d)
    m.set_direction(d, max(ivec1.in_direction(d), ivec2.in_direction(d)));
  return m;
}

geometric_volume::geometric_volume(const vec &vec1, const vec &vec2) {
  min_corner = min(vec1, vec2);
  max_corner = max(vec1, vec2);
  dim = vec1.dim; 
}

geometric_volume::geometric_volume(const vec &pt) {
  dim = pt.dim; 
  min_corner = pt;
  max_corner = pt;
}

double geometric_volume::computational_volume() const {
  double vol = 1.0; 
  LOOP_OVER_DIRECTIONS(dim,d) vol *= in_direction(d);
  return vol;
}

double geometric_volume::integral_volume() const {
  double vol = 1.0; 
  LOOP_OVER_DIRECTIONS(dim, d)
    if (in_direction(d) != 0.0) vol *= in_direction(d);
  if (dim == Dcyl) vol *= pi * (in_direction_max(R) + in_direction_min(R));
  return vol;
}

double geometric_volume::full_volume() const {
  double vol = computational_volume(); 
  if (dim == Dcyl) vol *= pi * (in_direction_max(R) + in_direction_min(R));
  return vol;
}

double geometric_volume::diameter() const {
  double diam = 0.0;
  LOOP_OVER_DIRECTIONS(dim,d) {
    diam = max(diam, in_direction(d));
  }
  return diam;
}

geometric_volume geometric_volume::intersect_with(const geometric_volume &a) const {
  if (a.dim != dim) abort("Can't intersect volumes of dissimilar dimensions.\n");
  geometric_volume result(dim);
  LOOP_OVER_DIRECTIONS(dim, d) {
    double minval = max(in_direction_min(d), a.in_direction_min(d));
    double maxval = min(in_direction_max(d), a.in_direction_max(d));
    if (minval > maxval)
      return geometric_volume(zero_vec(dim), zero_vec(dim));
    result.set_direction_min(d, minval);
    result.set_direction_max(d, maxval);
  }
  return result;
}

bool geometric_volume::intersects(const geometric_volume &a) const {
  if (a.dim != dim) abort("Can't intersect volumes of dissimilar dimensions.\n");
  LOOP_OVER_DIRECTIONS(dim, d) {
    double minval = max(in_direction_min(d), a.in_direction_min(d));
    double maxval = min(in_direction_max(d), a.in_direction_max(d));
    if (minval > maxval)
      return false;
  }
  return true;
}

// Return normal direction to volume, if the volume is dim-1 dimensional;
// otherwise, return NO_DIRECTION.
direction geometric_volume::normal_direction() const {
  direction d = NO_DIRECTION;
  switch (dim) {
  case D1: d = Z; break;
  case D2:
    if (in_direction(X) == 0 && in_direction(Y) > 0)
      d = X;
    else if (in_direction(X) > 0 && in_direction(Y) == 0)
      d = Y;
    break;
  case Dcyl:
    if (in_direction(R) == 0 && in_direction(Z) > 0)
      d = R;
    else if (in_direction(R) > 0 && in_direction(Z) == 0)
      d = Z;
    break;
  case D3: {
    bool zx = in_direction(X) == 0;
    bool zy = in_direction(Y) == 0;
    bool zz = in_direction(Z) == 0;
    if (zx && !zy && !zz) d = X;
    else if (!zx && zy && !zz) d = Y;
    else if (!zx && !zy && zz) d = Z;
    break;
  }
  }
  return d;
}

/* Used for n=0,1,2 nested loops in macros.  We should arrange
   the ordering so that this gives most efficient traversal of
   a field array, where n=2 is the innermost loop. */
static direction yucky_dir(ndim dim, int n) {
  if (dim == Dcyl)
    switch (n) {
    case 0: return P;
    case 1: return R;
    case 2: return Z;
    }
  else if (dim == D2)
    return (direction) ((n + 2) % 3); /* n = 0,1,2 gives Z, X, Y */
  return (direction) n ;
}

int ivec::yucky_val(int n) const {
  if (has_direction(dim, yucky_dir(dim, n)))
    return in_direction(yucky_dir(dim, n));
  return 0;
}

int volume::yucky_num(int n) const {
  if (has_direction(dim, yucky_dir(dim, n)))
    return num_direction(yucky_dir(dim, n));
  return 1;
}

direction volume::yucky_direction(int n) const {
  return yucky_dir(dim, n);
}

geometric_volume volume::surroundings() const {
  return geometric_volume(operator[](little_corner()), 
			  operator[](big_corner()));
}

geometric_volume volume::interior() const {
  return geometric_volume(operator[](little_corner()), 
			  operator[](big_corner() - one_ivec(dim) * 2));
}

void volume::update_ntot() {
  the_ntot = 1;
  LOOP_OVER_DIRECTIONS(dim, d) the_ntot *= num[d%3] + 1;
}

void volume::set_num_direction(direction d, int value) {
  num[d%3] = value; num_changed();
}

volume::volume(ndim td, double ta, int na, int nb, int nc) {
  dim = td; a = ta; inva = 1.0 / ta;
  num[0] = na;
  num[1] = nb;
  num[2] = nc;
  num_changed();
  set_origin(zero_vec(dim));
}

component volume::eps_component() const {
  switch (dim) {
  case D1: return Hy;
  case D2: return Hz;
  case D3: return Dielectric;
  case Dcyl: return Hp;
  }
  abort("Unsupported dimensionality eps.\n");
  return Ex;
}

vec volume::yee_shift(component c) const {
  return operator[](iyee_shift(c));
}

/* Return array offsets to average with a given array location of c in
   order to get c on the "dielectric" grid.  Then, to get the
   dielectric grid point i, you should average c over the four
   locations: i, i+offset1, i+offset2, i+offset1+offset2. 
   (offset2, and possibly offset1, may be zero if only 2 or 1
   locations need to be averaged). */
void volume::yee2diel_offsets(component c, int &offset1, int &offset2) {
  offset1 = offset2 = 0;
  LOOP_OVER_DIRECTIONS(dim,d) {
    if (!iyee_shift(c).in_direction(d)) {
      if (offset2) 
	abort("weird yee shift for component %s", component_name(c));
      if (offset1) offset2 = stride(d);
      else offset1 = stride(d);
    }
  }
}

bool geometric_volume::contains(const vec &p) const {
  LOOP_OVER_DIRECTIONS(dim,d) {
    if (p.in_direction(d) > in_direction_max(d) ||
        p.in_direction(d) < in_direction_min(d)) return false;
  }
  return true;
}

bool volume::contains(const ivec &p) const {
  // containts returns true if the volume has information about this grid
  // point.
  const ivec o = p - io;
  LOOP_OVER_DIRECTIONS(dim, d)
    if (o.in_direction(d) < 0 || o.in_direction(d) >= (num_direction(d)+1)*2)
      return false;
  return true;
}

bool volume::contains(const vec &p) const {
  // containts returns true if the volume has any information in it
  // relevant to the point p.  Basically has is like owns (see below)
  // except it is more lenient, in that more than one lattice may contain a
  // given point.
  const vec o = p - origin;
  LOOP_OVER_DIRECTIONS(dim, d)
    if (o.in_direction(d) < -inva || o.in_direction(d) > num_direction(d)*inva+inva)
      return false;
  return true;
}

/* Compute the corners (cs,ce) of the ib-th boundary for component c,
   returning true if ib is a valid index (ib = 0..#boundaries-1).  The
   boundaries are all the points that are in but not owned by the
   volume, and are a set of *disjoint* regions.  The main purpose of
   this function is currently to support the LOOP_OVER_NOT_OWNED
   macro.  (In the future, it may be used for other
   boundary-element-type computations, too.) */
bool volume::get_boundary_icorners(component c, int ib,
				   ivec *cs, ivec *ce) const {
  ivec cl(little_corner() + iyee_shift(c));
  ivec cb(big_corner() + iyee_shift(c));
  ivec clo(little_owned_corner(c));
  ivec cbo(big_corner() - iyee_shift(c));
  *cs = cl;
  *ce = cb;
  bool ib_found = false;
  int jb = 0;
  LOOP_OVER_DIRECTIONS(dim, d) {
    if (cl.in_direction(d) < clo.in_direction(d)) {
      if (jb == ib) {
	ce->set_direction(d, cs->in_direction(d));
	ib_found = true;
	break;
      }
      cs->set_direction(d, clo.in_direction(d));
      jb++;
    }
    if (cb.in_direction(d) > cbo.in_direction(d)) {
      if (jb == ib) {
	cs->set_direction(d, ce->in_direction(d));
	ib_found = true;
	break;
      }
      ce->set_direction(d, cbo.in_direction(d));
      jb++;
    }
  }
  if (!ib_found) { // yucky interaction here with LOOP_OVER_VOL_NOTOWNED
    *cs = one_ivec(dim);
    *ce = -one_ivec(dim);
  }
  return ib_found;
}

// first "owned" point for c in volume (see also volume::owns)
ivec volume::little_owned_corner(component c) const {
  ivec iloc(little_corner() + one_ivec(dim)*2 - iyee_shift(c));
  if (dim == Dcyl && origin.r() == 0.0 && iloc.r() == 2)
    iloc.set_direction(R, 0);
  return iloc;
}

int volume::nowned(component c) const {
  int n = 1;
  ivec v = big_corner() - little_owned_corner(c);
  LOOP_OVER_DIRECTIONS(dim, d) n *= v.in_direction(d) / 2 + 1;
  return n;
}

bool volume::owns(const ivec &p) const {
  // owns returns true if the point "owned" by this volume, meaning that it
  // is the volume that would timestep the point.
  const ivec o = p - io;
  if (dim == Dcyl) {
    if (origin.r() == 0.0 && o.z() > 0 && o.z() <= nz()*2 &&
        o.r() == 0) return true;
    return o.r() > 0 && o.z() > 0 &&
           o.r() <= nr()*2 && o.z() <= nz()*2;
  } else if (dim == D3) {
    return
      o.x() > 0 && o.x() <= nx()*2 &&
      o.y() > 0 && o.y() <= ny()*2 &&
      o.z() > 0 && o.z() <= nz()*2;
  } else if (dim == D2) {
    return
      o.x() > 0 && o.x() <= nx()*2 &&
      o.y() > 0 && o.y() <= ny()*2;
  } else if (dim == D1) {
    return o.z() > 0 && o.z() <= nz()*2;
  } else {
    abort("Unsupported dimension in owns.\n");
    return false;
  }
}

int volume::has_boundary(boundary_side b,direction d) const {
  switch (dim) {
  case Dcyl: return d == Z || (d == R && (b == High || get_origin().r() > 0));
  case D1: return d == Z;
  case D2: return d == X || d == Y;
  case D3: return d == X || d == Y || d == Z;
  }
  return 0; // This should never be reached.
}

int volume::index(component c, const ivec &p) const {
  const ivec offset = p - io - iyee_shift(c);
  int idx = 0;
  LOOP_OVER_DIRECTIONS(dim,d) idx += offset.in_direction(d)/2*stride(d);
  return idx;
}

void volume::set_strides() {
  FOR_DIRECTIONS(d) the_stride0[d] = 0; // Yuck yuck yuck.
  LOOP_OVER_DIRECTIONS(dim,d)
    switch(d) {
    case Z: the_stride0[d] = 1; break;
    case R: the_stride0[d] = nz()+1; break;
    case X: the_stride0[d] = (nz()+1)*(ny() + 1); break;
    case Y: the_stride0[d] = nz() + 1; break;
    case P: break; // There is no phi stride...
    case NO_DIRECTION: break; // no stride here, either
    }
  FOR_DIRECTIONS(d) the_stride[d] = the_stride0[d] ? the_stride0[d] : 1;
}

static inline void stupidsort(int *ind, double *w, int l) {
  while (l) {
    if (fabs(w[0]) < 2e-15) {
      w[0] = w[l-1];
      ind[0] = ind[l-1];
      w[l-1] = 0.0;
      ind[l-1] = 0;
    } else {