vec.cpp

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

      w += 1;
      ind += 1;
    }
    l -= 1;
  }
}

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

void volume::interpolate(component c, const vec &p,
                         int indices[8], double weights[8]) const {
  ivec locs[8];
  interpolate(c, p, locs, weights);
  for (int i=0;i<8&&weights[i];i++)
    if (!owns(locs[i])) weights[i] = 0.0;
  stupidsort(locs, weights, 8);
  for (int i=0;i<8&&weights[i];i++)
    indices[i] = index(c, locs[i]);
  if (!contains(p) && weights[0]) {
    printf("Error at point %g %g\n", p.r(), p.z());
    printf("Interpolated to point %d %d\n", locs[0].r(), locs[0].z());
    printf("Or in other words... %g %g\n",
           operator[](locs[0]).r(), operator[](locs[0]).z());
    printf("I %s own the interpolated point.\n",
           owns(locs[0])?"actually":"don't");
    print();
    abort("Error made in interpolation of %s--fix this bug!!!\n",
          component_name(c));
  }
  // Throw out out of range indices:
  for (int i=0;i<8&&weights[i];i++)
    if (indices[0] < 0 || indices[0] >= ntot()) weights[i] = 0.0;
  // Stupid very crude code to compactify arrays:
  stupidsort(indices, weights, 8);
  if (!contains(p) && weights[0]) {
    printf("Error at point %g %g\n", p.r(), p.z());
    printf("Interpolated to point %d %d\n", locs[0].r(), locs[0].z());
    print();
    abort("Error made in interpolation of %s--fix this bug!!!\n",
          component_name(c));
  }
}

void volume::interpolate(component c, const vec &pc,
                         ivec locs[8], double weights[8]) const {
  const double SMALL = 1e-13;
  const vec p = (pc - yee_shift(c))*a;
  ivec middle(dim);
  LOOP_OVER_DIRECTIONS(dim,d)
    middle.set_direction(d, ((int) floor(p.in_direction(d)))*2+1);
  middle += iyee_shift(c);
  const vec midv = operator[](middle);
  const vec dv = (pc - midv)*(2*a);
  int already_have = 1;
  for (int i=0;i<8;i++) {
    locs[i] = round_vec(midv);
    weights[i] = 1.0;
  }
  LOOP_OVER_DIRECTIONS(dim,d) {
    for (int i=0;i<already_have;i++) {
      locs[already_have+i] = locs[i];
      weights[already_have+i] = weights[i];
      locs[i].set_direction(d,middle.in_direction(d)-1);
      weights[i] *= 0.5*(1.0-dv.in_direction(d));
      locs[already_have+i].set_direction(d,middle.in_direction(d)+1);
      weights[already_have+i] *= 0.5*(1.0+dv.in_direction(d));
    }
    already_have *= 2;
  }
  for (int i=already_have;i<8;i++) weights[i] = 0.0;
  double total_weight = 0.0;
  for (int i=0;i<already_have;i++) total_weight += weights[i];
  for (int i=0;i<already_have;i++)
    weights[i] += (1.0 - total_weight)*(1.0/already_have);
  for (int i=0;i<already_have;i++) {
    if (weights[i] < 0.0) {
      if (-weights[i] >= SMALL * 1e5)
        abort("large negative interpolation weight[%d] = %e\n", i, weights[i]);
      weights[i] = 0.0;
    }
    else if (weights[i] < SMALL)
      weights[i] = 0.0;
  }
  stupidsort(locs, weights, already_have);
  // The rest of this code is a crude hack to get the weights right when we
  // are exactly between a few grid points.  i.e. to eliminate roundoff
  // error.
  bool all_same = true;
  for (int i=0;i<8&&weights[i];i++)
    if (weights[i] != weights[0]) all_same = false;
  if (all_same) {
    int num_weights = 0;
    for (int i=0;i<8&&weights[i];i++) num_weights++;
    for (int i=0;i<8&&weights[i];i++) weights[i] = 1.0/num_weights;    
  }
}

geometric_volume empty_volume(ndim dim) {
  geometric_volume out(dim);
  LOOP_OVER_DIRECTIONS(dim,d) {
    out.set_direction_max(d,0.0);
    out.set_direction_min(d,0.0);
  }
  return out;
}

