structure.cpp

采用FDTD时域有限差分法计算电磁波的传播问题等。

  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 ((boundary_side) b == High)
    pml_volume.origin_set_direction(d, (user_volume.num_direction(d)
                                        - pml_volume.num_direction(d))/user_volume.a);
  const int v_to_user_shift = (int)
    floor((user_volume.origin_in_direction(d) - v.origin_in_direction(d))*v.a + 0.5);
  if ((boundary_side) 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[] kerr[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->kerr[c]) {
      kerr[c] = new double[v.ntot()];
      if (kerr[c] == NULL) abort("Out of memory!\n");
      for (int i=0;i<v.ntot();i++) kerr[c][i] = o->kerr[c][i];
    } else {
      kerr[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];
          }
      }
}

// The following is defined in anisotropic_averaging.cpp:
double anisoaverage(component ec, direction d, material_function &eps,
                    const geometric_volume &vol, double minvol);

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

  FOR_ELECTRIC_COMPONENTS(c)
    if (inveps[c][component_direction(c)]) {
      if (!kerr[c]) kerr[c] = new double[v.ntot()];
      LOOP_OVER_VOL(v, c, i) {
	IVEC_LOOP_LOC(v, here);
        kerr[c][i] = epsilon.kerr(here);
        LOOP_OVER_DIRECTIONS(v.dim,d)
          if (d != component_direction(c)) {
            vec dx = zero_vec(v.dim);
            dx.set_direction(d,0.5/a);
            // Following TD3D, we set kerr coefficient to zero if any of
            // the adjoining epsilon points is linear.
            kerr[c][i] = min(kerr[c][i], epsilon.kerr(here + dx));
            kerr[c][i] = min(kerr[c][i], epsilon.kerr(here - dx));
          }
      }
    }
}

void structure_chunk::set_epsilon(material_function &epsilon,
				  bool use_anisotropic_averaging,
				  double minvol) {
  if (!is_mine()) return;

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

  if (!eps) eps = new double[v.ntot()];
  LOOP_OVER_VOL(v, v.eps_component(), i) {
    IVEC_LOOP_LOC(v, here);
    eps[i] = epsilon.eps(here);
  }
  
  if (!use_anisotropic_averaging) {
    FOR_ELECTRIC_COMPONENTS(c)
      if (v.has_field(c)) {
	bool have_other_direction = false;
	vec dxa = zero_vec(v.dim);
	vec dxb = zero_vec(v.dim);
	direction c_d = component_direction(c);
	LOOP_OVER_DIRECTIONS(v.dim,da)
	  if (da != c_d) {
	    dxa.set_direction(da,0.5/a);
	    LOOP_OVER_DIRECTIONS(v.dim,db)
	      if (db != c_d && db != da) {
		dxb.set_direction(db,0.5/a);
		have_other_direction = true;
	      }
	    break;
	  }
	if (!inveps[c][c_d]) inveps[c][c_d] = new double[v.ntot()];
	LOOP_OVER_VOL(v, c, i) {
	  IVEC_LOOP_LOC(v, here);
	  if (!have_other_direction)
	    inveps[c][c_d][i] =
	      2.0/(epsilon.eps(here + dxa) + epsilon.eps(here - dxa));
	  else
	    inveps[c][c_d][i] = 4.0/(epsilon.eps(here + dxa + dxb) +
				     epsilon.eps(here + dxa - dxb) +
				     epsilon.eps(here - dxa + dxb) +
				     epsilon.eps(here - dxa - dxb));
	}
      }
  } else {
    if (minvol == 0.0) minvol = v.dV(zero_ivec(v.dim)).full_volume()/100.0;
    FOR_ELECTRIC_COMPONENTS(c)
      if (v.has_field(c))
	FOR_ELECTRIC_COMPONENTS(c2) if (v.has_field(c2)) {
	  direction d = component_direction(c2);
          if (!inveps[c][d]) inveps[c][d] = new double[v.ntot()];
          if (!inveps[c][d]) abort("Memory allocation error.\n");
	  LOOP_OVER_VOL(v, c, i) {
	    IVEC_LOOP_ILOC(v, here);
	    inveps[c][d][i] = anisoaverage(c, d, epsilon, v.dV(here), minvol);
	  }
      }
  }

  update_pml_arrays(); // PML stuff depends on epsilon
}

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) kerr[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();
}

} // namespace meep