structure.cpp

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

structure_chunk::~structure_chunk() {
  FOR_COMPONENTS(c) {
    FOR_DIRECTIONS(d) {
      delete[] inveps[c][d];
      delete[] invmu[c][d];
      delete[] conductivity[c][d];
      delete[] condinv[c][d];
    }
    delete[] chi2[c];
    delete[] chi3[c];
  }
  delete[] eps;
  FOR_DIRECTIONS(d) { 
    delete[] sig[d];
    delete[] siginv[d];
  }
  if (pb) delete pb;
}

void structure_chunk::mix_with(const structure_chunk *n, double f) {
  for (int i=0;i<v.ntot();i++)
    eps[i] = 1.0/(1.0/eps[i] + f*(1.0/n->eps[i]-1.0/eps[i]));
  FOR_COMPONENTS(c) FOR_DIRECTIONS(d) {
    if (!inveps[c][d] && n->inveps[c][d]) {
      inveps[c][d] = new double[v.ntot()];
      if (component_direction(c) == d) // diagonal components = 1 by default
	for (int i=0;i<v.ntot();i++) inveps[c][d][i] = 1.0;
      else
	for (int i=0;i<v.ntot();i++) inveps[c][d][i] = 0.0;
    }
    if (!invmu[c][d] && n->invmu[c][d]) {
      invmu[c][d] = new double[v.ntot()];
      if (component_direction(c) == d) // diagonal components = 1 by default
	for (int i=0;i<v.ntot();i++) invmu[c][d][i] = 1.0;
      else
	for (int i=0;i<v.ntot();i++) invmu[c][d][i] = 0.0;
    }
    if (!conductivity[c][d] && n->conductivity[c][d]) {
      conductivity[c][d] = new double[v.ntot()];
      for (int i=0;i<v.ntot();i++) conductivity[c][d][i] = 0.0;
    }
    if (inveps[c][d]) {
      if (n->inveps[c][d])
	for (int i=0;i<v.ntot();i++)
	  inveps[c][d][i] += f*(n->inveps[c][d][i] - inveps[c][d][i]);
      else {
	double nval = component_direction(c) == d ? 1.0 : 0.0; // default
	for (int i=0;i<v.ntot();i++)
	  inveps[c][d][i] += f*(nval - inveps[c][d][i]);
      }
    }
    if (invmu[c][d]) {
      if (n->invmu[c][d])
	for (int i=0;i<v.ntot();i++)
	  invmu[c][d][i] += f*(n->invmu[c][d][i] - invmu[c][d][i]);
      else {
	double nval = component_direction(c) == d ? 1.0 : 0.0; // default
	for (int i=0;i<v.ntot();i++)
	  invmu[c][d][i] += f*(nval - invmu[c][d][i]);
      }
    }
    if (conductivity[c][d]) {
      if (n->conductivity[c][d])
	for (int i=0;i<v.ntot();i++)
	  conductivity[c][d][i] += f*(n->conductivity[c][d][i]
				      - conductivity[c][d][i]);
      else
	for (int i=0;i<v.ntot();i++)
	  conductivity[c][d][i] += f*(0.0 - conductivity[c][d][i]);
    }
    condinv_stale = true;
  }
  // Mix in the polarizability...
  polarizability *po = pb, *pn = n->pb;
  while (po && pn) {
    FOR_COMPONENTS(c)
      if (v.has_field(c) && is_electric(c))
        for (int i=0;i<v.ntot();i++)
          po->s[c][i] += f*(pn->s[c][i] - po->s[c][i]);
    po = po->next;
    pn = pn->next;
  }
}

