h5fields.cpp

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

    if (rank >= 3) abort("too many dimensions in output_hdf5");
    int n = (data.max_corner.in_direction(d)
	     - data.min_corner.in_direction(d)) / 2 + 1;
    if (n > 1) {
      data.ds[rank] = d;
      dims[rank++] = n;
    }
  }
  data.rank = rank;

  file->create_or_extend_data(dataname, rank, dims,
                              append_data, single_precision);

  data.buf = new double[data.bufsz];

  data.num_fields = num_fields;
  data.components = components;
  data.cS = new component[num_fields];
  data.ph = new complex<double>[num_fields];
  data.fields = new complex<double>[num_fields];
  data.sum = 0;
  data.fun = fun;
  data.fun_data_ = fun_data_;

  /* compute inverse-epsilon directions for computing Dielectric fields */
  data.ninveps = 0;
  bool needs_dielectric = false;
  for (int i = 0; i < num_fields; ++i)
    if (components[i] == Dielectric) { needs_dielectric = true; break; }
  if (needs_dielectric) 
    FOR_ELECTRIC_COMPONENTS(c) if (v.has_field(c)) {
      if (data.ninveps == 3) abort("more than 3 field components??");
      data.inveps_cs[data.ninveps] = c;
      data.inveps_ds[data.ninveps] = component_direction(c);
      ++data.ninveps;
    }
  
  /* compute inverse-mu directions for computing Permeability fields */
  data.ninvmu = 0;
  bool needs_permeability = false;
  for (int i = 0; i < num_fields; ++i)
    if (components[i] == Permeability) { needs_permeability = true; break; }
  if (needs_permeability) 
    FOR_MAGNETIC_COMPONENTS(c) if (v.has_field(c)) {
      if (data.ninvmu == 3) abort("more than 3 field components??");
      data.invmu_cs[data.ninvmu] = c;
      data.invmu_ds[data.ninvmu] = component_direction(c);
      ++data.ninvmu;
    }
  
  data.offsets = new int[2 * num_fields];
  for (int i = 0; i < 2 * num_fields; ++i)
    data.offsets[i] = 0;
  
  loop_in_chunks(h5_output_chunkloop, (void *) &data, 
		 where, Dielectric, true, true);

  delete[] data.offsets;
  delete[] data.fields;
  delete[] data.ph;
  delete[] data.cS;
  delete[] data.buf;
  file->done_writing_chunks();
  finished_working();
}

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

void fields::output_hdf5(const char *dataname,
                         int num_fields, const component *components,
                         field_function fun, void *fun_data_,
                         const geometric_volume &where,
			 h5file *file,
                         bool append_data,
                         bool single_precision,
			 const char *prefix,
			 bool real_part_only)
{
  bool delete_file;
  if ((delete_file = !file))
    file = open_h5file(dataname, h5file::WRITE, prefix, true);

  if (real_part_only) {
    output_hdf5(file, dataname, num_fields, components, fun, fun_data_, 
		0, where, append_data, single_precision);
  }
  else {
    int len = strlen(dataname) + 5;
    char *dataname2 = new char[len];
    snprintf(dataname2, len, "%s%s", dataname, ".r");
    output_hdf5(file, dataname2, num_fields, components, fun, fun_data_, 
		0, where, append_data, single_precision);
    snprintf(dataname2, len, "%s%s", dataname, ".i");
    output_hdf5(file, dataname2, num_fields, components, fun, fun_data_, 
		1, where, append_data, single_precision);
    delete[] dataname2;
  }
  if (delete_file) delete file;
}

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

typedef struct {
  field_rfunction fun;
  void *fun_data_;
} rintegrand_data;

static complex<double> rintegrand_fun(const complex<double> *fields,
                                     const vec &loc,
                                     void *data_)
{
  rintegrand_data *data = (rintegrand_data *) data_;
  return data->fun(fields, loc, data->fun_data_);
}

void fields::output_hdf5(const char *dataname,
                         int num_fields, const component *components,
                         field_rfunction fun, void *fun_data_,
                         const geometric_volume &where,
			 h5file *file,
                         bool append_data,
                         bool single_precision,
			 const char *prefix)
{
  bool delete_file;
  if ((delete_file = !file))
    file = open_h5file(dataname, h5file::WRITE, prefix, true);

  rintegrand_data data; data.fun = fun; data.fun_data_ = fun_data_;
  output_hdf5(file, dataname, num_fields, components, rintegrand_fun,
	      (void *) &data, 0, where, append_data, single_precision);

  if (delete_file) delete file;
}

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

static complex<double> component_fun(const complex<double> *fields,
				     const vec &loc,
				     void *data_)
{
     (void) loc; // unused
     (void) data_; // unused
     return fields[0];
}

void fields::output_hdf5(component c,
			 const geometric_volume &where,
			 h5file *file,
			 bool append_data,
                         bool single_precision,
			 const char *prefix) {
  if (is_derived(int(c))) {
    output_hdf5(derived_component(c), 
		where, file, append_data, single_precision, prefix);
    return;
  }

  if (coordinate_mismatch(v.dim, c)) return;

  char dataname[256];
  bool has_imag = !is_real && c != Dielectric && c != Permeability;

  bool delete_file;
  if ((delete_file = !file))
    file = open_h5file(component_name(c), h5file::WRITE, prefix, true);

  snprintf(dataname, 256, "%s%s", component_name(c), has_imag ? ".r" : "");
  output_hdf5(file, dataname, 1, &c, component_fun, 0, 0, where,
	      append_data, single_precision);
  if (has_imag) {
    snprintf(dataname, 256, "%s.i", component_name(c));
    output_hdf5(file, dataname, 1, &c, component_fun, 0, 1, where,
		append_data, single_precision);
  }

  if (delete_file) delete file;
}

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

void fields::output_hdf5(derived_component c,
			 const geometric_volume &where,
			 h5file *file,
			 bool append_data,
                         bool single_precision,
			 const char *prefix) {
  if (!is_derived(int(c))) {
    output_hdf5(component(c), 
		where, file, append_data, single_precision, prefix);
    return;
  }

  if (coordinate_mismatch(v.dim, c)) return;

  int nfields;
  component cs[12];
  field_rfunction fun = derived_component_func(c, v, nfields, cs);

  output_hdf5(component_name(c), nfields, cs, fun, &nfields, where,
	      file, append_data, single_precision, prefix);
}

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

const char *fields::h5file_name(const char *name, 
				const char *prefix, bool timestamp)
{
  const int buflen = 1024;
  static char filename[buflen];
  char time_step_string[32] = "";

  if (timestamp) {
    if (dt >= 0.01 && dt < 10)
      snprintf(time_step_string, 32, "-%09.2f", time());
    else
      snprintf(time_step_string, 32, "-%09d", t);
  } 

  snprintf(filename, buflen, "%s/" "%s%s" "%s" "%s" ".h5",
	   outdir,
	   prefix ? prefix : "", prefix && prefix[0] ? "-" : "",
	   name, time_step_string);
  return filename;
}

h5file *fields::open_h5file(const char *name, h5file::access_mode mode, 
			    const char *prefix, bool timestamp)
{
  const char *filename = h5file_name(name, prefix, timestamp);
  if (!quiet && mode == h5file::WRITE)
    master_printf("creating output file \"%s\"...\n", filename);
  return new h5file(filename, mode, true);
}

} // namespace meep