geometric_volume volume::dV(const ivec &here, double diameter) const {
  const double hinva = 0.5*inva * diameter;
  const volume &v = *this;
  const vec h = v[here];
  geometric_volume out(dim);
  LOOP_OVER_DIRECTIONS(dim,d) {
    out.set_direction_max(d,h.in_direction(d)+hinva);
    out.set_direction_min(d,h.in_direction(d)-hinva);
  }
  if (dim == Dcyl && here.r() == 0) {
    out.set_direction_min(R,0.0);
  }
  return out;
}

geometric_volume volume::dV(component c, int ind) const {
  if (!owns(iloc(c, ind))) return empty_volume(dim);
  return dV(iloc(c,ind));
}

double volume::xmax() const {
  const double qinva = 0.25*inva;
  return origin.x() + nx()*inva + qinva;
}

double volume::xmin() const {
  const double qinva = 0.25*inva;
  return origin.x() + qinva;
}

double volume::ymax() const {
  const double qinva = 0.25*inva;
  return origin.y() + ny()*inva + qinva;
}

double volume::ymin() const {
  const double qinva = 0.25*inva;
  return origin.y() + qinva;
}

double volume::zmax() const {
  const double qinva = 0.25*inva;
  return origin.z() + nz()*inva + qinva;
}

double volume::zmin() const {
  const double qinva = 0.25*inva;
  return origin.z() + qinva;
}

double volume::rmax() const {
  const double qinva = 0.25*inva;
  if (dim == Dcyl) return origin.r() + nr()*inva + qinva;
  abort("No rmax in these dimensions.\n");
  return 0.0; // This is never reached.
}

double volume::rmin() const {
  const double qinva = 0.25*inva;
  if (dim == Dcyl) {
    if (origin.r() == 0.0) {
      return 0.0;
    } else {
      return origin.r() + qinva;
    }
  }
  abort("No rmin in these dimensions.\n");
  return 0.0; // This is never reached.
}

double vec::project_to_boundary(direction d, double boundary_loc) {
  return fabs(boundary_loc - in_direction(d));
}

double volume::boundary_location(boundary_side b, direction d) const {
  // Returns the location of metallic walls...
  if (b == High) switch (d) {
  case X: return loc(Ez,ntot()-1).x();
  case Y: return loc(Ez,ntot()-1).y();
  case R: return loc(Ep,ntot()-1).r();
  case Z: if (dim == Dcyl) return loc(Ep,ntot()-1).z();
          else return loc(Ex,ntot()-1).z();
  case P: abort("P has no boundary!\n");
  case NO_DIRECTION: abort("NO_DIRECTION has no boundary!\n");
  }
  else switch (d) {
  case X: return loc(Ez,0).x();
  case Y: return loc(Ez,0).y();
  case R: return loc(Ep,0).r();
  case Z: if (dim == Dcyl) return loc(Ep,0).z();
          else return loc(Ex,0).z();
  case P: abort("P has no boundary!\n");
  case NO_DIRECTION: abort("NO_DIRECTION has no boundary!\n");
  }
  return 0.0;
}

ivec volume::big_corner() const {
  switch (dim) {
  case D1: return io + ivec(nz())*2;
  case D2: return io + ivec(nx(),ny())*2;
  case D3: return io + ivec(nx(),ny(),nz())*2;
  case Dcyl: return io + iveccyl(nr(),nz())*2;
  }
  return ivec(0); // This is never reached.
}

vec volume::corner(boundary_side b) const { 
  if (b == Low) return origin; // Low corner
  vec tmp = origin;
  LOOP_OVER_DIRECTIONS(dim, d)
    tmp.set_direction(d, tmp.in_direction(d) + num_direction(d) * inva);
  return tmp; // High corner
}

void volume::print() const {
  LOOP_OVER_DIRECTIONS(dim, d)
    printf("%s =%5g - %5g (%5g) \t", 
      direction_name(d), origin.in_direction(d), 
      origin.in_direction(d)+num_direction(d)/a, num_direction(d)/a); 
  printf("\n");
}

