slices.cpp

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

  const int buflen = 1024;
  char nname[buflen];
  if (*name) snprintf(nname, buflen, "%s-", name);
  else *nname = 0; // No additional name!
  char *n = new char[buflen];
  if (!n) abort("Allocation failure!\n");
  char time_step_string[buflen];
  snprintf(time_step_string, buflen, "%09.2f", time());
  
  snprintf(n, buflen, "%s/%senergy-%s.eps", outdir, nname, time_step_string);
  file *out = everyone_open_write(n);
  if (!out) abort("Unable to open file '%s' for slice output.\n", n);
  const double fmax = max(maxpolenergy_to_master(), -minpolenergy_to_master());
  if (am_master())
    output_complex_eps_header(v.eps_component(), fmax, user_volume,
                              what, out, n, v.eps_component());
  if (v.dim != D1) abort("Still only works in 1D.  :( \n");
  {
    bufprint buf(out);
    for (int i=0;i<num_chunks;i++)
      if (chunks[i]->is_mine())
        for (int sn=0;sn<S.multiplicity();sn++)
          for (int n=0;n<chunks[i]->v.ntot();n++) {
            const ivec here = chunks[i]->v.iloc(v.eps_component(), n);
            buf.printf("%g\t0\t%g\tP\n", chunks[i]->v[S.transform(here,sn)].z(),
                       chunks[i]->get_polarization_energy(p, here));
          }
  }
  all_wait();
  for (int i=0;i<num_chunks;i++)
    if (chunks[i]->is_mine())
      for (int sn=0;sn<S.multiplicity();sn++)
        for (int otherc=0;otherc<10;otherc++)
          if (S.transform((component)otherc,sn) == 0.5) // fixme -- what is this?
            eps_outline(v.eps_component(), chunks[i]->s->eps,
                        chunks[i]->v, what, S, sn, out);
  all_wait();
  if (am_master()) output_complex_eps_tail(out);
  everyone_close(out);
  delete[] n;
  finished_working();
}

void fields::eps_polarization_slice(const polarizability_identifier &p, const char *name) {
  if (v.dim == D3) abort("FIXME need to support 3D polarization slices...\n");
  geometric_volume what = user_volume.surroundings();
  if (v.dim == Dcyl) what.set_direction_min(R,-what.in_direction_max(R));
  eps_polarization_slice(p, what,name);
}

void fields::eps_polarization_slice(const polarizability_identifier &p,
                                    const geometric_volume &what, const char *name) {
  am_now_working_on(FieldOutput);
  const int buflen = 1024;
  char nname[buflen];
  if (*name) snprintf(nname, buflen, "%s-", name);
  else *nname = 0; // No additional name!
  char *n = new char[buflen];
  if (!n) abort("Allocation failure!\n");
  char time_step_string[buflen];
  snprintf(time_step_string, buflen, "%09.2f", time());

  FOR_ELECTRIC_COMPONENTS(c)
    if (v.has_field(c)) {
      snprintf(n, buflen, "%s/%sP-%s-%s.eps", outdir, nname,
               component_name(c), time_step_string);
      const double fmax = maxfieldmag_to_master(c);
      file *out = everyone_open_write(n);
      if (!out) {
        printf("Unable to open file '%s' for slice output.\n", n);
        return;
      }
      if (am_master())
        output_complex_eps_header(c, fmax, user_volume, what,
                                  out, n, v.eps_component());
      complex<double> phshift = optimal_phase_shift(c);
      all_wait();
      for (int i=0;i<num_chunks;i++)
        if (chunks[i]->is_mine())
          for (int sn=0;sn<S.multiplicity();sn++)
            FOR_COMPONENTS(otherc)
              if (S.transform(otherc,sn) == c)
                chunks[i]->output_eps_body(p, otherc,
                                           S, sn, what, out, phshift);
      all_wait();
      for (int i=0;i<num_chunks;i++)
        if (chunks[i]->is_mine())
          for (int sn=0;sn<S.multiplicity();sn++)
            FOR_COMPONENTS(otherc)
              if (S.transform(otherc,sn) == c)
                eps_outline(v.eps_component(), chunks[i]->s->eps,
                            chunks[i]->v, what, S, sn, out);
      all_wait();
      outline_chunks(out);
      all_wait();
      if (am_master()) output_complex_eps_tail(out);
      everyone_close(out);
    }
  delete[] n;
  finished_working();
}

void fields::eps_envelope(const char *name) {
  if (v.dim != D1) abort("Envelope only supported in 1D so far.\n");
  if (v.dim == D3) abort("Specify directions for EPS slices in 3D\n");
  geometric_volume what = user_volume.surroundings();
  if (v.dim == Dcyl) what.set_direction_min(R,-what.in_direction_max(R));
  eps_envelope(what,name);
}

