/ check that the coordinates and weights computed from the */ metadata match the correct values for all grid points in */ where. return true for zero mismatches, false otherwise. */ / if the environment variable MEEP_ARRAY_METADATA_LOGFILE is */ set, more detailed output is written to that file. */ /
| 69 | /* set, more detailed output is written to that file. */ |
| 70 | /***************************************************************/ |
| 71 | bool test_array_metadata(meep::fields &f, const volume &where) { |
| 72 | /***************************************************************/ |
| 73 | /* step 1: get coordinate grids and weights as reported by */ |
| 74 | /* get_array_metadata */ |
| 75 | /***************************************************************/ |
| 76 | size_t dims[3]; |
| 77 | direction dirs[3]; |
| 78 | int rank = f.get_array_slice_dimensions(where, dims, dirs); |
| 79 | std::vector<double> xyzw = f.get_array_metadata(where); |
| 80 | |
| 81 | // convert to a more convenient format |
| 82 | int offset = 0; |
| 83 | size_t nxyz[3], nw = 1; |
| 84 | vector<double> tics[3], weights; |
| 85 | for (int i = 0; i < 3; ++i) { |
| 86 | nxyz[i] = (size_t)xyzw[offset++]; |
| 87 | nw *= nxyz[i]; |
| 88 | for (size_t j = 0; j < nxyz[i]; ++j) |
| 89 | tics[i].push_back(xyzw[offset++]); |
| 90 | } |
| 91 | for (size_t j = 0; j < nw; ++j) |
| 92 | weights.push_back(xyzw[offset++]); |
| 93 | |
| 94 | size_t stride[3]; |
| 95 | stride[2] = 1; |
| 96 | stride[1] = nxyz[2]; |
| 97 | stride[0] = nxyz[1] * nxyz[2]; |
| 98 | |
| 99 | printf("Metadata: Rank=%i, dims=", rank); |
| 100 | for (int r = 0; r < rank; r++) |
| 101 | printf("%c %zu", r == 0 ? '{' : ',', dims[r]); |
| 102 | printf("}, "); |
| 103 | printf("xyz sizes={%zu, %zu, %zu}, ", nxyz[0], nxyz[1], nxyz[2]); |
| 104 | printf("strides={%zu, %zu, %zu}\n", stride[0], stride[1], stride[2]); |
| 105 | |
| 106 | /***************************************************************/ |
| 107 | /* step 2: initialize loop over grid points in the volume via */ |
| 108 | /* standard libmeep looping primitives */ |
| 109 | /***************************************************************/ |
| 110 | component cgrid = Centered; |
| 111 | grid_volume gv = f.gv; |
| 112 | vec yee_c(gv.yee_shift(Centered) - gv.yee_shift(cgrid)); |
| 113 | ivec iyee_c(gv.iyee_shift(Centered) - gv.iyee_shift(cgrid)); |
| 114 | volume wherec(where + yee_c); |
| 115 | ivec is(vec2diel_floor(wherec.get_min_corner(), gv.a, zero_ivec(gv.dim))); |
| 116 | ivec ie(vec2diel_ceil(wherec.get_max_corner(), gv.a, zero_ivec(gv.dim))); |
| 117 | |
| 118 | ivec imin = gv.little_corner() + one_ivec(gv.dim), imax = gv.big_corner() - one_ivec(gv.dim); |
| 119 | LOOP_OVER_DIRECTIONS(gv.dim, d) { |
| 120 | if (is.in_direction(d) < imin.in_direction(d)) is.set_direction(d, imin.in_direction(d)); |
| 121 | if (ie.in_direction(d) > imax.in_direction(d)) ie.set_direction(d, imax.in_direction(d)); |
| 122 | } |
| 123 | |
| 124 | bool snap_empty_dims = true; |
| 125 | vec s0(gv.dim), e0(gv.dim), s1(gv.dim), e1(gv.dim); |
| 126 | // this initialization step seems to be necessary here to avoid winding |
| 127 | // up with zero or undefined integration weights; I don't know why it |
| 128 | // seems to be unnecessary for loop_in_chunks above. |
no test coverage detected