meep.hpp

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

  void use_pml(direction d, boundary_side b, double dx, double strength);
  void add_to_effort_volumes(const volume &new_effort_volume, 
			     double extra_effort);
  void choose_chunkdivision(const volume &v, int num_chunks,
			    const boundary_region &br, const symmetry &s);
  void check_chunks();
  void changing_chunks();
};

class src_vol;
class bandsdata;
class fields;
class fields_chunk;
class flux_vol;

// Time-dependence of a current source, intended to be overridden by
// subclasses.  current() and dipole() are be related by
// current = d(dipole)/dt (or rather, the finite-difference equivalent).
class src_time {
 public:
  src_time() { current_time = nan; current_current = 0.0; next = NULL; }
  virtual ~src_time() { delete next; }
  src_time(const src_time &t) { 
       current_time = t.current_time;
       current_current = t.current_current;
       current_dipole = t.current_dipole;
       if (t.next) next = t.next->clone(); else next = NULL;
  }
  
  complex<double> dipole() const { return current_dipole; }
  complex<double> current() const { return current_current; }
  void update(double time, double dt) {
    if (time != current_time) {
      current_dipole = dipole(time);
      current_current = current(time, dt);
      current_time = time;
    }
  }

  complex<double> current(double time, double dt) const { 
    return ((dipole(time + dt) - dipole(time)) / dt);
  }

  double last_time_max() { return last_time_max(0.0); }
  double last_time_max(double after);
  
  src_time *add_to(src_time *others, src_time **added) const;
  src_time *next;

  // subclasses should override these methods:
  virtual complex<double> dipole(double time) const { (void)time; return 0; }
  virtual double last_time() const { return 0.0; }
  virtual src_time *clone() const { return new src_time(*this); }
  virtual bool is_equal(const src_time &t) const { (void)t; return 1; }
  virtual complex<double> frequency() const { return 0.0; }

 private:
  double current_time;
  complex<double> current_dipole, current_current;
};

bool src_times_equal(const src_time &t1, const src_time &t2);

// Gaussian-envelope source with given frequency, width, peak-time, cutoff
class gaussian_src_time : public src_time {
 public:
  gaussian_src_time(double f, double fwidth, double s = 5.0);
  gaussian_src_time(double f, double w, double start_time, double end_time);
  virtual ~gaussian_src_time() {}

  virtual complex<double> dipole(double time) const;
  virtual double last_time() const { return peak_time + cutoff; };
  virtual src_time *clone() const { return new gaussian_src_time(*this); }
  virtual bool is_equal(const src_time &t) const;
  virtual complex<double> frequency() const { return freq; }

 private:
  double freq, width, peak_time, cutoff;
};

// Continuous (CW) source with (optional) slow turn-on and/or turn-off.
class continuous_src_time : public src_time {
 public:
  continuous_src_time(complex<double> f, double w = 0.0, 
		      double st = 0.0, double et = infinity,
		      double s = 3.0) : freq(f), width(w), start_time(st),
					end_time(et), slowness(s) {}
  virtual ~continuous_src_time() {}
  
  virtual complex<double> dipole(double time) const;
  virtual double last_time() const { return end_time; };
  virtual src_time *clone() const { return new continuous_src_time(*this); }
  virtual bool is_equal(const src_time &t) const;
  virtual complex<double> frequency() const { return freq; }
  
 private:
  complex<double> freq;
  double width, start_time, end_time, slowness;
};

// user-specified source function with start and end times
class custom_src_time : public src_time {
 public:
  custom_src_time(complex<double> (*func)(double t, void *), void *data,
		  double st = -infinity, double et = infinity)
    : func(func), data(data), start_time(st), end_time(et) {}
  virtual ~custom_src_time() {}
  
  virtual complex<double> dipole(double time) const { return func(time,data); }
  virtual double last_time() const { return end_time; };
  virtual src_time *clone() const { return new custom_src_time(*this); }
  virtual bool is_equal(const src_time &t) const;
  
 private:
  complex<double> (*func)(double t, void *);
  void *data;
  double start_time, end_time;
};

class monitor_point {
 public:
  monitor_point();
  ~monitor_point();
  vec loc;
  double t;
  complex<double> f[NUM_FIELD_COMPONENTS];
  monitor_point *next;

  complex<double> get_component(component);
  double poynting_in_direction(direction d);
  double poynting_in_direction(vec direction_v);

