diff --git a/include/openmc/plot.h b/include/openmc/plot.h index be9030fa7..689a156a7 100644 --- a/include/openmc/plot.h +++ b/include/openmc/plot.h @@ -78,6 +78,8 @@ void output_ppm(ObjectPlot* pl, extern "C" void create_ppm(ObjectPlot* pl); +extern "C" void create_voxel(ObjectPlot *pl); + extern "C" void position_rgb(Particle* p, ObjectPlot* pl, int rgb[3], int &id); extern "C" void voxel_init(hid_t file_id, const hsize_t* dims, hid_t* dspace, diff --git a/src/plot.F90 b/src/plot.F90 index 8e41dcb1b..79b286924 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -40,6 +40,11 @@ module plot type(ObjectPlot), intent(in) :: pl end subroutine create_ppm + subroutine create_voxel(pl) bind(C) + import ObjectPlot + type(ObjectPlot), intent(in) :: pl + end subroutine create_voxel + end interface contains @@ -71,130 +76,4 @@ contains err = 0 end function openmc_plot_geometry -!=============================================================================== -! CREATE_VOXEL outputs a binary file that can be input into silomesh for 3D -! geometry visualization. It works the same way as create_ppm by dragging a -! particle across the geometry for the specified number of voxels. The first -! 3 int(4)'s in the binary are the number of x, y, and z voxels. The next 3 -! real(8)'s are the widths of the voxels in the x, y, and z directions. The next -! 3 real(8)'s are the x, y, and z coordinates of the lower left point. Finally -! the binary is filled with entries of four int(4)'s each. Each 'row' in the -! binary contains four int(4)'s: 3 for x,y,z position and 1 for cell or material -! id. For 1 million voxels this produces a file of approximately 15MB. -!=============================================================================== - - subroutine create_voxel(pl) - type(ObjectPlot), intent(in) :: pl - - integer(C_INT) :: x, y, z ! voxel location indices - integer :: rgb(3) ! colors (red, green, blue) from 0-255 - integer :: id ! id of cell or material - integer(C_INT), target :: data(pl%pixels(3),pl%pixels(2)) - integer(HID_T) :: file_id - integer(HID_T) :: dspace - integer(HID_T) :: memspace - integer(HID_T) :: dset - integer(HSIZE_T) :: dims(3) - real(8) :: vox(3) ! x, y, and z voxel widths - real(8) :: ll(3) ! lower left starting point for each sweep direction - type(Particle) :: p - type(ProgressBar) :: progress - - interface - subroutine voxel_init(file_id, dims, dspace, dset, memspace) bind(C) - import HID_T, HSIZE_T - integer(HID_T), value :: file_id - integer(HSIZE_T), intent(in) :: dims(*) - integer(HID_T), intent(out) :: dspace - integer(HID_T), intent(out) :: dset - integer(HID_T), intent(out) :: memspace - end subroutine voxel_init - - subroutine voxel_write_slice(x, dspace, dset, memspace, buf) bind(C) - import C_INT, HID_T, C_PTR - integer(C_INT), value :: x - integer(HID_T), value :: dspace - integer(HID_T), value :: dset - integer(HID_T), value :: memspace - type(C_PTR), value :: buf - end subroutine voxel_write_slice - - subroutine voxel_finalize(dspace, dset, memspace) bind(C) - import HID_T - integer(HID_T), value :: dspace - integer(HID_T), value :: dset - integer(HID_T), value :: memspace - end subroutine voxel_finalize - end interface - - ! compute voxel widths in each direction - vox = pl % width/dble(pl % pixels) - - ! initial particle position - ll = pl % origin - pl % width / TWO - - ! allocate and initialize particle - call particle_initialize(p) - p % coord(1) % xyz = ll - p % coord(1) % uvw = [ HALF, HALF, HALF ] - p % coord(1) % universe = root_universe - - ! Open binary plot file for writing - file_id = file_open(pl%path_plot, 'w') - - ! write header info - call write_attribute(file_id, "filetype", 'voxel') - call write_attribute(file_id, "version", VERSION_VOXEL) - call write_attribute(file_id, "openmc_version", VERSION) -#ifdef GIT_SHA1 - call write_attribute(file_id, "git_sha1", GIT_SHA1) -#endif - - ! Write current date and time - call write_attribute(file_id, "date_and_time", time_stamp()) - - call write_attribute(file_id, "num_voxels", pl%pixels) - call write_attribute(file_id, "voxel_width", vox) - call write_attribute(file_id, "lower_left", ll) - - ! Create dataset for voxel data -- note that the dimensions are reversed - ! since we want the order in the file to be z, y, x - dims(:) = pl % pixels - call voxel_init(file_id, dims, dspace, dset, memspace) - - ! move to center of voxels - ll = ll + vox / TWO - - do x = 1, pl % pixels(1) - call progress % set_value(dble(x)/dble(pl % pixels(1))*100) - do y = 1, pl % pixels(2) - do z = 1, pl % pixels(3) - ! get voxel color - call position_rgb(p, pl, rgb, id) - ! write to plot file - data(z,y) = id - - ! advance particle in z direction - p % coord(1) % xyz(3) = p % coord(1) % xyz(3) + vox(3) - end do - - ! advance particle in y direction - p % coord(1) % xyz(2) = p % coord(1) % xyz(2) + vox(2) - p % coord(1) % xyz(3) = ll(3) - end do - - ! advance particle in y direction - p % coord(1) % xyz(1) = p % coord(1) % xyz(1) + vox(1) - p % coord(1) % xyz(2) = ll(2) - p % coord(1) % xyz(3) = ll(3) - - ! Write to HDF5 dataset - call voxel_write_slice(x, dspace, dset, memspace, c_loc(data)) - end do - - call voxel_finalize(dspace, dset, memspace) - call file_close(file_id) - - end subroutine create_voxel - end module plot diff --git a/src/plot.cpp b/src/plot.cpp index 6f48eec22..1042ccd78 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -376,33 +376,27 @@ void create_voxel(ObjectPlot *pl) { #endif // Write current date and time - // write_attribute(file_id, "date_and_time", time_stamp()); + // write_attribute(file_id, "date_and_time", time_stamp()); <-- RE-ADD!!! hsize_t three = 3; write_attr_int(file_id, 1, &three, "num_voxels", pl->pixels); write_attr_double(file_id, 1, &three, "voxel_width", vox); write_attr_double(file_id, 1, &three, "lower_left", ll); + // Create dataset for voxel data -- note that the dimensions are reversed + // since we want the order in the file to be z, y, x hsize_t dims[3]; dims[0] = pl->pixels[0]; dims[1] = pl->pixels[1]; dims[2] = pl->pixels[2]; hid_t dspace, dset, memspace; voxel_init(file_id, &(dims[0]), &dspace, &dset, &memspace); - + + // move to center of voxels ll[0] = ll[0] + vox[0] / TWO; ll[1] = ll[1] + vox[1] / TWO; ll[2] = ll[2] + vox[2] / TWO; - ImageData data; - - - data.resize(pl->pixels[2]); - for (auto & i : data) { - i.resize(pl->pixels[1]); - for (auto & j : i) { - j.resize(1); - } - } + int data[pl->pixels[1]][pl->pixels[2]]; int rgb[3], id; for (int x = 0; x < pl->pixels[0]; x++) { @@ -411,6 +405,8 @@ void create_voxel(ObjectPlot *pl) { for(int z = 0; z < pl->pixels[2]; z++) { // get voxel color position_rgb(p, pl, rgb, id); + // write to plot data + data[y][z] = id; // advance particle in z direction p->coord[0].xyz[2] = p->coord[0].xyz[2] + vox[2]; } @@ -419,12 +415,11 @@ void create_voxel(ObjectPlot *pl) { p->coord[0].xyz[2] = ll[2]; } // advance particle in x direction - p->coord[0].xyz[0] = p->coord[0].xyz[0] + vox[0]; - p->coord[0].xyz[1] = ll[1]; - p->coord[0].xyz[2] = ll[2]; - - // Write to HDF5 dataset - voxel_write_slice(x, dspace, dset, memspace, &(data[0])); + p->coord[0].xyz[0] = p->coord[0].xyz[0] + vox[0]; + p->coord[0].xyz[1] = ll[1]; + p->coord[0].xyz[2] = ll[2]; + // Write to HDF5 dataset + voxel_write_slice(x, dspace, dset, memspace, &(data[0])); } voxel_finalize(dspace, dset, memspace); @@ -456,7 +451,7 @@ voxel_init(hid_t file_id, const hsize_t* dims, hid_t* dspace, hid_t* dset, void voxel_write_slice(int x, hid_t dspace, hid_t dset, hid_t memspace, void* buf) { - hssize_t offset[3] {x - 1, 0, 0}; + hssize_t offset[3] {x, 0, 0}; H5Soffset_simple(dspace, offset); H5Dwrite(dset, H5T_NATIVE_INT, memspace, dspace, H5P_DEFAULT, buf); }