structure.cpp

来自mit的fdtd开放性源代码meep

}
#endif

// fallback meaneps using libctl's adaptive cubature routine
void geom_epsilon::fallback_meaneps(double &meps, double &minveps,
				    const meep::geometric_volume &gv,
				    double tol, int maxeval)
{
  number esterr;
  integer errflag, n;
  number xmin[3], xmax[3];
  vector3 gvmin, gvmax;
  gvmin = vec_to_vector3(gv.get_min_corner());
  gvmax = vec_to_vector3(gv.get_max_corner());
  xmin[0] = gvmin.x; xmax[0] = gvmax.x; 
  if (dim == meep::Dcyl) {
    xmin[1] = gvmin.z; xmin[2] = gvmin.y; xmax[1] = gvmax.z; xmax[2] = gvmax.y;
  }
  else{
    xmin[1] = gvmin.y; xmin[2] = gvmin.z; xmax[1] = gvmax.y; xmax[2] = gvmax.z;
  }
  if (xmin[2] == xmax[2])
    n = xmin[1] == xmax[1] ? 1 : 2;
  else
    n = 3;
  double vol = 1;
  for (int i = 0; i < n; ++i) vol *= xmax[i] - xmin[i];
  if (dim == meep::Dcyl) vol *= (xmin[0] + xmax[0]) * 0.5;
  eps_ever_negative = 0;
#ifdef CTL_HAS_COMPLEX_INTEGRATION
  cnumber ret = cadaptive_integration(ceps_func, xmin, xmax, n, (void*) this,
				      0, tol, maxeval, &esterr, &errflag);
  meps = ret.re / vol;
  minveps = ret.im / vol;
#else
  meps = adaptive_integration(eps_func, xmin, xmax, n, (void*) this,
			      0, tol, maxeval, &esterr, &errflag) / vol;
  minveps = adaptive_integration(inveps_func, xmin, xmax, n, (void*) this,
				 0, tol, maxeval, &esterr, &errflag) / vol;
#endif
  if (eps_ever_negative) // averaging negative eps causes instability
    minveps = 1.0 / (meps = eps(gv.center()));
}

bool geom_epsilon::has_chi3()
{
  for (int i = 0; i < geometry.num_items; ++i) {
    if (geometry.items[i].material.which_subclass == MTS::MEDIUM) {
      if (geometry.items[i].material.subclass.medium_data->chi3 != 0)
	return true; 
    }
  }
    /* FIXME: what to do about material-functions?
       Currently, we require that at least *one* ordinary material
       property have non-zero chi3 for Kerr to be enabled.   It might
       be better to have set_chi3 automatically delete chi3[] if the
       chi3's are all zero. */
  return (default_material.which_subclass == MTS::MEDIUM &&
	  default_material.subclass.medium_data->chi3 != 0);
}

double geom_epsilon::chi3(const meep::vec &r) {
  material_type material;
  bool destroy_material = get_material_pt(material, r);
  
  double chi3_val;
  switch (material.which_subclass) {
  case MTS::MEDIUM:
    chi3_val = material.subclass.medium_data->chi3;
    break;
  default:
    chi3_val = 0;
  }
  
  if (destroy_material)
    material_type_destroy(material);
  
  return chi3_val;
}

bool geom_epsilon::has_chi2()
{
  for (int i = 0; i < geometry.num_items; ++i) {
    if (geometry.items[i].material.which_subclass == MTS::MEDIUM) {
      if (geometry.items[i].material.subclass.medium_data->chi2 != 0)
	return true; 
    }
  }
    /* FIXME: what to do about material-functions?
       Currently, we require that at least *one* ordinary material
       property have non-zero chi2 for Kerr to be enabled.   It might
       be better to have set_chi2 automatically delete chi2[] if the
       chi2's are all zero. */
  return (default_material.which_subclass == MTS::MEDIUM &&
	  default_material.subclass.medium_data->chi2 != 0);
}

double geom_epsilon::chi2(const meep::vec &r) {
  material_type material;
  bool destroy_material = get_material_pt(material, r);
  
  double chi2_val;
  switch (material.which_subclass) {
  case MTS::MEDIUM:
    chi2_val = material.subclass.medium_data->chi2;
    break;
  default:
    chi2_val = 0;
  }
  
  if (destroy_material)
    material_type_destroy(material);
  
  return chi2_val;
}

