symmetry.cpp

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

  }
  return 1;
}

int test_3D_metal_odd_zmirror(double eps(const vec &)) {
  double a = 8.0;
  double ttot = 3.0;

  const volume v = vol3d(1.1, 0.6, 1.0, a);
  const symmetry S = -mirror(Z,v);
  structure s(v, eps, no_pml(), S);
  structure s1(v, eps);
  master_printf("Testing odd Z mirror symmetry in 3D...\n");

  fields f1(&s1);
  f1.add_point_source(Ez, 0.7, 2.5, 0.0, 4.0, vec(0.55,0.51,0.5));
  fields f(&s);
  f.add_point_source(Ez, 0.7, 2.5, 0.0, 4.0, vec(0.55,0.51,0.5));
  check_unequal_layout(f, f1);
  double total_energy_check_time = 1.0;
  while (f.time() < ttot) {
    f.step();
    f1.step();
    if (!compare_point(f, f1, vec(0.5  , 0.01 , 0.75))) return 0;
    if (!compare_point(f, f1, vec(0.5  , 0.21 , 0.15))) return 0;
    if (!compare_point(f, f1, vec(0.5  , 0.501, 0.5))) return 0;
    if (!compare_point(f, f1, vec(0.46 , 0.33 , 0.51))) return 0;
    if (!compare_point(f, f1, vec(0.2  , 0.2  , 0.05))) return 0;
    if (f.time() >= total_energy_check_time) {
      if (!compare(f.electric_energy_in_box(v.surroundings()),
                   f1.electric_energy_in_box(v.surroundings()),
                   "electric energy")) return 0;
      if (!compare(f.magnetic_energy_in_box(v.surroundings()),
                   f1.magnetic_energy_in_box(v.surroundings()),
                   "magnetic energy")) return 0;
      if (!compare(f.total_energy(), f1.total_energy(),
                   "   total energy")) return 0;
      total_energy_check_time += 1.0;
    }
  }
  return 1;
}

int test_3D_metal_rot4z(double eps(const vec &)) {
  double a = 8.0;
  double ttot = 3.0;

  const volume v = vol3d(1.0, 1.0, 1.0, a);
  const symmetry S = rotate4(Z,v);
  structure s(v, eps, no_pml(), S);
  structure s1(v, eps);
  master_printf("Testing Z fourfold rotational symmetry in 3D...\n");

  fields f1(&s1);
  f1.add_point_source(Ez, 0.7, 2.5, 0.0, 4.0, vec(0.5,0.5,0.52));
  f1.add_point_source(Hz, 0.8, 0.6, 0.0, 4.0, vec(0.5,0.5,0.43));
  fields f(&s);
  f.add_point_source(Ez, 0.7, 2.5, 0.0, 4.0, vec(0.5,0.5,0.52));
  f.add_point_source(Hz, 0.8, 0.6, 0.0, 4.0, vec(0.5,0.5,0.43));
  check_unequal_layout(f, f1);
  double total_energy_check_time = 1.0;
  while (f.time() < ttot) {
    f.step();
    f1.step();
    if (!compare_point(f, f1, vec(0.5  , 0.01 , 0.75))) return 0;
    if (!compare_point(f, f1, vec(0.5  , 0.21 , 0.15))) return 0;
    if (!compare_point(f, f1, vec(0.5  , 0.501, 0.5))) return 0;
    if (!compare_point(f, f1, vec(0.46 , 0.33 , 0.51))) return 0;
    if (!compare_point(f, f1, vec(0.2  , 0.2  , 0.05))) return 0;
    if (f.time() >= total_energy_check_time) {
      if (!compare(f.electric_energy_in_box(v.surroundings()),
                   f1.electric_energy_in_box(v.surroundings()),
                   "electric energy")) return 0;
      if (!compare(f.magnetic_energy_in_box(v.surroundings()),
                   f1.magnetic_energy_in_box(v.surroundings()),
                   "magnetic energy")) return 0;
      if (!compare(f.total_energy(), f1.total_energy(),
                   "   total energy")) return 0;
      total_energy_check_time += 1.0;
    }
  }
  return 1;
}