void structure_chunk::use_pml(direction d, double dx, double bloc,
			      double Rasymptotic,
			      pml_profile_func pml_profile,
			      void *pml_profile_data,
			      double pml_profile_integral) {
  if (dx <= 0.0) return;
  const double prefac = (-log(Rasymptotic))/(4*dx*pml_profile_integral);
  // Don't bother with PML if we don't even overlap with the PML region...
  if (bloc > v.boundary_location(High,d) + dx + 1.0/a - 1e-10 ||
      bloc < v.boundary_location(Low,d) - dx - 1.0/a + 1e-10) return;
  if (is_mine()) {
    if (sig[d]) {
      delete[] sig[d]; 
      delete[] siginv[d];
      sig[d] = NULL;
      siginv[d] = NULL;
    }
    LOOP_OVER_FIELD_DIRECTIONS(v.dim, dd) {
      if (!sig[dd]) {
	int spml = (dd==d)?(2*v.num_direction(d)+2):1;
	sigsize[dd] = spml;
	sig[dd] = new double[spml];
	siginv[dd] = new double[spml];
	for (int i=0;i<spml;++i) {
	  sig[dd][i] = 0.0;
	  siginv[dd][i] = 1.0;
	}
      }
    }
    for (int i=v.little_corner().in_direction(d);
	 i<=v.big_corner().in_direction(d)+1;++i) {
      int idx = i - v.little_corner().in_direction(d);
      double here = i * 0.5/a;
      const double x =
	0.5/a*((int)(dx*(2*a)+0.5) - (int)(fabs(bloc-here)*(2*a)+0.5));
      if (x > 0) {
	sig[d][idx]=0.5*dt*prefac*pml_profile(x/dx, pml_profile_data);
	siginv[d][idx] = 1/(1+sig[d][idx]);	
      }
    }
  }
  condinv_stale = true;
}

void structure_chunk::update_condinv() {
  if (!condinv_stale || !is_mine()) return;
  FOR_COMPONENTS(c) {
    direction d = component_direction(c);
    if (conductivity[c][d]) {
      if (!condinv[c][d]) condinv[c][d] = new double[v.ntot()];
      const direction dsig = cycle_direction(v.dim,d,1);
      const bool have_pml = sigsize[dsig] > 1;
      if (!have_pml) {
	LOOP_OVER_VOL(v, c, i)
	  condinv[c][d][i] = 1 / (1 + conductivity[c][d][i] * dt * 0.5);
      }
      else { // include PML conductivity in condinv
	int k0 = v.little_corner().in_direction(dsig);
	LOOP_OVER_VOL(v, c, i) {
	  IVEC_LOOP_ILOC(v, iloc);
          int k = iloc.in_direction(dsig) - k0;
	  condinv[c][d][i] = 1 / (1 + conductivity[c][d][i] * dt * 0.5
				  + sig[dsig][k]);
	}
      }
    }
    else if (condinv[c][d]) { // condinv not needed
      delete[] condinv[c][d];
      condinv[c][d] = NULL;
    }
  }
  condinv_stale = false;
}