bool geom_epsilon::has_mu()
{
  for (int i = 0; i < geometry.num_items; ++i) {
    if (geometry.items[i].material.which_subclass == MTS::MEDIUM) {
      if (geometry.items[i].material.subclass.medium_data->mu != 1)
	return true; 
    }
  }
    /* FIXME: what to do about material-functions?
       Currently, we require that at least *one* ordinary material
       property have non-zero chi2 for Kerr to be enabled.   It might
       be better to have set_chi2 automatically delete chi2[] if the
       chi2's are all zero. */
  return (default_material.which_subclass == MTS::MEDIUM &&
	  default_material.subclass.medium_data->mu != 1);
}

double geom_epsilon::mu(const meep::vec &r) {
  material_type material;
  bool destroy_material = get_material_pt(material, r);

  double mu_val;
  switch (material.which_subclass) {
  case MTS::MEDIUM:
    mu_val = material.subclass.medium_data->mu;
    break;
  default:
    mu_val = 1;
  }
  
  if (destroy_material)
    material_type_destroy(material);
  
  return mu_val;
}

static double get_cnd(meep::component c, const medium *m) {
  switch (c) {
  case meep::Dr: case meep::Dx: return m->D_conductivity_diag.x;
  case meep::Dp: case meep::Dy: return m->D_conductivity_diag.y;
  case meep::Dz: return m->D_conductivity_diag.z;
  case meep::Br: case meep::Bx: return m->B_conductivity_diag.x;
  case meep::Bp: case meep::By: return m->B_conductivity_diag.y;
  case meep::Bz: return m->B_conductivity_diag.z;
  default: return 0;
  }
}

bool geom_epsilon::has_conductivity(meep::component c)
{
  for (int i = 0; i < geometry.num_items; ++i) {
    if (geometry.items[i].material.which_subclass == MTS::MEDIUM) {
      if (get_cnd(c, geometry.items[i].material.subclass.medium_data) != 0)
	return true; 
    }
  }
    /* FIXME: what to do about material-functions?
       Currently, we require that at least *one* ordinary material
       property have non-zero chi2 for Kerr to be enabled.   It might
       be better to have set_chi2 automatically delete chi2[] if the
       chi2's are all zero. */
  return (default_material.which_subclass == MTS::MEDIUM &&
	  get_cnd(c, default_material.subclass.medium_data) != 0);
}

double geom_epsilon::conductivity(meep::component c, const meep::vec &r) {
  material_type material;
  bool destroy_material = get_material_pt(material, r);

  double cond_val;
  switch (material.which_subclass) {
  case MTS::MEDIUM:
    cond_val = get_cnd(c, material.subclass.medium_data);
    break;
  default:
    cond_val = 0;
  }
  
  if (destroy_material)
    material_type_destroy(material);
  
  return cond_val;
}

double geom_epsilon::sigma(const meep::vec &r) {
  vector3 p = vec_to_vector3(r);

  boolean inobject;
  material_type material =
    material_of_unshifted_point_in_tree_inobject(p, restricted_tree, &inobject);
  
  int destroy_material = 0;
  if (material.which_subclass == MTS::MATERIAL_TYPE_SELF) {
    material = default_material;
  }
  if (material.which_subclass == MTS::MATERIAL_FUNCTION) {
    material = eval_material_func(material.subclass.
				  material_function_data->material_func,
				  p);
    destroy_material = 1;
  }
  
  double sigma = 0;
  if (material.which_subclass == MTS::MEDIUM) {
    polarizability_list plist = 
      material.subclass.medium_data->polarizations;
    for (int j = 0; j < plist.num_items; ++j)
      if (plist.items[j].omega == omega &&
	  plist.items[j].gamma == gamma &&
	  plist.items[j].delta_epsilon == deps &&
	  plist.items[j].energy_saturation == energy_sat) {
	sigma = plist.items[j].sigma;
	break;
      }
  }
  
  if (destroy_material)
    material_type_destroy(material);

  return sigma;
}