int test_3D_metal_rot4z_mirror(double eps(const vec &)) {
  double a = 8.0;
  double ttot = 3.0;

  const volume v = vol3d(1.0, 1.0, 1.0, a);
  const symmetry S = rotate4(Z,v) + mirror(Z,v);
  structure s(v, eps, no_pml(), S);
  structure s1(v, eps);
  master_printf("Testing Z fourfold rotational symmetry in 3D with horizontal mirror...\n");

  fields f1(&s1);
  f1.add_point_source(Hz, 0.7, 2.5, 0.0, 4.0, vec(0.5,0.5,0.5));
  fields f(&s);
  f.add_point_source(Hz, 0.7, 2.5, 0.0, 4.0, vec(0.5,0.5,0.5));
  check_unequal_layout(f, f1);
  double total_energy_check_time = 1.0;
  while (f.time() < ttot) {
    f.step();
    f1.step();
    if (!compare_point(f, f1, vec(0.5  , 0.01 , 0.75))) return 0;
    if (!compare_point(f, f1, vec(0.5  , 0.21 , 0.15))) return 0;
    if (!compare_point(f, f1, vec(0.5  , 0.501, 0.5))) return 0;
    if (!compare_point(f, f1, vec(0.46 , 0.33 , 0.51))) return 0;
    if (!compare_point(f, f1, vec(0.2  , 0.2  , 0.05))) return 0;
    if (f.time() >= total_energy_check_time) {
      if (!compare(f.electric_energy_in_box(v.surroundings()),
                   f1.electric_energy_in_box(v.surroundings()),
                   "electric energy")) return 0;
      if (!compare(f.magnetic_energy_in_box(v.surroundings()),
                   f1.magnetic_energy_in_box(v.surroundings()),
                   "magnetic energy")) return 0;
      if (!compare(f.total_energy(), f1.total_energy(),
                   "   total energy")) return 0;
      total_energy_check_time += 1.0;
    }
  }
  return 1;
}

int test_3D_metal_3mirror(double eps(const vec &)) {
  double a = 8.0;
  double ttot = 3.0;

  const volume v = vol3d(1.0, 1.0, 1.0, a);
  const symmetry S = mirror(Z,v) - mirror(Y,v) - mirror(X,v);
  structure s(v, eps, no_pml(), S);
  structure s1(v, eps);
  master_printf("Testing three mirror planes in 3D...\n");

  fields f1(&s1);
  f1.add_point_source(Hz, 0.7, 2.5, 0.0, 4.0, vec(0.5,0.5,0.5));
  fields f(&s);
  f.add_point_source(Hz, 0.7, 2.5, 0.0, 4.0, vec(0.5,0.5,0.5));
  check_unequal_layout(f, f1);
  double total_energy_check_time = 1.0;
  while (f.time() < ttot) {
    f.step();
    f1.step();
    if (!compare_point(f, f1, vec(0.5  , 0.01 , 0.75))) return 0;
    if (!compare_point(f, f1, vec(0.5  , 0.21 , 0.15))) return 0;
    if (!compare_point(f, f1, vec(0.5  , 0.501, 0.5))) return 0;
    if (!compare_point(f, f1, vec(0.46 , 0.33 , 0.51))) return 0;
    if (!compare_point(f, f1, vec(0.2  , 0.2  , 0.05))) return 0;
    if (f.time() >= total_energy_check_time) {
      if (!compare(f.electric_energy_in_box(v.surroundings()),
                   f1.electric_energy_in_box(v.surroundings()),
                   "electric energy")) return 0;
      if (!compare(f.magnetic_energy_in_box(v.surroundings()),
                   f1.magnetic_energy_in_box(v.surroundings()),
                   "magnetic energy")) return 0;
      if (!compare(f.total_energy(), f1.total_energy(),
                   "   total energy")) return 0;
      total_energy_check_time += 1.0;
    }
  }
  return 1;
}

int test_metal_ymirror(double eps(const vec &)) {
  double a = 8.0;
  double ttot = 5.0;

  const volume v = voltwo(1.0, 1.0, a);
  the_center = v.center();
  const symmetry S = mirror(Y,v);
  structure s(v, eps, no_pml(), S);
  structure s1(v, eps);
  master_printf("Testing Y mirror symmetry...\n");

  fields f1(&s1);
  f1.add_point_source(Ex, 0.7, 2.5, 0.0, 4.0, vec(0.85 ,0.5));
  f1.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(0.401,0.5));
  fields f(&s);
  f.add_point_source(Ex, 0.7, 2.5, 0.0, 4.0, vec(0.85 ,0.5));
  f.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(0.401,0.5));
  check_unequal_layout(f, f1);
  double total_energy_check_time = 1.0;
  while (f.time() < ttot) {
    f.step();
    f1.step();
    if (!compare_point(f, f1, vec(0.01 ,  0.5))) return 0;
    if (!compare_point(f, f1, vec(0.21 ,  0.5))) return 0;
    if (!compare_point(f, f1, vec(0.46 , 0.33))) return 0;
    if (!compare_point(f, f1, vec(0.2  , 0.2 ))) return 0;
    if (f.time() >= total_energy_check_time) {
      if (!compare(f.electric_energy_in_box(v.surroundings()),
                   f1.electric_energy_in_box(v.surroundings()),
                   "electric energy")) return 0;
      if (!compare(f.magnetic_energy_in_box(v.surroundings()),
                   f1.magnetic_energy_in_box(v.surroundings()),
                   "magnetic energy")) return 0;
      if (!compare(f.total_energy(), f1.total_energy(),
                   "   total energy")) return 0;
      total_energy_check_time += 1.0;
    }
  }
  return 1;
}

