structure.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 <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

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

namespace meep {

structure::structure()
  : Courant(0.5), gv(D1) // Aaack, this is very hokey.
{
  num_chunks = 0;
  num_effort_volumes = 0;
  effort_volumes = NULL;
  effort = NULL;
  outdir = ".";
  S = identity();
  a = 1; dt = Courant/a;
}

typedef structure_chunk *structure_chunk_ptr;

structure::structure(const volume &thev, material_function &eps,
		     const boundary_region &br,
		     const symmetry &s,
		     int num, double Courant, bool use_anisotropic_averaging,
		     double tol, int maxeval) :
  Courant(Courant), gv(D1) // Aaack, this is very hokey.
{
  outdir = ".";
  choose_chunkdivision(thev, num == 0 ? count_processors() : num, br, s);
  set_materials(eps, use_anisotropic_averaging, tol, maxeval);
}

structure::structure(const volume &thev, double eps(const vec &), 
		     const boundary_region &br,
		     const symmetry &s,
		     int num, double Courant, bool use_anisotropic_averaging,
		     double tol, int maxeval) :
  Courant(Courant), gv(D1) // Aaack, this is very hokey.
{
  outdir = ".";
  choose_chunkdivision(thev, num == 0 ? count_processors() : num, br, s);
  simple_material_function epsilon(eps);
  set_materials(epsilon, use_anisotropic_averaging, tol, maxeval);
}

void structure::choose_chunkdivision(const volume &thev, 
				     int desired_num_chunks, 
				     const boundary_region &br,
				     const symmetry &s) {
  user_volume = thev;
  if (thev.dim == Dcyl && thev.get_origin().r() < 0)
    abort("r < 0 origins are not supported");
  v = thev;
  gv = v.surroundings();
  S = s;
  a = v.a;
  dt = Courant/a;

  // First, reduce overall volume v by symmetries:
  if (S.multiplicity() > 1) {
    bool break_this[3];
    for (int dd=0;dd<3;dd++) {
      const direction d = (direction) dd;
      break_this[dd] = false;
      for (int n=0;n<S.multiplicity();n++)
        if (has_direction(thev.dim,d) &&
            (S.transform(d,n).d != d || S.transform(d,n).flipped)) {
          if (thev.num_direction(d) & 1 && !break_this[d])
            master_printf("Padding %s to even number of grid points.\n",
			  direction_name(d));
          break_this[dd] = true;
        }
    }
    int break_mult = 1;
    for (int d=0;d<3;d++) {
      if (break_mult == S.multiplicity()) break_this[d] = false;
      if (break_this[d]) {
	break_mult *= 2;
	master_printf("Halving computational cell along direction %s\n",
		      direction_name(direction(d)));
	v = v.halve((direction)d);
      }
    }
    // Before padding, find the corresponding geometric volume.
    gv = v.surroundings();
    // Pad the little cell in any direction that we've shrunk:
    for (int d=0;d<3;d++)
      if (break_this[d]) v = v.pad((direction)d);
  }

  // initialize effort volumes
  num_effort_volumes = 1;
  effort_volumes = new volume[num_effort_volumes];
  effort_volumes[0] = v;
  effort = new double[num_effort_volumes];
  effort[0] = 1.0;

  // Next, add effort volumes for PML boundary regions:
  br.apply(this);

  // Finally, create the chunks:
  num_chunks = 0;
  chunks = new structure_chunk_ptr[desired_num_chunks * num_effort_volumes];
  for (int i = 0; i < desired_num_chunks; i++) {
    const int proc = i * count_processors() / desired_num_chunks;
    volume vi = v.split_by_effort(desired_num_chunks, i,
                                  num_effort_volumes, effort_volumes, effort);
    for (int j = 0; j < num_effort_volumes; j++) {
      volume vc;
      if (vi.intersect_with(effort_volumes[j], &vc)) {
	chunks[num_chunks] = new structure_chunk(vc, gv, Courant, proc);
	br.apply(this, chunks[num_chunks++]);
      }
    }
  }

  check_chunks();
}

void boundary_region::apply(structure *s) const {
  if (has_direction(s->v.dim, d) && s->user_volume.has_boundary(side, d)
      && s->user_volume.num_direction(d) > 1) {
    switch (kind) {
    case NOTHING_SPECIAL: break;
    case PML: s->use_pml(d, side, thickness, strength); break;
    default: abort("unknown boundary region kind");
    }
  }
  if (next)
    next->apply(s);
}

void boundary_region::apply(const structure *s, structure_chunk *sc) const {
  if (has_direction(s->v.dim, d) && s->user_volume.has_boundary(side, d)
      && s->user_volume.num_direction(d) > 1) {
    switch (kind) {
    case NOTHING_SPECIAL: break;
    case PML: 
      sc->use_pml(d, thickness, s->user_volume.boundary_location(side, d),
		  strength); 
      break;
    default: abort("unknown boundary region kind");
    }
  }
  if (next)
    next->apply(s, sc);
}

boundary_region pml(double thickness, direction d, boundary_side side) {
  return boundary_region(boundary_region::PML, thickness, 1.0, d, side, NULL);
}
boundary_region pml(double thickness, direction d) {
  return (pml(thickness, d, Low) + pml(thickness, d, High));
}
boundary_region pml(double thickness) {
  boundary_region r;
  for (int id = 0; id < 5; ++id)
    r = r + pml(thickness, (direction) id);
  return r;
}

// First check that the chunk volumes do not intersect and that they add
// up to the total volume
void structure::check_chunks() {
  volume vol_intersection;
  for (int i=0; i<num_chunks; i++)
    for (int j=i+1; j<num_chunks; j++)
      if (chunks[i]->v.intersect_with(chunks[j]->v, &vol_intersection))
        abort("chunks[%d] intersects with chunks[%d]\n", i, j);
  // FIXME: should use 'long long' else will fail if grid > 2e9 points
  int sum = 0;
  for (int i=0; i<num_chunks; i++) {
    int grid_points = 1;
    LOOP_OVER_DIRECTIONS(chunks[i]->v.dim, d)
      grid_points *= chunks[i]->v.num_direction(d);
    sum += grid_points;
  }
  int v_grid_points = 1;
  LOOP_OVER_DIRECTIONS(v.dim, d) v_grid_points *= v.num_direction(d);
  if (sum != v_grid_points)
    abort("v_grid_points = %d, sum(chunks) = %d\n", v_grid_points, sum);
}

void structure::add_to_effort_volumes(const volume &new_effort_volume,
                                double extra_effort) {
  volume *temp_volumes =
    new volume[(2*number_of_directions(v.dim)+1)*num_effort_volumes]; 
  double *temp_effort =
    new double[(2*number_of_directions(v.dim)+1)*num_effort_volumes];
  // Intersect previous mat_volumes with this new_effort_volume
  int counter = 0;
  for (int j=0; j<num_effort_volumes; j++) {
    volume intersection, others[6];
    int num_others;
    if (effort_volumes[j].intersect_with(new_effort_volume, &intersection,
                                         others, &num_others)) {
      if (num_others > 1) {
        printf("effort_volumes[%d]  ", j);
        effort_volumes[j].print();
        printf("new_effort_volume  ");
        new_effort_volume.print();
        // NOTE: this may not be a bug if this function is used for
        // something other than PML.
        abort("Did not expect num_others > 1 in add_to_effort_volumes\n");
      }
      temp_effort[counter] = extra_effort + effort[j];
      temp_volumes[counter] = intersection;
      counter++;
      for (int k = 0; k<num_others; k++) {
        temp_effort[counter] = effort[j];	  
        temp_volumes[counter] = others[k];
        counter++;
      }
    } else {
      temp_effort[counter] = effort[j];	  
      temp_volumes[counter] = effort_volumes[j];
      counter++;
    }
  }

  delete[] effort_volumes;
  delete[] effort;
  effort_volumes = temp_volumes;
  effort = temp_effort;
  num_effort_volumes = counter;
}

structure::structure(const structure *s) : gv(s->gv) {
  num_chunks = s->num_chunks;
  outdir = s->outdir;
  v = s->v;
  S = s->S;
  user_volume = s->user_volume;
  chunks = new structure_chunk_ptr[num_chunks];
  for (int i=0;i<num_chunks;i++)
    chunks[i] = new structure_chunk(s->chunks[i]);
  num_effort_volumes = s->num_effort_volumes;
  effort_volumes = new volume[num_effort_volumes];
  effort = new double[num_effort_volumes];
  for (int i=0;i<num_effort_volumes;i++) {
    effort_volumes[i] = s->effort_volumes[i];
    effort[i] = s->effort[i];
  }
  a = s->a;
  Courant = s->Courant;
  dt = s->dt;
}

structure::structure(const structure &s) : gv(s.gv) {
  num_chunks = s.num_chunks;
  outdir = s.outdir;
  v = s.v;
  S = s.S;
  user_volume = s.user_volume;
  chunks = new structure_chunk_ptr[num_chunks];
  for (int i=0;i<num_chunks;i++) {
    chunks[i] = new structure_chunk(s.chunks[i]);
  }
  num_effort_volumes = s.num_effort_volumes;
  effort_volumes = new volume[num_effort_volumes];
  effort = new double[num_effort_volumes];
  for (int i=0;i<num_effort_volumes;i++) {
    effort_volumes[i] = s.effort_volumes[i];
    effort[i] = s.effort[i];
  }  
  a = s.a;
  Courant = s.Courant;
  dt = s.dt;
}

structure::~structure() {
  for (int i=0;i<num_chunks;i++) {
    if (chunks[i]->refcount-- <= 1) delete chunks[i];
    chunks[i] = NULL; // Just to be sure...
  }
  delete[] chunks;
  delete[] effort_volumes;
  delete[] effort;
}

/* To save memory, the structure chunks are shared with the
   fields_chunk objects instead of making a copy.  However, to
   preserve the illusion that the structure and fields are
   independent objects, we implement copy-on-write semantics. */
void structure::changing_chunks() { // call this whenever chunks are modified
  for (int i=0; i<num_chunks; i++)
    if (chunks[i]->refcount > 1) { // this chunk is shared, so make a copy
      chunks[i]->refcount--;
      chunks[i] = new structure_chunk(chunks[i]);
    }
}

void structure::set_materials(material_function &mat, 
			      bool use_anisotropic_averaging,
			      double tol, int maxeval) {
  set_epsilon(mat, use_anisotropic_averaging, tol, maxeval);
  if (mat.has_chi3()) set_chi3(mat);
  if (mat.has_chi2()) set_chi2(mat);
}

void structure::set_epsilon(material_function &eps, 
			    bool use_anisotropic_averaging,
			    double tol, int maxeval) {
  double tstart = wall_time();
  changing_chunks();
  for (int i=0;i<num_chunks;i++)
    if (chunks[i]->is_mine())
      chunks[i]->set_epsilon(eps, use_anisotropic_averaging, tol, maxeval);
  if (!quiet)
    master_printf("time for set_epsilon = %g s\n", wall_time() - tstart);
}

void structure::set_epsilon(double eps(const vec &),
                            bool use_anisotropic_averaging,
			    double tol, int maxeval) {
  simple_material_function epsilon(eps);
  set_epsilon(epsilon, use_anisotropic_averaging, tol, maxeval);
}

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