  // When called with only its first four arguments, fourier_transform
  // performs an FFT on its monitor points, putting the frequencies in f
  // and the amplitudes in a.  Yes, the frequencies are trivial and
  // redundant, but this saves you the risk of making a mistake in
  // converting your units.  Note also, that in this case f is always a
  // real number, although it's stored in a complex.
  //
  // Note that in either case, fourier_transform assumes that the monitor
  // points are all equally spaced in time.
  void fourier_transform(component w,
                         complex<double> **a, complex<double> **f, int *numout,
                         double fmin=0.0, double fmax=0.0, int maxbands=100);
  // harminv works much like fourier_transform, except that it is not yet
  // implemented.
  void harminv(component w,
               complex<double> **a, complex<double> **f,
               int *numout, double fmin, double fmax,
               int maxbands);
};

// dft.cpp
// this should normally only be created with fields::add_dft
class dft_chunk {
public:
  dft_chunk(fields_chunk *fc_,
	    ivec is_, ivec ie_,
	    vec s0_, vec s1_, vec e0_, vec e1_,
	    double dV0_, double dV1_,
	    complex<double> scale_,
	    component c_,
	    const void *data_);
  ~dft_chunk();
  
  void update_dft(double time);

  void scale_dft(complex<double> scale);

  void operator-=(const dft_chunk &chunk);

  // the frequencies to loop_in_chunks
  double omega_min, domega;
  int Nomega;

  component c; // component to DFT (possibly transformed by symmetry)

  int N; // number of spatial points (on epsilon grid)
  complex<double> *dft; // N x Nomega array of DFT values.

  struct dft_chunk *next_in_chunk; // per-fields_chunk list of DFT chunks
  struct dft_chunk *next_in_dft; // next for this particular DFT vol./component

private:
  // parameters passed from field_integrate:
  fields_chunk *fc;
  ivec is, ie;
  vec s0, s1, e0, e1;
  double dV0, dV1;
  complex<double> scale; // scale factor * phase from shift and symmetry

  // cache of exp(iwt) * scale, of length Nomega
  complex<double> *dft_phase;

  int avg1, avg2; // index offsets for average to get epsilon grid
};

void save_dft_hdf5(dft_chunk *dft_chunks, component c, h5file *file,
		   const char *dprefix = 0);
void load_dft_hdf5(dft_chunk *dft_chunks, component c, h5file *file,
		   const char *dprefix = 0);

// dft.cpp (normally created with fields::add_dft_flux)
class dft_flux {
public:
  dft_flux(const component cE_, const component cH_,
	   dft_chunk *E_, dft_chunk *H_,
	   double fmin, double fmax, int Nf);
  dft_flux(const dft_flux &f);

  double *flux();

  void save_hdf5(h5file *file, const char *dprefix = 0);
  void load_hdf5(h5file *file, const char *dprefix = 0);

  void operator-=(const dft_flux &fl) { if (E && fl.E) *E -= *fl.E; if (H && fl.H) *H -= *fl.H; }

  void save_hdf5(fields &f, const char *fname, const char *dprefix = 0,
		 const char *prefix = 0);
  void load_hdf5(fields &f, const char *fname, const char *dprefix = 0,
		 const char *prefix = 0);

  void scale_dfts(complex<double> scale);

  void remove() { delete E; delete H; E = H = 0; }

  double freq_min, dfreq;
  int Nfreq;
  dft_chunk *E, *H;
  component cE, cH;
};

enum in_or_out { Incoming=0, Outgoing };
enum connect_phase { CONNECT_PHASE = 0, CONNECT_NEGATE=1, CONNECT_COPY=2 };

class fields_chunk {
 public:
  double *f[NUM_FIELD_COMPONENTS][2];
  double *f_backup[NUM_FIELD_COMPONENTS][2];
  double *f_p_pml[NUM_FIELD_COMPONENTS][2];
  double *f_m_pml[NUM_FIELD_COMPONENTS][2];
  double *f_backup_p_pml[NUM_FIELD_COMPONENTS][2];
  double *f_backup_m_pml[NUM_FIELD_COMPONENTS][2];

  bool have_d_minus_p;
  double *d_minus_p[NUM_FIELD_COMPONENTS][2];

  dft_chunk *dft_chunks;