int test_yperiodic_ymirror(double eps(const vec &)) {
  double a = 8.0;
  double ttot = 5.0;

  const volume v = voltwo(1.0, 1.0, a);
  the_center = v.center();
  const symmetry S = mirror(Y,v);
  structure s(v, eps, no_pml(), S);
  structure s1(v, eps);
  s.set_output_directory(mydirname);
  s1.set_output_directory(mydirname);
  master_printf("Testing Y periodic with mirror symmetry...\n");

  fields f1(&s1);
  f1.use_bloch(vec(0.1*pi/2,0.0));
  //f1.add_point_source(Ex, 0.7, 2.5, 0.0, 4.0, vec(0.85 ,0.5));
  f1.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(0.401,0.5));
  fields f(&s);
  f.use_bloch(vec(0.1*pi/2,0.0));
  //f.add_point_source(Ex, 0.7, 2.5, 0.0, 4.0, vec(0.85 ,0.5));
  f.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(0.401,0.5));
  check_unequal_layout(f, f1);
  double total_energy_check_time = 1.0;
  while (f.time() < ttot) {
    f.step();
    f1.step();
    if (!compare_point(f, f1, vec(0.951 ,  0.5))) return 0;
    if (!compare_point(f, f1, vec(0.01 ,  0.5))) return 0;
    if (!compare_point(f, f1, vec(0.21 ,  0.5))) return 0;
    if (!compare_point(f, f1, vec(0.46 , 0.33))) return 0;
    if (!compare_point(f, f1, vec(0.2  , 0.2 ))) return 0;
    if (f.time() >= total_energy_check_time) {
      if (!compare(f.electric_energy_in_box(v.surroundings()),
                   f1.electric_energy_in_box(v.surroundings()),
                   "electric energy")) {
        f.output_real_imaginary_slices("multi");
        f1.output_real_imaginary_slices("single");
        return 0;
      }
      if (!compare(f.magnetic_energy_in_box(v.surroundings()),
                   f1.magnetic_energy_in_box(v.surroundings()),
                   "magnetic energy")) return 0;
      if (!compare(f.total_energy(), f1.total_energy(),
                   "   total energy")) return 0;
      total_energy_check_time += 1.0;
    }
  }
  return 1;
}

int test_metal_rot2y(double eps(const vec &)) {
  double a = 16.0;
  double ttot = 5.0;

  const volume v = voltwo(1.0, 1.0, a);
  the_center = v.center();
  const symmetry S = rotate2(Y,v);
  structure s(v, eps, no_pml(), S);
  structure s1(v, eps);
  master_printf("Testing Y twofold rotational symmetry...\n");

  fields f1(&s1);
  f1.add_point_source(Hz, 0.7, 2.5, 0.0, 4.0, vec(0.25, 0.875), 1.0);
  f1.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(0.25,0.375), 1.0);
  f1.add_point_source(Hz, 0.7, 2.5, 0.0, 4.0, vec(0.75, 0.875),-1.0);
  f1.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(0.75,0.375),-1.0);
  fields f(&s);
  f.add_point_source(Hz, 0.7, 2.5, 0.0, 4.0, vec(0.25,0.875 ), 1.0);
  f.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(0.25,0.375), 1.0);
  f.add_point_source(Hz, 0.7, 2.5, 0.0, 4.0, vec(0.75,0.875 ),-1.0);
  f.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(0.75,0.375),-1.0);
  check_unequal_layout(f, f1);
  double total_energy_check_time = 1.0;
  while (f.time() < ttot) {
    f.step();
    f1.step();
    if (!compare_point(f, f1, vec(0.01 ,  0.5))) return 0;
    if (!compare_point(f, f1, vec(0.21 ,  0.5))) return 0;
    if (!compare_point(f, f1, vec(0.46 , 0.33))) return 0;
    if (!compare_point(f, f1, vec(0.2  , 0.2 ))) return 0;
    if (f.time() >= total_energy_check_time) {
      if (!compare(f.electric_energy_in_box(v.surroundings()),
                   f1.electric_energy_in_box(v.surroundings()),
                   "electric energy")) return 0;
      if (!compare(f.magnetic_energy_in_box(v.surroundings()),
                   f1.magnetic_energy_in_box(v.surroundings()),
                   "magnetic energy")) return 0;
      if (!compare(f.total_energy(), f1.total_energy(),
                   "   total energy")) return 0;
      total_energy_check_time += 1.0;
    }
  }
  return 1;
}

