fields.cpp

来自mit的fdtd开放性源代码meep

    else if (thef.f[c][cmp]) {
      f[c][cmp] = new double[v.ntot()];
      memcpy(f[c][cmp], thef.f[c][cmp], sizeof(double) * v.ntot());
    }
  }
  FOR_FIELD_TYPES(ft) {
    for (int ip=0;ip<3;ip++)
      num_connections[ft][ip][Incoming] 
	= num_connections[ft][ip][Outgoing] = 0;
    connection_phases[ft] = 0;
    for (int ip=0;ip<3;ip++) for (int io=0;io<2;io++)
      connections[ft][ip][io] = NULL;
    zeroes[ft] = NULL;
    num_zeroes[ft] = 0;
  }
  FOR_COMPONENTS(c) DOCMP2 
    if (thef.f_minus_p[c][cmp]) {
      f_minus_p[c][cmp] = new double[v.ntot()];
      memcpy(f_minus_p[c][cmp], thef.f_minus_p[c][cmp], 
	     sizeof(double) * v.ntot());
    }
  f_rderiv_int = NULL;
  figure_out_step_plan();
}

static inline bool cross_negative(direction a, direction b) {
  if (a >= R) a = direction(a - 3);
  if (b >= R) b = direction(b - 3);
  return ((3+b-a)%3) == 2;
}

static inline direction cross(direction a, direction b) {
  if (a == b) abort("bug - cross expects different directions");
  bool dcyl = a >= R || b >= R;
  if (a >= R) a = direction(a - 3);
  if (b >= R) b = direction(b - 3);
  direction c = direction((3+2*a-b)%3);
  if (dcyl && c < Z) return direction(c + 3);
  return c;
}

/* Call this whenever we modify the structure_chunk (fields_chunk::s) to
   implement copy-on-write semantics.  See also structure::changing_chunks. */
void fields_chunk::changing_structure() {
  if (s->refcount > 1) { // this chunk is shared, so make a copy
    s->refcount--;
    s = new structure_chunk(s);
  }
}

void fields_chunk::figure_out_step_plan() {
  FOR_COMPONENTS(cc)
    have_minus_deriv[cc] = have_plus_deriv[cc] = false;
  FOR_COMPONENTS(c1)
    if (f[c1][0]) {
      const direction dc1 = component_direction(c1);
      // Figure out which field components contribute.
      FOR_COMPONENTS(c2)
        if ((is_electric(c1) && is_magnetic(c2)) ||
            (is_D(c1) && is_magnetic(c2)) ||
            (is_magnetic(c1) && is_electric(c2)) ||
	    (is_B(c1) && is_electric(c2))) {
          const direction dc2 = component_direction(c2);
          if (dc1 != dc2 && v.has_field(c2) && v.has_field(c1) &&
              (has_direction(v.dim,cross(dc1,dc2)) ||
	       (v.dim == Dcyl && has_field_direction(v.dim,cross(dc1,dc2))))) {
            direction d_deriv = cross(dc1,dc2);
            if (cross_negative(dc2, dc1)) {
              minus_component[c1] = c2;
              have_minus_deriv[c1] = true;
              minus_deriv_direction[c1] = d_deriv;
            } else {
              plus_component[c1] = c2;
              have_plus_deriv[c1] = true;
              plus_deriv_direction[c1] = d_deriv;
            }
          }
        }
    }
  for (int i=0;i<3;i++) {
    num_each_direction[i] = v.yucky_num(i);
    stride_each_direction[i] = v.stride(v.yucky_direction(i));
  }
  FOR_DIRECTIONS(d) {
    num_any_direction[d] = 1;
    stride_any_direction[d] = 0;
    for (int i=0;i<3;i++)
      if (d == v.yucky_direction(i)) {
        num_any_direction[d] = v.yucky_num(i);
        stride_any_direction[d] = v.stride(v.yucky_direction(i));
      }
  }
}

static bool is_tm(component c) {
  switch (c) {
  case Hx: case Hy: case Bx: case By: case Ez: case Dz: return true;
  default: return false;
  }
  return false;
}

static bool is_like(ndim d, component c1, component c2) {
  if (d != D2) return true;
  return !(is_tm(c1) ^ is_tm(c2));
}

void fields_chunk::alloc_f(component the_c) {
  FOR_COMPONENTS(c)
    if (is_mine() && v.has_field(c) && is_like(v.dim, the_c, c)
	&& !is_magnetic(c))
      DOCMP {
        if (!f[c][cmp]) {
          f[c][cmp] = new double[v.ntot()];
          for (int i=0;i<v.ntot();i++) f[c][cmp][i] = 0.0;
	}
    }
  /* initially, we just set H == B ... later on, we lazily allocate
     H fields if needed (if mu != 1 or in PML) in update_h_from_b */
  FOR_H_AND_B(hc,bc) DOCMP 
    if (!f[hc][cmp] && f[bc][cmp]) f[hc][cmp] = f[bc][cmp];
  figure_out_step_plan();
}