  double **zeroes[NUM_FIELD_TYPES]; // Holds pointers to metal points.
  int num_zeroes[NUM_FIELD_TYPES];
  double **connections[NUM_FIELD_TYPES][CONNECT_COPY+1][Outgoing+1];
  int num_connections[NUM_FIELD_TYPES][CONNECT_COPY+1][Outgoing+1];
  complex<double> *connection_phases[NUM_FIELD_TYPES];

  polarization *pol, *olpol;
  double a, Courant, dt; // res. a, Courant num., and timestep dt=Courant/a
  volume v;
  geometric_volume gv;
  double m, rshift;
  int is_real;
  bandsdata *bands;
  src_vol *e_sources, *h_sources;
  const structure_chunk *new_s;
  structure_chunk *s;
  const char *outdir;

  fields_chunk(structure_chunk *, const char *outdir, double m);
  fields_chunk(const fields_chunk &);
  ~fields_chunk();

  // step.cpp
  double peek_field(component, const vec &);

  void use_real_fields();
  bool have_component(component c, bool is_complex = false) {
    switch (c) {
    case Dielectric:
      return !is_complex;
    default:
      return (f[c][0] && f[c][is_complex]);
    }
  }

  double last_source_time();
  // monitor.cpp
  complex<double> get_field(component, const ivec &) const;

  // for non-collective interpolation:
  geometric_volume get_field_gv(component) const;
  complex<double> get_field(component, const vec &) const;

  complex<double> get_polarization_field(const polarizability_identifier &p,
                                         component c, const ivec &iloc) const;
  double get_polarization_energy(const ivec &) const;
  double my_polarization_energy(const ivec &) const;
  double get_polarization_energy(const polarizability_identifier &, const ivec &) const;
  double my_polarization_energy(const polarizability_identifier &, const ivec &) const;
  double get_inveps(component, direction, const ivec &iloc) const;
  double get_eps(const ivec &iloc) const;
  complex<double> analytic_epsilon(double freq, const vec &) const;
  
  void backup_h();
  void restore_h();

  void set_output_directory(const char *name);
  void verbose(int v=1) { verbosity = v; }

  double count_volume(component);
  friend class fields;

  int n_proc() const { return s->n_proc(); };
  int is_mine() const { return s->is_mine(); };
  // boundaries.cpp
  void zero_metal(field_type);
  // polarization.cpp
  void initialize_polarization_energy(const polarizability_identifier &,
                                      double energy(const vec &));

  // fields.cpp
  void alloc_f(component c);
  void remove_sources();
  void remove_polarizabilities();
  void zero_fields();

 private: 
  int verbosity; // Turn on verbosity for debugging purposes...
  // fields.cpp
  void figure_out_step_plan();
  bool have_plus_deriv[NUM_FIELD_COMPONENTS], have_minus_deriv[NUM_FIELD_COMPONENTS];
  component plus_component[NUM_FIELD_COMPONENTS], minus_component[NUM_FIELD_COMPONENTS];
  direction plus_deriv_direction[NUM_FIELD_COMPONENTS],
            minus_deriv_direction[NUM_FIELD_COMPONENTS];
  int num_each_direction[3], stride_each_direction[3];
  int num_any_direction[5], stride_any_direction[5];
  // bands.cpp
  void record_bands(int tcount);
  // step.cpp
  void phase_in_material(const structure_chunk *s);
  void phase_material(int phasein_time);
  void step_h();
  void step_h_source(src_vol *);
  void step_d();
  void step_d_source(src_vol *);
  void update_e_from_d();
  void update_from_e();
  void calc_sources(double time);

  // initialize.cpp
  void initialize_field(component, complex<double> f(const vec &));
  void initialize_polarizations(polarization *op=NULL, polarization *np=NULL);
  void initialize_with_nth_te(int n, double kz);
  void initialize_with_nth_tm(int n, double kz);
  // boundaries.cpp
  void alloc_extra_connections(field_type, connect_phase, in_or_out, int);
  // dft.cpp
  void update_dfts(double timeE, double timeH);

  void changing_structure();
};

enum boundary_condition { Periodic=0, Metallic, Magnetic, None };
enum time_sink { Connecting, Stepping, Boundaries, MpiTime,
                 FieldOutput, FourierTransforming, Other };

typedef void (*field_chunkloop)(fields_chunk *fc, int ichunk, component cgrid,
				ivec is, ivec ie,
				vec s0, vec s1, vec e0, vec e1,
				double dV0, double dV1,
				ivec shift, complex<double> shift_phase,