void fields::eps_envelope(const geometric_volume &what, const char *name) {
  am_now_working_on(FieldOutput);
  const int buflen = 1024;
  char nname[buflen];
  if (*name) snprintf(nname, buflen, "%s-", name);
  else *nname = 0; // No additional name!
  char *n = (char *)malloc(buflen);
  if (!n) abort("Allocation failure!\n");
  char time_step_string[buflen];
  snprintf(time_step_string, buflen, "%09.2f", time());
  FOR_COMPONENTS(c)
    if (v.has_field(c)) {
      snprintf(n, buflen, "%s/%s%s-%s.eps", outdir, nname,
               component_name((component)c), time_step_string);
      const double fmax = maxfieldmag_to_master(c);
      file *out = everyone_open_write(n);
      if (!out) {
        printf("Unable to open file '%s' for slice output.\n", n);
        return;
      }
      if (am_master())
        output_complex_eps_header(c, fmax,
                                  user_volume, what,
                                  out, n, v.eps_component());
      for (double z = 0.0 + v.inva; z < user_volume.ntot(); z += v.inva)
        if (what.contains(vec(z))) {
          const double fhere = real(get_field(c, vec(z)));
          if (fhere > 0.0 && fhere > real(get_field(c, vec(z)-v.inva)) &&
              fhere > real(get_field(c, vec(z)+v.inva)))
            master_fprintf(out, "%g\t0\t%g\tP\n", z, fhere);
        }
      all_wait();
      for (int i=0;i<num_chunks;i++)
        if (chunks[i]->is_mine())
          for (int sn=0;sn<S.multiplicity();sn++)
            FOR_COMPONENTS(otherc)
              if (S.transform(otherc,sn) == c)
                eps_outline(v.eps_component(), chunks[i]->s->eps,
                            chunks[i]->v, what, S, sn, out);
      all_wait();
      outline_chunks(out);
      all_wait();
      if (am_master()) output_complex_eps_tail(out);
      everyone_close(out);
    }
  free(n);
  finished_working();
}

void fields::eps_slices(const char *name) {
  if (v.dim == D3) abort("Specify directions for EPS slices in 3D\n");
  geometric_volume what = user_volume.surroundings();
  if (v.dim == Dcyl) what.set_direction_min(R,-what.in_direction_max(R));
  eps_slices(what,name);
}

void fields::eps_slices(const vec &origin, const vec &xside, const vec &yside,
                        const double dx, const char *name) {
  am_now_working_on(FieldOutput);
  const int buflen = 1024;
  char nname[buflen];
  if (*name) snprintf(nname, buflen, "%s-", name);
  else *nname = 0; // No additional name!
  char *n = (char *)malloc(buflen);
  if (!n) abort("Allocation failure!\n");
  char time_step_string[buflen];
  snprintf(time_step_string, buflen, "%09.2f", time());
  // Define some convenient grid-related values:
  const vec xhat = xside / abs(xside);
  const vec yhat = yside / abs(yside);
  const double xlen = abs(xside);
  const double ylen = abs(yside);
  const double xmin = origin & xhat;
  const double ymin = origin & yhat;
  // Now start outputing pretty pictures.
  FOR_COMPONENTS(c)
    if (v.has_field(c)) {
      snprintf(n, buflen, "%s/%s%s-%s.eps", outdir, nname,
               component_name(c), time_step_string);
      const double fmax = maxfieldmag_to_master(c);
      file *out = everyone_open_write(n);
      if (!out) {
        printf("Unable to open file '%s' for slice output.\n", n);
        return;
      }
      if (am_master())
        output_eps_header(fmax, dx, xmin, xmin + xlen, ymin, ymin + ylen, out, n);
      complex<double> phshift = optimal_phase_shift(c);
      for (double x = xmin; x <= xmin + xlen + dx; x += dx)
        for (double y = ymin; y <= ymin + ylen + dx; y += dx)
          master_fprintf(out, "%g\t%g\t%g\tP\n", x, y,
                         real(phshift*
                              get_field(c, origin + xhat*(x-xmin) + yhat*(y-ymin))));
      for (double x = xmin; x <= xmin + xlen + dx; x += 1.0/v.a)
        for (double y = ymin; y <= ymin + ylen + dx; y += 1.0/v.a) {
          vec loc = origin + xhat*(x-xmin) + yhat*(y-ymin);
          if (has_eps_interface(&loc))
            master_fprintf(out, "%g\t%g\tB\n", loc & xhat, loc & yhat);
        }
      if (am_master()) output_complex_eps_tail(out);
      everyone_close(out);
    }
  free(n);
  finished_working();
}