structure_chunk::structure_chunk(const structure_chunk *o) : gv(o->gv) {
  refcount = 1;
  if (o->pb) pb = new polarizability(o->pb);
  else pb = NULL;
  a = o->a;
  Courant = o->Courant;
  dt = o->dt;
  v = o->v;
  the_proc = o->the_proc;
  the_is_mine = my_rank() == n_proc();
  if (is_mine() && o->eps) {
    eps = new double[v.ntot()];
    if (eps == NULL) abort("Out of memory!\n");
    for (int i=0;i<v.ntot();i++) eps[i] = o->eps[i];
  } else {
    eps = NULL;
  }
  FOR_COMPONENTS(c) {
    if (is_mine() && o->chi3[c]) {
      chi3[c] = new double[v.ntot()];
      if (chi3[c] == NULL) abort("Out of memory!\n");
      for (int i=0;i<v.ntot();i++) chi3[c][i] = o->chi3[c][i];
    } else {
      chi3[c] = NULL;
    }
    if (is_mine() && o->chi2[c]) {
      chi2[c] = new double[v.ntot()];
      if (chi2[c] == NULL) abort("Out of memory!\n");
      for (int i=0;i<v.ntot();i++) chi2[c][i] = o->chi2[c][i];
    } else {
      chi2[c] = NULL;
    }
  }
  FOR_COMPONENTS(c) FOR_DIRECTIONS(d) if (is_mine()) {
    if (o->inveps[c][d]) {
      inveps[c][d] = new double[v.ntot()];
      memcpy(inveps[c][d], o->inveps[c][d], v.ntot()*sizeof(double));
    } else inveps[c][d] = NULL;
    if (o->invmu[c][d]) {
      invmu[c][d] = new double[v.ntot()];
      memcpy(invmu[c][d], o->invmu[c][d], v.ntot()*sizeof(double));
    } else invmu[c][d] = NULL;
    if (o->conductivity[c][d]) {
      conductivity[c][d] = new double[v.ntot()];
      memcpy(conductivity[c][d], o->conductivity[c][d],
	     v.ntot()*sizeof(double));
      condinv[c][d] = new double[v.ntot()];
      memcpy(condinv[c][d], o->condinv[c][d], v.ntot()*sizeof(double));
    } else conductivity[c][d] = condinv[c][d] = NULL;
  }
  condinv_stale = o->condinv_stale;
  // Allocate the PML conductivity arrays:
  FOR_DIRECTIONS(d) { 
    sig[d] = NULL; 
    siginv[d] = NULL;
    sigsize[d] = 0;
  }
  for (int i=0;i<5;++i) sigsize[i] = 0;
  // Copy over the PML conductivity arrays:
  if (is_mine())
    FOR_DIRECTIONS(d) 
      if (o->sig[d]) {
	sig[d] = new double[2*v.num_direction(d)+1];
	siginv[d] = new double[2*v.num_direction(d)+1];
	sigsize[d] = o->sigsize[d];
	for (int i=0;i<2*v.num_direction(d)+1;i++) {
	  sig[d][i] = o->sig[d][i];
	  siginv[d][i] = o->sig[d][i];
	}
      }
}

void structure_chunk::set_chi3(material_function &epsilon) {
  if (!is_mine()) return;
  
  epsilon.set_volume(v.pad().surroundings());

  FOR_ELECTRIC_COMPONENTS(c)
    if (inveps[c][component_direction(c)]) {
      if (!chi3[c]) chi3[c] = new double[v.ntot()];
      bool trivial = true;
      LOOP_OVER_VOL(v, c, i) {
	IVEC_LOOP_LOC(v, here);
        chi3[c][i] = epsilon.chi3(here);
	trivial = trivial && (chi3[c][i] == 0.0);
      }
      
      /* currently, our update_e_from_d routine requires that
	 chi2 be present if chi3 is, and vice versa */
      if (!chi2[c]) {
	if (!trivial) {
	  chi2[c] = new double[v.ntot()]; 
	  memset(chi2[c], 0, v.ntot() * sizeof(double)); // chi2 = 0
	}
	else { // no chi3, and chi2 is trivial (== 0), so delete
	  delete[] chi3[c];
	  chi3[c] = NULL;
	}
      }
    }
}

void structure_chunk::set_chi2(material_function &epsilon) {
  if (!is_mine()) return;
  
  epsilon.set_volume(v.pad().surroundings());

  FOR_ELECTRIC_COMPONENTS(c)
    if (inveps[c][component_direction(c)]) {
      if (!chi2[c]) chi2[c] = new double[v.ntot()];
      bool trivial = true;
      LOOP_OVER_VOL(v, c, i) {
	IVEC_LOOP_LOC(v, here);
        chi2[c][i] = epsilon.chi2(here);
	trivial = trivial && (chi2[c][i] == 0.0);
      }

      /* currently, our update_e_from_d routine requires that
	 chi3 be present if chi2 is, and vice versa */
      if (!chi3[c]) {
	if (!trivial) {
	  chi3[c] = new double[v.ntot()]; 
	  memset(chi3[c], 0, v.ntot() * sizeof(double)); // chi3 = 0
	}
	else { // no chi2, and chi3 is trivial (== 0), so delete 
	  delete[] chi2[c];
          chi2[c] = NULL;
	}
      }
    }
}