bool volume::intersect_with(const volume &vol_in, volume *intersection, volume *others, int *num_others) const {
  int temp_num[3] = {0,0,0};
  ivec new_io(dim);
  LOOP_OVER_DIRECTIONS(dim, d) {
    int minval = max(little_corner().in_direction(d), vol_in.little_corner().in_direction(d));
    int maxval = min(big_corner().in_direction(d), vol_in.big_corner().in_direction(d));
    if (minval >= maxval)
      return false;
    temp_num[d%3] = (maxval - minval)/2;
    new_io.set_direction(d, minval);
  }
  if (intersection != NULL) {
    *intersection = volume(dim, a, temp_num[0], temp_num[1], temp_num[2]); // fix me : ugly, need new constructor
    intersection->set_origin(new_io);
  }
  if (others != NULL) {
    int counter = 0;
    volume vol_containing = *this;
    LOOP_OVER_DIRECTIONS(dim, d) {
      if (vol_containing.little_corner().in_direction(d)
	  < vol_in.little_corner().in_direction(d)) {
	// shave off lower slice from vol_containing and add it to others
	volume other = vol_containing;
	const int thick = (vol_in.little_corner().in_direction(d)
			   - vol_containing.little_corner().in_direction(d))/2;
	other.set_num_direction(d, thick);
	others[counter] = other;
	counter++;
	vol_containing.shift_origin(d, thick*2);
	vol_containing.set_num_direction(d, vol_containing.num_direction(d)
					 - thick);
	if (vol_containing.little_corner().in_direction(d)
	    < vol_in.little_corner().in_direction(d))
	  abort("intersect_with: little corners differ by odd integer?");
      }
      if (vol_containing.big_corner().in_direction(d)
	  > vol_in.big_corner().in_direction(d)) {
	// shave off upper slice from vol_containing and add it to others
	volume other = vol_containing;
	const int thick = (vol_containing.big_corner().in_direction(d)
			   - vol_in.big_corner().in_direction(d))/2;
	other.set_num_direction(d, thick);
	other.shift_origin(d, (vol_containing.num_direction(d) - thick)*2);
	others[counter] = other;
	counter++;
	vol_containing.set_num_direction(d, vol_containing.num_direction(d) 
					 - thick);
	if (vol_containing.big_corner().in_direction(d)
	    < vol_in.big_corner().in_direction(d))
	  abort("intersect_with: big corners differ by odd integer?");
      }
    }
    *num_others = counter;
    
    int initial_points = 1;
    LOOP_OVER_DIRECTIONS(dim, d) initial_points *= num_direction(d);
    int final_points , temp = 1;
    LOOP_OVER_DIRECTIONS(dim, d) temp *= intersection->num_direction(d);    
    final_points = temp;
    for (int j=0; j<*num_others; j++) {
      temp = 1;
      LOOP_OVER_DIRECTIONS(dim, d) temp *= others[j].num_direction(d);
      final_points += temp;
    }
    if (initial_points != final_points)
      abort("intersect_with: initial_points != final_points,  %d, %d\n", 
	    initial_points, final_points);
  }
  return true;
}

vec volume::loc_at_resolution(int index, double res) const {
  vec where = origin;
  for (int dd=X;dd<=R;dd++) {
    const direction d = (direction) dd;
    if (has_boundary(High,d)) {
      const double dist = boundary_location(High,d)-boundary_location(Low,d);
      const int nhere = max(1,(int)floor(dist*res+0.5));
      where.set_direction(d,origin.in_direction(d) +
                          ((index % nhere)+0.5)*(1.0/res));
      index /= nhere;
    }
  }
  return where;
}

int volume::ntot_at_resolution(double res) const {
  int mytot = 1;
  for (int d=X;d<=R;d++)
    if (has_boundary(High,(direction)d)) {
      const double dist = boundary_location(High,(direction)d)
                        - boundary_location(Low,(direction)d);
      mytot *= max(1,(int)(dist*res+0.5));
    }
  return mytot;
}

vec volume::loc(component c, int ind) const {
  return operator[](iloc(c,ind));
}

ivec volume::iloc(component c, int ind) const {
  ivec out(dim);
  LOOP_OVER_DIRECTIONS(dim,d) {
    int ind_over_stride = ind/stride(d);
    while (ind_over_stride < 0) ind_over_stride += num_direction(d)+1;
    out.set_direction(d, 2*(ind_over_stride%(num_direction(d)+1)));
  }
  return out + iyee_shift(c) + io;
}

vec volume::dr() const {
  switch (dim) {
  case Dcyl: return veccyl(inva, 0.0);
  case D1: case D2: case D3: abort("Error in dr\n");
  }
  return vec(0); // This is never reached.
}

vec volume::dx() const {
  switch (dim) {
  case D3: return vec(inva,0,0);
  case D2: return vec(inva,0);
  case D1: case Dcyl: abort("Error in dx.\n");
  }
  return vec(0); // This is never reached.
}

vec volume::dy() const {
  switch (dim) {
  case D3: return vec(0,inva,0);
  case D2: return vec(0,inva);
  case D1: case Dcyl: abort("Error in dy.\n");
  }
  return vec(0); // This is never reached.
}

