structure.cpp

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

      chunks[i]->set_chi3(eps);
}

void structure::set_chi3(double eps(const vec &)) {
  simple_material_function epsilon(eps);
  set_chi3(epsilon);
}

void structure::set_chi2(material_function &eps) {
  changing_chunks();
  for (int i=0;i<num_chunks;i++)
    if (chunks[i]->is_mine())
      chunks[i]->set_chi2(eps);
}

void structure::set_chi2(double eps(const vec &)) {
  simple_material_function epsilon(eps);
  set_chi2(epsilon);
}

void structure::use_pml(direction d, boundary_side b, double dx,
			double strength) {
  if (strength == 0.0 || dx <= 0.0) return;
  volume pml_volume = v;
  pml_volume.set_num_direction(d, int(dx*user_volume.a + 1 + 0.5)); //FIXME: exact value?
  if (b == High)
    pml_volume.set_origin(d, user_volume.big_corner().in_direction(d)
  			  - pml_volume.num_direction(d) * 2);
  const int v_to_user_shift = (user_volume.little_corner().in_direction(d) 
			       - v.little_corner().in_direction(d)) / 2;
  if (b == Low && v_to_user_shift != 0)
    pml_volume.set_num_direction(d, pml_volume.num_direction(d) + v_to_user_shift);
  add_to_effort_volumes(pml_volume, 0.60); // FIXME: manual value for pml effort
}

void structure::mix_with(const structure *oth, double f) {
  if (num_chunks != oth->num_chunks)
    abort("You can't phase materials with different chunk topologies...\n");
  changing_chunks();
  for (int i=0;i<num_chunks;i++)
    if (chunks[i]->is_mine())
      chunks[i]->mix_with(oth->chunks[i], f);
}

structure_chunk::~structure_chunk() {
  FOR_ELECTRIC_COMPONENTS(c) {
    FOR_DIRECTIONS(d)
      delete[] inveps[c][d];
    delete[] chi3[c];
    delete[] chi2[c];
  }
  delete[] eps;

  FOR_COMPONENTS(c) FOR_DIRECTIONS(d) delete[] C[d][c];
  FOR_COMPONENTS(c)
    FOR_DIRECTIONS(d) FOR_DIRECTIONS(d2)
        delete[] Cdecay[d][c][d2];
  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_ELECTRIC_COMPONENTS(c) FOR_DIRECTIONS(d)
    if (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]);
  // 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;
  }
}

const double Cmax = 0.5;

void structure_chunk::use_pml(direction d, double dx, double bloc,
			      double strength) {
  if (strength == 0.0 || dx <= 0.0) return;
  const double prefac = strength * Cmax/(dx*dx);
  // 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 ||
      bloc < v.boundary_location(Low,d) - dx - 1.0/a) return;
  if (is_mine()) {
    FOR_COMPONENTS(c)
      if (v.has_field(c) && component_direction(c) != d) {
        if (!C[d][c]) {
          C[d][c] = new double[v.ntot()];
          for (int i=0;i<v.ntot();i++) C[d][c][i] = 0.0;
        }
	LOOP_OVER_VOL(v, (component) c, i) {
	  IVEC_LOOP_LOC(v, here);
          const double x =
            0.5/a*((int)(dx*(2*a)+0.5) -
                   (int)(2*a*fabs(bloc-here.in_direction(d))+0.5));
          if (x > 0) C[d][c][i] = prefac*x*x;
        }
      }
  }
}

void structure_chunk::update_pml_arrays() {
  if (!is_mine()) return;
  FOR_DIRECTIONS(d)
    FOR_COMPONENTS(c) 
      if (C[d][c] != NULL) {
	bool all_zeros = true;
	LOOP_OVER_VOL_OWNED(v, c, i) {
	  if (C[d][c][i] != 0.0) {
	    all_zeros = false;
	    goto done; // can't use 'break' since LOOP has nested loops
	  }
	}
      done:
	if (all_zeros) {
	  delete[] C[d][c];
	  C[d][c] = NULL;
	}
      }
  update_Cdecay();
}

void structure_chunk::update_Cdecay() {
  FOR_DIRECTIONS(d) FOR_COMPONENTS(c) FOR_DIRECTIONS(d2) {
    delete[] Cdecay[d][c][d2];
    Cdecay[d][c][d2] = NULL;
  }
  FOR_DIRECTIONS(d) FOR_COMPONENTS(c) 
    if (C[d][c] != NULL)
      FOR_DIRECTIONS(d2) 
	if ((inveps[c][d2] || d2 == component_direction(c)) && d2 != d) {
	  Cdecay[d][c][d2] = new double[v.ntot()];
	  for (int i=0;i<v.ntot();i++) {
	    if (is_magnetic(c)) Cdecay[d][c][d2][i] = 1.0/(1.0+0.5*C[d][c][i]);
	    else if (is_electric(c))
          Cdecay[d][c][d2][i] = inveps[c][d2][i]/(1.0+0.5*C[d][c][i]*inveps[c][d2][i]);
        else Cdecay[d][c][d2][i] = 1.0/(1.0+0.5*C[d][c][i]); // FIXME: Is this right?
	    /*
	      TODO: THIS IS NOT CORRECT FOR NON-DIAGONAL INVEPS...maybe
              above code is not correct either????
	    if (Cdecay[d][c][d2][i] == 0.0)
	      abort("In update_Cdecay: Cdecay == 0\n");
	    */
	  }
	}
}

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() && o->inveps[c][d]) {
      inveps[c][d] = new double[v.ntot()];
      for (int i=0;i<v.ntot();i++) inveps[c][d][i] = o->inveps[c][d][i];
    } else {
      inveps[c][d] = NULL;
    }
  // Allocate the conductivity arrays:
  FOR_DIRECTIONS(d) FOR_COMPONENTS(c) C[d][c] = NULL;
  FOR_DIRECTIONS(d) FOR_COMPONENTS(c) FOR_DIRECTIONS(d2)
    Cdecay[d][c][d2] = NULL;
  // Copy over the conductivity arrays:
  if (is_mine())
    FOR_DIRECTIONS(d) FOR_COMPONENTS(c)
      if (o->C[d][c]) {
        C[d][c] = new double[v.ntot()];
        for (int i=0;i<v.ntot();i++) C[d][c][i] = o->C[d][c][i];
        FOR_DIRECTIONS(d2)
          if (o->Cdecay[d][c][d2]) {
            Cdecay[d][c][d2] = new double[v.ntot()];
            for (int i=0;i<v.ntot();i++)
              Cdecay[d][c][d2][i] = o->Cdecay[d][c][d2][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()];
      LOOP_OVER_VOL(v, c, i) {
	IVEC_LOOP_LOC(v, here);
        chi3[c][i] = epsilon.chi3(here);
      }
      
      /* currently, our update_e_from_d routine requires that
	 chi2 be present if chi3 is, and vice versa */
      if (!chi2[c]) {
	chi2[c] = new double[v.ntot()]; 
	memset(chi2[c], 0, v.ntot() * sizeof(double)); // chi2 = 0 by default
      }
    }
}

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()];
      LOOP_OVER_VOL(v, c, i) {
	IVEC_LOOP_LOC(v, here);
        chi2[c][i] = epsilon.chi2(here);
      }

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

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;
  FOR_DIRECTIONS(d) FOR_COMPONENTS(c) C[d][c] = NULL;
  FOR_DIRECTIONS(d) FOR_DIRECTIONS(d2) FOR_COMPONENTS(c) Cdecay[d][c][d2] = NULL;
}

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