void structure_chunk::set_conductivity(component c, material_function &C) {
  if (!is_mine() || !v.has_field(c)) return;

  C.set_volume(v.pad().surroundings());

  if (!is_electric(c) && !is_magnetic(c) && !is_D(c) && !is_B(c))
    abort("invalid component for conductivity");

  direction c_d = component_direction(c);
  component c_C = is_electric(c) ? direction_component(Dx, c_d) :
    (is_magnetic(c) ? direction_component(Bx, c_d) : c);
  double *multby = is_electric(c) ? inveps[c][c_d] :
    (is_magnetic(c) ? invmu[c][c_d] : NULL);
  if (!conductivity[c_C][c_d]) 
    conductivity[c_C][c_d] = new double[v.ntot()]; 
  if (!conductivity[c_C][c_d]) abort("Memory allocation error.\n");
  bool trivial = true;
  double *cnd = conductivity[c_C][c_d];
  if (multby) {
    LOOP_OVER_VOL(v, c_C, i) {
      IVEC_LOOP_LOC(v, here);
      cnd[i] = C.conductivity(c, here) * multby[i];
      trivial = trivial && (cnd[i] == 0.0);
    }
  }
  else {
    LOOP_OVER_VOL(v, c_C, i) {
      IVEC_LOOP_LOC(v, here);
      cnd[i] = C.conductivity(c, here);
      trivial = trivial && (cnd[i] == 0.0);
    }
  }
  if (trivial) { // skip conductivity computations if conductivity == 0
    delete[] conductivity[c_C][c_d];
    conductivity[c_C][c_d] = NULL;
  }
  condinv_stale = true;
}

structure_chunk::structure_chunk(const volume &thev, 
				 const geometric_volume &vol_limit, 
				 double Courant, int pr)
  : Courant(Courant), gv(thev.surroundings() & vol_limit) {
  refcount = 1;
  pml_fmin = 0.2;
  pb = NULL;
  v = thev;
  a = thev.a;
  dt = Courant/a;
  the_proc = pr;
  the_is_mine = n_proc() == my_rank();

  // initialize materials arrays to NULL
  eps = NULL;
  FOR_COMPONENTS(c) chi3[c] = NULL;
  FOR_COMPONENTS(c) chi2[c] = NULL;
  FOR_COMPONENTS(c) FOR_DIRECTIONS(d) {
    inveps[c][d] = NULL;
    invmu[c][d] = NULL;
    conductivity[c][d] = NULL;
    condinv[c][d] = NULL;
  }
  condinv_stale = false;
  FOR_DIRECTIONS(d) { 
    sig[d] = NULL; 
    siginv[d] = NULL;
    sigsize[d] = 0;
  }
}

double structure::max_eps() const {
  double themax = 0.0;
  for (int i=0;i<num_chunks;i++)
    if (chunks[i]->is_mine())
      themax = max(themax,chunks[i]->max_eps());
  return max_to_all(themax);
}

double fields::max_eps() const {
  double themax = 0.0;
  for (int i=0;i<num_chunks;i++)
    if (chunks[i]->is_mine())
      themax = max(themax,chunks[i]->s->max_eps());
  return max_to_all(themax);
}

double structure_chunk::max_eps() const {
  double themax = 0.0;
  for (int i=0;i<v.ntot();i++) themax = max(themax,eps[i]);
  return themax;
}

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

void structure_chunk::remove_polarizabilities() {
  delete pb;
  pb = NULL;
}

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

// for debugging, display the chunk layout
void structure::print_layout(void) const {
  direction d0 = v.yucky_direction(0);
  direction d1 = v.yucky_direction(1);
  direction d2 = v.yucky_direction(2);
  for (int i = 0; i < num_chunks; ++i) {
    master_printf("chunk[%d] on process %d, resolution %g (%s,%s,%s):"
		  " (%d,%d,%d) - (%d,%d,%d)\n",
		  i, chunks[i]->n_proc(), chunks[i]->a,
		  direction_name(d0),direction_name(d1),direction_name(d2),
		  chunks[i]->v.little_corner().yucky_val(0),
		  chunks[i]->v.little_corner().yucky_val(1),
		  chunks[i]->v.little_corner().yucky_val(2),
		  chunks[i]->v.big_corner().yucky_val(0),
		  chunks[i]->v.big_corner().yucky_val(1),
		  chunks[i]->v.big_corner().yucky_val(2));
  }
}

} // namespace meep