struct pol {
  double omega, gamma, deps, esat;
  struct pol *next;
};

// add a polarization to the list if it is not already there
static pol *add_pol(pol *pols,
		    double omega, double gamma, double deps, double esat)
{
  struct pol *p = pols;
  while (p && !(p->omega == omega && p->gamma == gamma
		&& p->deps == deps && p->esat == esat))
    p = p->next;
  if (!p) {
    p = new pol;
    p->omega = omega;
    p->gamma = gamma;
    p->deps = deps;
    p->esat = esat;
    p->next = pols;
    pols = p;
  }
  return pols;
}

static pol *add_pols(pol *pols, const polarizability_list plist) {
  for (int j = 0; j < plist.num_items; ++j) {
    pols = add_pol(pols,
		   plist.items[j].omega, plist.items[j].gamma,
		   plist.items[j].delta_epsilon,
		   plist.items[j].energy_saturation);
  }
  return pols;
}

void geom_epsilon::add_polarizabilities(meep::structure *s) {
  pol *pols = 0;

  // construct a list of the unique polarizabilities in the geometry:
  for (int i = 0; i < geometry.num_items; ++i) {
    if (geometry.items[i].material.which_subclass == MTS::MEDIUM)
      pols = add_pols(pols, geometry.items[i].material
		      .subclass.medium_data->polarizations);
  }
  if (default_material.which_subclass == MTS::MEDIUM)
    pols = add_pols(pols, default_material
		    .subclass.medium_data->polarizations);
    
  for (struct pol *p = pols; p; p = p->next) {
    master_printf("polarizability: omega=%g, gamma=%g, deps=%g, esat=%g\n",
		  p->omega, p->gamma, p->deps, p->esat);
    s->add_polarizability(*this, p->omega, p->gamma, p->deps, p->esat);
  }
  
  while (pols) {
    struct pol *p = pols;
    pols = pols->next;
    delete p;
  }
}

/***********************************************************************/

// wrapper around Scheme function for PML profile
static double scm_pml_profile(double u, void *f_)
{
  SCM f = (SCM) f_;
  return ctl_convert_number_to_c(gh_call1(f, ctl_convert_number_to_scm(u)));
}

// for passing to multidimensional integration routine
static double scm_pml_profile2(int dim, double *u, void *f_)
{
  SCM f = (SCM) f_; (void) dim;
  return ctl_convert_number_to_c(gh_call1(f, ctl_convert_number_to_scm(*u)));
}