void fields_chunk::remove_sources() {
  delete b_sources;
  delete d_sources;
  d_sources = b_sources = NULL;
}

void fields::remove_sources() {
  delete sources;
  sources = NULL;
  for (int i=0;i<num_chunks;i++) 
    chunks[i]->remove_sources();
}

void fields_chunk::remove_polarizabilities() {
  delete pol;
  delete olpol;
  pol = olpol = NULL;
  changing_structure();
  s->remove_polarizabilities();
}

void fields::remove_polarizabilities() {
  for (int i=0;i<num_chunks;i++) 
    chunks[i]->remove_polarizabilities();
}

void fields::remove_fluxes() {
  delete fluxes;
  fluxes = NULL;
}

void fields_chunk::zero_fields() {
  FOR_COMPONENTS(c) DOCMP {
    if (f[c][cmp]) for (int i=0;i<v.ntot();i++) f[c][cmp][i] = 0.0;
    if (f_backup[c][cmp]) for (int i=0;i<v.ntot();i++) f_backup[c][cmp][i] = 0.0;
    if (f_prev[c][cmp]) for (int i=0;i<v.ntot();i++) f_prev[c][cmp][i] = 0.0;
  }
  if (is_mine() && pol) pol->zero_fields();
  if (is_mine() && olpol) olpol->zero_fields();
}

void fields::zero_fields() {
  for (int i=0;i<num_chunks;i++)
    chunks[i]->zero_fields();
}

void fields::reset() {
  remove_sources();
  remove_fluxes();
  zero_fields();
  t = 0;
}

void fields_chunk::use_real_fields() {
  is_real = 1;
  // for mu=1 non-PML regions, H==B to save space/time - don't delete twice!
  FOR_H_AND_B(hc,bc) if (f[hc][1] == f[bc][1]) f[bc][1] = NULL;
  FOR_COMPONENTS(c) if (f[c][1]) {
    delete[] f[c][1];
    f[c][1] = 0;
  }
  if (is_mine() && pol) pol->use_real_fields();
  if (is_mine() && olpol) olpol->use_real_fields();
}

int fields::phase_in_material(const structure *snew, double time) {
  if (snew->num_chunks != num_chunks)
    abort("Can only phase in similar sets of chunks: %d vs %d\n", 
	  snew->num_chunks, num_chunks);
  for (int i=0;i<num_chunks;i++)
    if (chunks[i]->is_mine())
      chunks[i]->phase_in_material(snew->chunks[i]);
  phasein_time = (int) (time/dt);
  return phasein_time;
}

void fields_chunk::phase_in_material(const structure_chunk *snew) {
  new_s = snew;
}

int fields::is_phasing() {
  return phasein_time > 0;
}

// This is used for phasing the *radial origin* of a cylindrical structure
void fields::set_rshift(double rshift) {
  if (v.dim != Dcyl) abort("set_rshift is only for cylindrical coords");
  if (gv.in_direction_min(R) <= 0 && gv.in_direction_max(R) >= 0)
    abort("set_rshift is invalid if volume contains r=0");
  for (int i = 0; i < num_chunks; ++i)
    chunks[i]->rshift = rshift;
}

bool fields::equal_layout(const fields &f) const {
  if (a != f.a || 
      num_chunks != f.num_chunks ||
      gv != f.gv ||
      S != f.S)
    return false;
  for (int d=0;d<5;d++)
    if (k[d] != f.k[d])
      return false;
  for (int i = 0; i < num_chunks; ++i)
    if (chunks[i]->a != f.chunks[i]->a ||
	chunks[i]->gv != f.chunks[i]->gv)
      return false;
  return true;
}

// total computational volume, including regions redundant by symmetry
geometric_volume fields::total_volume(void) const {
  geometric_volume gv0 = v.interior();
  geometric_volume gv = gv0;
  for (int n = 1; n < S.multiplicity(); ++n)
    gv = gv | S.transform(gv0, n);
  if (gv.dim == Dcyl && gv.in_direction_min(R) < 0)
    gv.set_direction_min(R, 0);
  return gv;
}

/* One-pixel periodic dimensions are used almost exclusively to
   emulate lower-dimensional computations, so if the user passes an
   empty size in that direction, they probably really intended to
   specify that whole dimension.  This function detects that case. */
bool fields::nosize_direction(direction d) const {
  return (v.has_boundary(Low, d) && v.has_boundary(High, d) &&
	  boundaries[Low][d] == Periodic && boundaries[High][d] == Periodic
	  && v.num_direction(d) == 1);
}

} // namespace meep