bool fields::has_eps_interface(vec *loc) const {
  ivec ilocs[8];
  double w[8];
  double val[8];
  for (int i=0;i<8;i++) val[i] = 0.0;
  v.interpolate(v.eps_component(), *loc, ilocs, w);
  for (int argh=0;argh<8&&w[argh];argh++)
    val[argh] = get_eps(ilocs[argh]);
  double epshere = 0.0, epsmin = 1e100, epsmax = -1e100;
  for (int argh=0;argh<8&&w[argh];argh++) {
    epsmin = min(epsmin, val[argh]);
    epsmax = max(epsmax, val[argh]);
    epshere += w[argh] * val[argh];
  }
  vec grad = zero_vec(v.dim);
  vec center = zero_vec(v.dim);
  double norm = 0.0, mean = 0.0;
  for (int argh=0;argh<8&&w[argh];argh++) {
    center += v[ilocs[argh]];
    mean += val[argh];
    norm++;
  }
  center = center/norm;
  mean *= 1.0/norm;
  for (int argh=0;argh<8&&w[argh];argh++) {
    const vec dist = v[ilocs[argh]] - center;
    grad += dist*(val[argh] - mean)/(dist & dist);
  }
  for (int argh=0;argh<8&&w[argh];argh++)
    if (val[argh] != val[0]){
      *loc = *loc + grad*(0.5*(epsmax - epshere) - mean)/abs(grad & grad);
      return true;
    }
  return false;
}

void fields::eps_slices(const geometric_volume &what, const char *name) {
  am_now_working_on(FieldOutput);
  const int buflen = 1024;
  char nname[buflen];
  if (*name) snprintf(nname, buflen, "%s-", name);
  else *nname = 0; // No additional name!
  char *n = (char *)malloc(buflen);
  if (!n) abort("Allocation failure!\n");
  char time_step_string[buflen];
  snprintf(time_step_string, buflen, "%09.2f", time());
  FOR_COMPONENTS(c)
    if (v.has_field(c)) {
      snprintf(n, buflen, "%s/%s%s-%s.eps", outdir, nname,
               component_name(c), time_step_string);
      const double fmax = maxfieldmag_to_master(c);
      file *out = everyone_open_write(n);
      if (!out) {
        printf("Unable to open file '%s' for slice output.\n", n);
        return;
      }
      if (am_master())
        output_complex_eps_header(c, fmax, user_volume, what,
                                  out, n, v.eps_component());
      complex<double> phshift = optimal_phase_shift(c);
      all_wait();
      for (int i=0;i<num_chunks;i++)
        if (chunks[i]->is_mine())
          for (int sn=0;sn<S.multiplicity();sn++)
            FOR_COMPONENTS(otherc)
              if (S.transform(otherc,sn) == c)
                chunks[i]->output_eps_body(otherc,
                                           S, sn, what, out, phshift);
      all_wait();
      for (int i=0;i<num_chunks;i++)
        if (chunks[i]->is_mine())
          for (int sn=0;sn<S.multiplicity();sn++)
            FOR_COMPONENTS(otherc)
              if (S.transform(otherc,sn) == c)
                eps_outline(v.eps_component(), chunks[i]->s->eps,
                            chunks[i]->v, what, S, sn, out);
      all_wait();
      outline_chunks(out);
      all_wait();
      if (am_master()) output_complex_eps_tail(out);
      everyone_close(out);
    }
  free(n);
  finished_working();
}

double fields::minpolenergy_to_master() const {
  double themin = 1e300;
  for (int i=0;i<num_chunks;i++)
    if (chunks[i]->is_mine())
      themin = min(themin,chunks[i]->minpolenergy());
  return -max_to_master(-themin);
}

double fields::maxpolenergy_to_master() const {
  double themax = -1e300;
  for (int i=0;i<num_chunks;i++)
    if (chunks[i]->is_mine())
      themax = max(themax,chunks[i]->maxpolenergy());
  return max_to_master(themax);
}

double fields::maxfieldmag_to_master(component c) const {
  double themax = 0.0;
  for (int i=0;i<num_chunks;i++)
    if (chunks[i]->is_mine())
      themax = max(themax,chunks[i]->maxfieldmag(c));
  return max_to_master(themax);
}

double fields_chunk::maxfieldmag(component c) const {
  double themax = 0.0;
  if (f[c][0])
    for (int i=0;i<v.ntot();i++) {
      double norm = 0.0;
      DOCMP norm += f[c][cmp][i]*f[c][cmp][i];
      themax = max(themax,sqrt(norm));
    }
  return themax;
}

double fields_chunk::minpolenergy() const {
  double themin = 1e300;
  const component c = v.eps_component();
  for (int i=0;i<v.ntot();i++)
    themin = min(themin,my_polarization_energy(v.iloc(c, i)));
  return themin;
}

double fields_chunk::maxpolenergy() const {
  double themax = -1e300;
  const component c = v.eps_component();
  for (int i=0;i<v.ntot();i++)
    themax = max(themax,my_polarization_energy(v.iloc(c, i)));
  return themax;
}

} // namespace meep