meep::structure *make_structure(int dims, vector3 size, vector3 center,
				double resolution, bool enable_averaging,
				double subpixel_tol, int subpixel_maxeval,
				bool ensure_periodicity_p,
				geometric_object_list geometry,
				material_type default_mat,
				pml_list pml_layers,
				symmetry_list symmetries,
				int num_chunks, double Courant)
{
  master_printf("-----------\nInitializing structure...\n");
  
  // only cartesian lattices are currently allowed
  geom_initialize();
  geometry_center = center;
  
  number no_size = 2.0 / ctl_get_number("infinity");
  if (size.x <= no_size)
    size.x = 0.0;
  if (size.y <= no_size)
    size.y = 0.0;
  if (size.z <= no_size)
    size.z = 0.0;
  
  set_dimensions(dims);
  
  geometry_lattice.size = size;

  master_printf("Working in %s dimensions.\n", meep::dimension_name(dim));
  
  meep::volume v;
  switch (dims) {
  case 0: case 1:
    v = meep::vol1d(size.z, resolution);
    break;
  case 2:
    v = meep::vol2d(size.x, size.y, resolution);
    break;
  case 3:
    v = meep::vol3d(size.x, size.y, size.z, resolution);
    break;
  case CYLINDRICAL:
    v = meep::volcyl(size.x, size.z, resolution);
    break;
  default:
    CK(0, "unsupported dimensionality");
  }
  v.center_origin();
  v.shift_origin(vector3_to_vec(center));
  
  meep::symmetry S;
  for (int i = 0; i < symmetries.num_items; ++i) 
    switch (symmetries.items[i].which_subclass) {
    case symmetry::SYMMETRY_SELF: break; // identity
    case symmetry::MIRROR_SYM:
      S = S + meep::mirror(meep::direction(symmetries.items[i].direction), v)
	* complex<double>(symmetries.items[i].phase.re,
			  symmetries.items[i].phase.im);
      break;
    case symmetry::ROTATE2_SYM:
      S = S + meep::rotate2(meep::direction(symmetries.items[i].direction), v)
	* complex<double>(symmetries.items[i].phase.re,
			  symmetries.items[i].phase.im);
      break;
    case symmetry::ROTATE4_SYM:
      S = S + meep::rotate4(meep::direction(symmetries.items[i].direction), v)
	* complex<double>(symmetries.items[i].phase.re,
			  symmetries.items[i].phase.im);
      break;
    }

  meep::boundary_region br;
  for (int i = 0; i < pml_layers.num_items; ++i) {
    double umin = 0, umax = 1, esterr;
    int errflag;
    using namespace meep;
    if (pml_layers.items[i].direction == -1) {
      LOOP_OVER_DIRECTIONS(v.dim, d) {
	if (pml_layers.items[i].side == -1) {
	  FOR_SIDES(b)
	    br = br + meep::boundary_region
	    (meep::boundary_region::PML,
	     pml_layers.items[i].thickness,
	     pow(pml_layers.items[i].R_asymptotic,
		 pml_layers.items[i].strength),
	     scm_pml_profile, pml_layers.items[i].pml_profile,
	     adaptive_integration(scm_pml_profile2, &umin, &umax, 1,
				  (void*) pml_layers.items[i].pml_profile,
				  1e-9, 1e-4, 50000, &esterr, &errflag),
	     d, b);
	}
	else
	  br = br + meep::boundary_region
	    (meep::boundary_region::PML,
	     pml_layers.items[i].thickness,
	     pow(pml_layers.items[i].R_asymptotic,
		 pml_layers.items[i].strength),
	     scm_pml_profile, pml_layers.items[i].pml_profile,
	     adaptive_integration(scm_pml_profile2, &umin, &umax, 1,
				  (void*) pml_layers.items[i].pml_profile,
				  1e-9, 1e-4, 50000, &esterr, &errflag),
	     d,
	     (meep::boundary_side) pml_layers.items[i].side);
      }
    }
    else {
	if (pml_layers.items[i].side == -1) {
	  FOR_SIDES(b)
	    br = br + meep::boundary_region
	    (meep::boundary_region::PML,
	     pml_layers.items[i].thickness,
	     pow(pml_layers.items[i].R_asymptotic,
		 pml_layers.items[i].strength),
	     scm_pml_profile, pml_layers.items[i].pml_profile,
	     adaptive_integration(scm_pml_profile2, &umin, &umax, 1,
				  (void*) pml_layers.items[i].pml_profile,
				  1e-9, 1e-4, 50000, &esterr, &errflag),
	     (meep::direction) pml_layers.items[i].direction,
	     b);
	}
	else
	  br = br + meep::boundary_region
	    (meep::boundary_region::PML,
	     pml_layers.items[i].thickness,
	     pow(pml_layers.items[i].R_asymptotic,
		 pml_layers.items[i].strength),
	     scm_pml_profile, pml_layers.items[i].pml_profile,
	     adaptive_integration(scm_pml_profile2, &umin, &umax, 1,
				  (void*) pml_layers.items[i].pml_profile,
				  1e-9, 1e-4, 50000, &esterr, &errflag),
	     (meep::direction) pml_layers.items[i].direction,
	     (meep::boundary_side) pml_layers.items[i].side);
    }
  }
  
  ensure_periodicity = ensure_periodicity_p;
  default_material = default_mat;
  geom_epsilon geps(geometry, v.pad().surroundings());

  if (subpixel_maxeval < 0) subpixel_maxeval = 0; // no limit

  meep::structure *s = new meep::structure(v, geps, br, S, 
					   num_chunks, Courant,
					   enable_averaging,
					   subpixel_tol,
					   subpixel_maxeval);

  geps.add_polarizabilities(s);

  master_printf("-----------\n");
  
  return s;
}

/*************************************************************************/