vec volume::dz() const {
  switch (dim) {
  case Dcyl: return veccyl(0.0,inva);
  case D3: return vec(0,0,inva);
  case D1: return vec(inva);
  case D2: abort("dz doesn't exist in 2D\n");
  }
  return vec(0); // This is never reached.
}

volume volone(double zsize, double a) {
  return volume(D1, a, 0, 0, (int) (zsize*a + 0.5));
}

volume voltwo(double xsize, double ysize, double a) {
  return volume(D2, a, (xsize==0)?1:(int) (xsize*a + 0.5),
                       (ysize==0)?1:(int) (ysize*a + 0.5),0);
}

volume vol1d(double zsize, double a) {
  return volone(zsize, a);
}

volume vol2d(double xsize, double ysize, double a) {
  return voltwo(xsize, ysize, a);
}

volume vol3d(double xsize, double ysize, double zsize, double a) {
  return volume(D3, a,(xsize==0)?1:(int) (xsize*a + 0.5),
                      (ysize==0)?1:(int) (ysize*a + 0.5),
                      (zsize==0)?1:(int) (zsize*a + 0.5));
}

volume volcyl(double rsize, double zsize, double a) {
  if (zsize == 0.0) return volume(Dcyl, a, (int) (rsize*a + 0.5), 0, 1);
  else return volume(Dcyl, a, (int) (rsize*a + 0.5), 0, (int) (zsize*a + 0.5));
}

volume volume::split(int n, int which) const {
  if (n > nowned_min())
    abort("Cannot split %d grid points into %d parts\n", nowned_min(), n);
  if (n == 1) return *this;

  // Try to get as close as we can...
  int biglen = 0;
  for (int i=0;i<3;i++) if (num[i] > biglen) biglen = num[i];
  const int split_point = (int)(biglen*(n/2)/(double)n + 0.5);
  const int num_low = (int)(split_point*n/(double)biglen + 0.5);
  if (which < num_low)
    return split_at_fraction(false, split_point).split(num_low,which);
  else
    return split_at_fraction(true, split_point).split(n-num_low,which-num_low);
}

volume volume::split_by_effort(int n, int which, int Ngv, const volume *gv, double *effort) const {
  const int grid_points_owned = nowned_min();
  if (n > grid_points_owned)
    abort("Cannot split %d grid points into %d parts\n", nowned_min(), n);
  if (n == 1) return *this;
  int biglen = 0;
  direction splitdir = NO_DIRECTION;
  LOOP_OVER_DIRECTIONS(dim, d) if (num_direction(d) > biglen) { biglen = num_direction(d); splitdir = d; } 
  double best_split_measure = 1e20, left_effort_fraction = 0;
  int best_split_point = 0;
  vec corner = zero_vec(dim);
  LOOP_OVER_DIRECTIONS(dim, d) corner.set_direction(d, origin.in_direction(d) + num_direction(d)/a); 

  for (int split_point = 1; split_point < biglen; split_point+=1) {
    volume v_left = *this;
    v_left.set_num_direction(splitdir, split_point);
    volume v_right = *this;
    v_right.set_num_direction(splitdir, num_direction(splitdir) - split_point);
    v_right.shift_origin(splitdir, split_point*2);

    double total_left_effort = 0, total_right_effort = 0;
    volume vol;
    if (Ngv == 0) {
      total_left_effort = v_left.ntot();
      total_right_effort = v_right.ntot();
    }
    else {
      for (int j = 0; j<Ngv; j++) {
	if (v_left.intersect_with(gv[j], &vol))
	  total_left_effort += effort[j] * vol.ntot();
	if (v_right.intersect_with(gv[j], &vol))
	  total_right_effort += effort[j] * vol.ntot();
      }
    }
    double split_measure = max(total_left_effort/(n/2), total_right_effort/(n-n/2));
    if (split_measure < best_split_measure) {
      best_split_measure = split_measure;
      best_split_point = split_point;
      left_effort_fraction = total_left_effort/(total_left_effort + total_right_effort);
    }
  }
  const int split_point = best_split_point;
    
  const int num_low = (int)(left_effort_fraction *n + 0.5);
  // Revert to split() when effort method gives less grid points than chunks
  if (num_low > best_split_point*(grid_points_owned/biglen) || 
      (n-num_low) > (grid_points_owned - best_split_point*(grid_points_owned/biglen)))
    return split(n, which);