int exact_metal_rot2y(double eps(const vec &)) {
  double a = 16.0;
  double ttot = 5.0;

  const volume v = voltwo(1.0, 1.5, a);
  the_center = v.center();
  const symmetry S = rotate2(Y,v);
  structure s(v, eps, no_pml(), S);
  structure s1(v, eps);
  master_printf("Testing exact Y twofold rotational symmetry...\n");

  fields f1(&s1);
  f1.add_point_source(Ey, 0.7, 2.5, 0.0, 4.0, vec(0.5, 0.875));
  f1.add_point_source(Hy, 0.8, 0.6, 0.0, 4.0, vec(0.5,0.375));
  fields f(&s);
  f.add_point_source(Ey, 0.7, 2.5, 0.0, 4.0, vec(0.5, 0.875));
  f.add_point_source(Hy, 0.8, 0.6, 0.0, 4.0, vec(0.5,0.375));
  check_unequal_layout(f, f1);
  double total_energy_check_time = 1.0;
  while (f.time() < ttot) {
    f.step();
    f1.step();
    if (!compare_point(f, f1, vec(0.01 ,  0.5))) return 0;
    if (!compare_point(f, f1, vec(0.21 ,  0.5))) return 0;
    if (!compare_point(f, f1, vec(0.46 , 0.33))) return 0;
    if (!compare_point(f, f1, vec(0.2  , 0.2 ))) return 0;
    if (f.time() >= total_energy_check_time) {
      if (!compare(f.electric_energy_in_box(v.surroundings()),
                   f1.electric_energy_in_box(v.surroundings()),
                   "electric energy")) return 0;
      if (!compare(f.magnetic_energy_in_box(v.surroundings()),
                   f1.magnetic_energy_in_box(v.surroundings()),
                   "magnetic energy")) return 0;
      if (!compare(f.total_energy(), f1.total_energy(),
                   "   total energy")) return 0;
      total_energy_check_time += 1.0;
    }
  }
  return 1;
}

int pml_twomirrors(double eps(const vec &)) {
  double a = 16.0;
  double ttot = 10.0;

  const volume v = voltwo(2.0, 2.0, a);
  the_center = v.center();
  const symmetry S = mirror(X,v) + mirror(Y,v);

  structure s_mm(v, eps, pml(0.5), S);
  structure s1(v, eps, pml(0.5), identity());
  structure ss[2] = { s1, s_mm };

  master_printf("Testing two mirrors with PML...\n");

  fields fs[2] = { fields(&ss[0]), fields(&ss[1]) };
  for (int i=0;i<2;i++) {
    fs[i].add_point_source(Ez, 0.7, 2.5, 0.0, 4.0, vec(1.0,1.0),-1.5);
    fs[i].add_point_source(Ez, 0.7, 2.5, 0.0, 4.0, vec(0.75,0.75));
    fs[i].add_point_source(Ez, 0.7, 2.5, 0.0, 4.0, vec(0.75,1.25));
    fs[i].add_point_source(Ez, 0.7, 2.5, 0.0, 4.0, vec(1.25,0.75));
    fs[i].add_point_source(Ez, 0.7, 2.5, 0.0, 4.0, vec(1.25,1.25));
  }
  check_unequal_layout(fs[0], fs[1]);
  double total_energy_check_time = 3.0;
  while (fs[0].time() < ttot) {
    for (int i=0;i<2;i++) fs[i].step();
    if (!compare_point(fs[1], fs[0], vec(0.01 ,  0.5))) return 0;
    if (!compare_point(fs[1], fs[0], vec(0.21 ,  0.5))) return 0;
    if (!compare_point(fs[1], fs[0], vec(0.46 , 0.33))) return 0;
    if (!compare_point(fs[1], fs[0], vec(0.2  , 0.2 ))) return 0;
    if (fs[0].time() >= total_energy_check_time) {
      if (!compare(fs[0].electric_energy_in_box(v.surroundings()),
                   fs[1].electric_energy_in_box(v.surroundings()),