boundaries.cpp

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

/* Copyright (C) 2006 Massachusetts Institute of Technology
%
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2, or (at your option)
%  any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software Foundation,
%  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
*/

#include <stdlib.h>
#include <complex>

#include "meep.hpp"
#include "meep_internals.hpp"

#define UNUSED(x) (void) x // silence compiler warnings

namespace meep {

void fields::set_boundary(boundary_side b,direction d,boundary_condition cond){
  if (boundaries[b][d] != cond) {
    boundaries[b][d] = cond;
    chunk_connections_valid = false;
  }
}

void fields::use_bloch(direction d, complex<double> kk) {
  k[d] = kk;
  for (int b=0;b<2;b++) set_boundary(boundary_side(b), d, Periodic);
  if (real(kk) * v.num_direction(d) == 0.5 * a) // check b.z. edge exactly
    eikna[d] = -exp(-imag(kk) * ((2*pi/a)*v.num_direction(d)));
  else {
    const complex<double> I = complex<double>(0.0,1.0);
    eikna[d] = exp(I*kk*((2*pi/a)*v.num_direction(d)));
  }
  coskna[d] = real(eikna[d]);
  sinkna[d] = imag(eikna[d]);
  if (is_real && kk != 0.0) // FIXME: allow real phases (c.f. CONNECT_PHASE)
    abort("Can't use real fields with bloch boundary conditions!\n");
  chunk_connections_valid = false; // FIXME: we don't always need to invalidate
}

void fields::use_bloch(const vec &k) {
  // Note that I allow a 1D k input when in cylindrical, since in that case
  // it is unambiguous.
  if (k.dim != v.dim && !(k.dim == D1 && v.dim == Dcyl))
    abort("Aaaack, k has wrong dimensions!\n");
  LOOP_OVER_DIRECTIONS(v.dim, d)
    if (v.has_boundary(Low,d) && d != R)
      use_bloch(d, k.in_direction(d));
}

ivec fields::ilattice_vector(direction d) const {
  switch (user_volume.dim) {
  case D1: return ivec(2*user_volume.nz());
  case Dcyl: return iveccyl(0,2*user_volume.nz()); // Only Z direction here
  case D2:
    switch (d) {
    case X: return ivec(user_volume.nx()*2,0);
    case Y: return ivec(0,user_volume.ny()*2);
    case Z: case R: case P: case NO_DIRECTION: break;
    }
  case D3:
    switch (d) {
    case X: return ivec(user_volume.nx()*2,0,0);
    case Y: return ivec(0,user_volume.ny()*2,0);
    case Z: return ivec(0,0,user_volume.nz()*2);
    case R: case P: case NO_DIRECTION: break;
    }
  }
  abort("Aaack in ilattice_vector.\n");
  return ivec(0);
}

vec fields::lattice_vector(direction d) const {
  return v[ilattice_vector(d)];
}

void fields::disconnect_chunks() {
  chunk_connections_valid = false;
  for (int i=0;i<num_chunks;i++) {
    DOCMP {
      FOR_FIELD_TYPES(f)
	for (int ip=0;ip<3;++ip)
	  for (int io=0;io<2;io++) {
	    delete[] chunks[i]->connections[f][ip][io];
	    chunks[i]->connections[f][ip][io] = NULL;
	  }
    }
    FOR_FIELD_TYPES(f) {
      delete[] chunks[i]->connection_phases[f];
      chunks[i]->connection_phases[f] = NULL;
      for (int ip=0;ip<3;++ip)
	for (int io=0;io<2;io++)
	  chunks[i]->num_connections[f][ip][io] = 0;
    }
  }
  FOR_FIELD_TYPES(ft)
    for (int i=0;i<num_chunks*num_chunks;i++) {
      delete[] comm_blocks[ft][i];
      comm_blocks[ft][i] = 0;
      for (int ip=0;ip<3;++ip)
	comm_sizes[ft][ip][i] = 0;
    }
}

void fields::connect_chunks() {
  if (!chunk_connections_valid) {
    am_now_working_on(Connecting);
    disconnect_chunks();
    find_metals();
    connect_the_chunks();
    finished_working();
    chunk_connections_valid = true;
  }
}

inline bool fields::on_metal_boundary(const ivec &here) {
  LOOP_OVER_DIRECTIONS(v.dim, d) {
    if (user_volume.has_boundary(High, d) &&
        here.in_direction(d) == user_volume.big_corner().in_direction(d)) {
      if (boundaries[High][d] == Metallic) return true;
    }
    if (boundaries[Low][d] == Magnetic &&
        here.in_direction(d) ==
        user_volume.little_corner().in_direction(d)+1)
      return true;
    if (boundaries[Low][d] == Metallic &&
        here.in_direction(d) ==
        user_volume.little_corner().in_direction(d))
      return true;
  }
  return false;
}

bool fields::locate_point_in_user_volume(ivec *there, complex<double> *phase) const {
  // Check if a translational symmetry is needed to bring the point in...
  if (!user_volume.owns(*there))
    FOR_DIRECTIONS(d) {
      if (boundaries[High][d] == Periodic &&
          there->in_direction(d) <= user_volume.little_corner().in_direction(d)) {
        while (there->in_direction(d) <=
               user_volume.little_corner().in_direction(d)) {
          *there += ilattice_vector(d);
          *phase *= conj(eikna[d]);
        }
      } else if (boundaries[High][d] == Periodic &&
                 there->in_direction(d)-ilattice_vector(d).in_direction(d)
                 > user_volume.little_corner().in_direction(d)) {
        while (there->in_direction(d)-ilattice_vector(d).in_direction(d)
               > user_volume.little_corner().in_direction(d)) {
          *there -= ilattice_vector(d);
          *phase *= eikna[d];
        }
      }
    }
  return user_volume.owns(*there);
}

void fields::locate_volume_source_in_user_volume(const vec p1, const vec p2, vec newp1[8], vec newp2[8],
                                                  complex<double> kphase[8], int &ncopies) const {
  // For periodic boundary conditions, 
  // this function locates up to 8 translated copies of the initial volume specified by (p1,p2)
  // First bring center of volume inside
  ncopies = 1;
  newp1[0] = p1;
  newp2[0] = p2;
  kphase[0] = 1;
  vec cen = (newp1[0] + newp2[0]) * 0.5;
  LOOP_OVER_DIRECTIONS(v.dim, d) 
    if (boundaries[High][d] == Periodic)  {
      while (cen.in_direction(d) < v.boundary_location(Low, d)) {
        newp1[0] += lattice_vector(d);
        newp2[0] += lattice_vector(d);
        kphase[0] *= conj(eikna[d]);
        cen = (newp1[0] + newp2[0]) * 0.5;
      }
      while (cen.in_direction(d) > v.boundary_location(High, d)) {
        newp1[0] -= lattice_vector(d);
        newp2[0] -= lattice_vector(d);
        kphase[0] *= eikna[d];
        cen = (newp1[0] + newp2[0]) * 0.5;
      }
    }
  
  // if volume extends outside user_volume in any direction, we need to duplicate already existing copies
  LOOP_OVER_DIRECTIONS(v.dim, d) 
    if (boundaries[High][d] == Periodic) {
      if (newp1[0].in_direction(d) < v.boundary_location(Low, d) ||
          newp2[0].in_direction(d) < v.boundary_location(Low, d)) {
        for (int j=0; j<ncopies; j++) {
          newp1[ncopies+j] = newp1[j] + lattice_vector(d);
          newp2[ncopies+j] = newp2[j] + lattice_vector(d);
          kphase[ncopies+j] = kphase[j] * conj(eikna[d]);
        }
        ncopies *= 2;
      }
      else if (newp1[0].in_direction(d) > v.boundary_location(High, d) ||
               newp2[0].in_direction(d) > v.boundary_location(High, d)) {
        for (int j=0; j<ncopies; j++) {
          newp1[ncopies+j] = newp1[j] - lattice_vector(d);
          newp2[ncopies+j] = newp2[j] - lattice_vector(d);
          kphase[ncopies+j] = kphase[j] * eikna[d];
        }
        ncopies *= 2;
      }
    }
}

bool fields::locate_component_point(component *c, ivec *there,
                                    complex<double> *phase) const {
  // returns true if this point and component exist in the user_volume.  If
  // that is the case, on return *c and *there store the component and
  // location of where the point actually is, and *phase determines holds
  // the phase needed to get the true field.  If the point is not located,
  // *c and *there will hold undefined values.

  // Check if nothing tricky is needed...
  *phase = 1.0;
  if (!locate_point_in_user_volume(there, phase)) return false;
  // Check if a rotation or inversion brings the point in...
  if (user_volume.owns(*there))