• Home
  • Features
  • Pricing
  • Docs
  • Announcements
  • Sign In

openmc-dev / openmc / 22642461683

03 Mar 2026 08:59PM UTC coverage: 81.553% (+0.01%) from 81.539%
22642461683

Pull #3806

github

web-flow
Merge 148550442 into 53d98ce71
Pull Request #3806: Introduce new C API function for slice plots

17635 of 25352 branches covered (69.56%)

Branch coverage included in aggregate %.

198 of 245 new or added lines in 4 files covered. (80.82%)

1 existing line in 1 file now uncovered.

58018 of 67414 relevant lines covered (86.06%)

44724342.19 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

70.89
/src/plot.cpp
1
#include "openmc/plot.h"
2

3
#include <algorithm>
4
#define _USE_MATH_DEFINES // to make M_PI declared in Intel and MSVC compilers
5
#include <cmath>
6
#include <cstdio>
7
#include <fstream>
8
#include <sstream>
9

10
#include "openmc/tensor.h"
11
#include <fmt/core.h>
12
#include <fmt/ostream.h>
13
#ifdef USE_LIBPNG
14
#include <png.h>
15
#endif
16

17
#include "openmc/cell.h"
18
#include "openmc/constants.h"
19
#include "openmc/container_util.h"
20
#include "openmc/dagmc.h"
21
#include "openmc/error.h"
22
#include "openmc/file_utils.h"
23
#include "openmc/geometry.h"
24
#include "openmc/hdf5_interface.h"
25
#include "openmc/material.h"
26
#include "openmc/mesh.h"
27
#include "openmc/message_passing.h"
28
#include "openmc/openmp_interface.h"
29
#include "openmc/output.h"
30
#include "openmc/particle.h"
31
#include "openmc/progress_bar.h"
32
#include "openmc/random_lcg.h"
33
#include "openmc/settings.h"
34
#include "openmc/simulation.h"
35
#include "openmc/string_utils.h"
36
#include "openmc/tallies/filter.h"
37

38
namespace openmc {
39

40
//==============================================================================
41
// Constants
42
//==============================================================================
43

44
constexpr int PLOT_LEVEL_LOWEST {-1}; //!< lower bound on plot universe level
45
constexpr int32_t NOT_FOUND {-2};
46
constexpr int32_t OVERLAP {-3};
47

48
IdData::IdData(size_t h_res, size_t v_res, bool /*include_filter*/)
5,151✔
49
  : data_({v_res, h_res, 3}, NOT_FOUND)
5,151✔
50
{}
5,151✔
51

52
void IdData::set_value(size_t y, size_t x, const Particle& p, int level,
38,743,747✔
53
  Filter* /*filter*/, FilterMatch* /*match*/)
54
{
55
  // set cell data
56
  if (p.n_coord() <= level) {
38,743,747!
57
    data_(y, x, 0) = NOT_FOUND;
×
58
    data_(y, x, 1) = NOT_FOUND;
×
59
  } else {
60
    data_(y, x, 0) = model::cells.at(p.coord(level).cell())->id_;
38,743,747!
61
    data_(y, x, 1) = level == p.n_coord() - 1
38,743,747✔
62
                       ? p.cell_instance()
38,743,747!
63
                       : cell_instance_at_level(p, level);
×
64
  }
65

66
  // set material data
67
  Cell* c = model::cells.at(p.lowest_coord().cell()).get();
38,743,747✔
68
  if (p.material() == MATERIAL_VOID) {
38,743,747✔
69
    data_(y, x, 2) = MATERIAL_VOID;
29,804,276✔
70
  } else if (c->type_ == Fill::MATERIAL) {
8,939,471!
71
    Material* m = model::materials.at(p.material()).get();
8,939,471✔
72
    data_(y, x, 2) = m->id_;
8,939,471✔
73
  }
74
}
38,743,747✔
75

76
void IdData::set_overlap(size_t y, size_t x)
374,308✔
77
{
78
  for (size_t k = 0; k < data_.shape(2); ++k)
2,994,464!
79
    data_(y, x, k) = OVERLAP;
1,122,924✔
80
}
374,308✔
81

82
PropertyData::PropertyData(size_t h_res, size_t v_res, bool /*include_filter*/)
11✔
83
  : data_({v_res, h_res, 2}, NOT_FOUND)
11✔
84
{}
11✔
85

86
void PropertyData::set_value(size_t y, size_t x, const Particle& p, int level,
99✔
87
  Filter* /*filter*/, FilterMatch* /*match*/)
88
{
89
  Cell* c = model::cells.at(p.lowest_coord().cell()).get();
99!
90
  data_(y, x, 0) = (p.sqrtkT() * p.sqrtkT()) / K_BOLTZMANN;
99✔
91
  if (c->type_ != Fill::UNIVERSE && p.material() != MATERIAL_VOID) {
99!
92
    Material* m = model::materials.at(p.material()).get();
99✔
93
    data_(y, x, 1) = m->density_gpcc_;
99✔
94
  }
95
}
99✔
96

97
void PropertyData::set_overlap(size_t y, size_t x)
×
98
{
99
  data_(y, x) = OVERLAP;
×
100
}
×
101

102
//==============================================================================
103
// RasterData implementation
104
//==============================================================================
105

106
RasterData::RasterData(size_t h_res, size_t v_res, bool include_filter)
110✔
107
  : id_data_({v_res, h_res, include_filter ? 4u : 3u}, NOT_FOUND),
198✔
108
    property_data_({v_res, h_res, 2}, static_cast<double>(NOT_FOUND)),
110✔
109
    include_filter_(include_filter)
110✔
110
{}
110✔
111

112
void RasterData::set_value(size_t y, size_t x, const Particle& p, int level,
275,000✔
113
  Filter* filter, FilterMatch* match)
114
{
115
  // set cell data
116
  if (p.n_coord() <= level) {
275,000!
NEW
117
    id_data_(y, x, 0) = NOT_FOUND;
×
NEW
118
    id_data_(y, x, 1) = NOT_FOUND;
×
119
  } else {
120
    id_data_(y, x, 0) = model::cells.at(p.coord(level).cell())->id_;
275,000!
121
    id_data_(y, x, 1) = level == p.n_coord() - 1
275,000✔
122
                          ? p.cell_instance()
275,000!
NEW
123
                          : cell_instance_at_level(p, level);
×
124
  }
125

126
  // set material data
127
  Cell* c = model::cells.at(p.lowest_coord().cell()).get();
275,000!
128
  if (p.material() == MATERIAL_VOID) {
275,000!
NEW
129
    id_data_(y, x, 2) = MATERIAL_VOID;
×
130
  } else if (c->type_ == Fill::MATERIAL) {
275,000!
131
    Material* m = model::materials.at(p.material()).get();
275,000✔
132
    id_data_(y, x, 2) = m->id_;
275,000✔
133
  }
134

135
  // set filter index (only if filter is being used)
136
  if (include_filter_ && filter) {
275,000!
137
    filter->get_all_bins(p, TallyEstimator::COLLISION, *match);
55,000✔
138
    if (match->bins_.empty()) {
55,000!
NEW
139
      id_data_(y, x, 3) = -1;
×
140
    } else {
141
      id_data_(y, x, 3) = match->bins_[0];
55,000✔
142
    }
143
    match->bins_.clear();
55,000!
144
    match->weights_.clear();
55,000!
145
  }
146

147
  // set temperature (in K)
148
  property_data_(y, x, 0) = (p.sqrtkT() * p.sqrtkT()) / K_BOLTZMANN;
275,000✔
149

150
  // set density (g/cm³)
151
  if (c->type_ != Fill::UNIVERSE && p.material() != MATERIAL_VOID) {
275,000!
152
    Material* m = model::materials.at(p.material()).get();
275,000✔
153
    property_data_(y, x, 1) = m->density_gpcc_;
275,000✔
154
  }
155
}
275,000✔
156

NEW
157
void RasterData::set_overlap(size_t y, size_t x)
×
158
{
159
  // Set cell, instance, and material to OVERLAP, but preserve filter bin
NEW
160
  id_data_(y, x, 0) = OVERLAP;
×
NEW
161
  id_data_(y, x, 1) = OVERLAP;
×
NEW
162
  id_data_(y, x, 2) = OVERLAP;
×
163
  // Note: id_data_(y, x, 3) is NOT overwritten - preserves filter bin for tally
164
  // plotting
165

NEW
166
  property_data_(y, x, 0) = OVERLAP;
×
NEW
167
  property_data_(y, x, 1) = OVERLAP;
×
NEW
168
}
×
169

170
//==============================================================================
171
// Global variables
172
//==============================================================================
173

174
namespace model {
175

176
std::unordered_map<int, int> plot_map;
177
vector<std::unique_ptr<PlottableInterface>> plots;
178
uint64_t plotter_seed = 1;
179

180
} // namespace model
181

182
//==============================================================================
183
// RUN_PLOT controls the logic for making one or many plots
184
//==============================================================================
185

186
extern "C" int openmc_plot_geometry()
121✔
187
{
188

189
  for (auto& pl : model::plots) {
407✔
190
    write_message(5, "Processing plot {}: {}...", pl->id(), pl->path_plot());
286✔
191
    pl->create_output();
286✔
192
  }
193

194
  return 0;
121✔
195
}
196

197
void PlottableInterface::write_image(const ImageData& data) const
231✔
198
{
199
#ifdef USE_LIBPNG
200
  output_png(path_plot(), data);
231✔
201
#else
202
  output_ppm(path_plot(), data);
203
#endif
204
}
231✔
205

206
void Plot::create_output() const
198✔
207
{
208
  if (PlotType::slice == type_) {
198✔
209
    // create 2D image
210
    ImageData image = create_image();
143✔
211
    write_image(image);
143✔
212
  } else if (PlotType::voxel == type_) {
198!
213
    // create voxel file for 3D viewing
214
    create_voxel();
55✔
215
  }
216
}
198✔
217

218
void Plot::print_info() const
154✔
219
{
220
  // Plot type
221
  if (PlotType::slice == type_) {
154✔
222
    fmt::print("Plot Type: Slice\n");
121✔
223
  } else if (PlotType::voxel == type_) {
33!
224
    fmt::print("Plot Type: Voxel\n");
33✔
225
  }
226

227
  // Plot parameters
228
  fmt::print("Origin: {} {} {}\n", origin_[0], origin_[1], origin_[2]);
154✔
229

230
  if (PlotType::slice == type_) {
154✔
231
    fmt::print("Width: {:4} {:4}\n", width_[0], width_[1]);
121✔
232
  } else if (PlotType::voxel == type_) {
33!
233
    fmt::print("Width: {:4} {:4} {:4}\n", width_[0], width_[1], width_[2]);
33✔
234
  }
235

236
  if (PlotColorBy::cells == color_by_) {
154✔
237
    fmt::print("Coloring: Cells\n");
88✔
238
  } else if (PlotColorBy::mats == color_by_) {
66!
239
    fmt::print("Coloring: Materials\n");
66✔
240
  }
241

242
  if (PlotType::slice == type_) {
154✔
243
    switch (basis_) {
121!
244
    case PlotBasis::xy:
77✔
245
      fmt::print("Basis: XY\n");
77✔
246
      break;
77✔
247
    case PlotBasis::xz:
22✔
248
      fmt::print("Basis: XZ\n");
22✔
249
      break;
22✔
250
    case PlotBasis::yz:
22✔
251
      fmt::print("Basis: YZ\n");
22✔
252
      break;
22✔
253
    }
254
    fmt::print("Pixels: {} {}\n", pixels()[0], pixels()[1]);
121✔
255
  } else if (PlotType::voxel == type_) {
33!
256
    fmt::print("Voxels: {} {} {}\n", pixels()[0], pixels()[1], pixels()[2]);
33✔
257
  }
258
}
154✔
259

260
void read_plots_xml()
1,339✔
261
{
262
  // Check if plots.xml exists; this is only necessary when the plot runmode is
263
  // initiated. Otherwise, we want to read plots.xml because it may be called
264
  // later via the API. In that case, its ok for a plots.xml to not exist
265
  std::string filename = settings::path_input + "plots.xml";
1,339✔
266
  if (!file_exists(filename) && settings::run_mode == RunMode::PLOTTING) {
1,339!
267
    fatal_error(fmt::format("Plots XML file '{}' does not exist!", filename));
×
268
  }
269

270
  write_message("Reading plot XML file...", 5);
1,339✔
271

272
  // Parse plots.xml file
273
  pugi::xml_document doc;
1,339✔
274
  doc.load_file(filename.c_str());
1,339✔
275

276
  pugi::xml_node root = doc.document_element();
1,339✔
277

278
  read_plots_xml(root);
1,339✔
279
}
1,339✔
280

281
void read_plots_xml(pugi::xml_node root)
1,803✔
282
{
283
  for (auto node : root.children("plot")) {
2,751✔
284
    std::string plot_desc = "<auto>";
957✔
285
    if (check_for_node(node, "id")) {
957!
286
      plot_desc = get_node_value(node, "id", true);
957✔
287
    }
288

289
    if (check_for_node(node, "type")) {
957!
290
      std::string type_str = get_node_value(node, "type", true);
957✔
291
      if (type_str == "slice") {
957✔
292
        model::plots.emplace_back(
805✔
293
          std::make_unique<Plot>(node, Plot::PlotType::slice));
1,619✔
294
      } else if (type_str == "voxel") {
143✔
295
        model::plots.emplace_back(
55✔
296
          std::make_unique<Plot>(node, Plot::PlotType::voxel));
110✔
297
      } else if (type_str == "wireframe_raytrace") {
88✔
298
        model::plots.emplace_back(
55✔
299
          std::make_unique<WireframeRayTracePlot>(node));
110✔
300
      } else if (type_str == "solid_raytrace") {
33!
301
        model::plots.emplace_back(std::make_unique<SolidRayTracePlot>(node));
33✔
302
      } else {
303
        fatal_error(fmt::format(
×
304
          "Unsupported plot type '{}' in plot {}", type_str, plot_desc));
305
      }
306
      model::plot_map[model::plots.back()->id()] = model::plots.size() - 1;
948✔
307
    } else {
948✔
308
      fatal_error(fmt::format("Must specify plot type in plot {}", plot_desc));
×
309
    }
310
  }
948✔
311
}
1,794✔
312

313
void free_memory_plot()
8,349✔
314
{
315
  model::plots.clear();
8,349✔
316
  model::plot_map.clear();
8,349✔
317
}
8,349✔
318

319
// creates an image based on user input from a plots.xml <plot>
320
// specification in the PNG/PPM format
321
ImageData Plot::create_image() const
143✔
322
{
323
  size_t width = pixels()[0];
143✔
324
  size_t height = pixels()[1];
143✔
325

326
  ImageData data({width, height}, not_found_);
143✔
327

328
  // generate ids for the plot
329
  auto ids = get_map<IdData>();
143✔
330

331
  // assign colors
332
  for (size_t y = 0; y < height; y++) {
30,063✔
333
    for (size_t x = 0; x < width; x++) {
7,622,120✔
334
      int idx = color_by_ == PlotColorBy::cells ? 0 : 2;
7,592,200✔
335
      auto id = ids.data_(y, x, idx);
7,592,200✔
336
      // no setting needed if not found
337
      if (id == NOT_FOUND) {
7,592,200✔
338
        continue;
1,082,532✔
339
      }
340
      if (id == OVERLAP) {
6,537,916✔
341
        data(x, y) = overlap_color_;
28,248✔
342
        continue;
28,248✔
343
      }
344
      if (PlotColorBy::cells == color_by_) {
6,509,668✔
345
        data(x, y) = colors_[model::cell_map[id]];
3,011,668✔
346
      } else if (PlotColorBy::mats == color_by_) {
3,498,000!
347
        if (id == MATERIAL_VOID) {
3,498,000!
348
          data(x, y) = WHITE;
×
349
          continue;
×
350
        }
351
        data(x, y) = colors_[model::material_map[id]];
3,498,000✔
352
      } // color_by if-else
353
    }
354
  }
355

356
  // draw mesh lines if present
357
  if (index_meshlines_mesh_ >= 0) {
143✔
358
    draw_mesh_lines(data);
33✔
359
  }
360

361
  return data;
143✔
362
}
143✔
363

364
void PlottableInterface::set_id(pugi::xml_node plot_node)
957✔
365
{
366
  int id {C_NONE};
957✔
367
  if (check_for_node(plot_node, "id")) {
957!
368
    id = std::stoi(get_node_value(plot_node, "id"));
957✔
369
  }
370

371
  try {
957✔
372
    set_id(id);
957✔
373
  } catch (const std::runtime_error& e) {
×
374
    fatal_error(e.what());
×
375
  }
×
376
}
957✔
377

378
void PlottableInterface::set_id(int id)
968✔
379
{
380
  if (id < 0 && id != C_NONE) {
968!
381
    throw std::runtime_error {fmt::format("Invalid plot ID: {}", id)};
×
382
  }
383

384
  if (id == C_NONE) {
968✔
385
    id = 1;
11✔
386
    for (const auto& p : model::plots) {
22✔
387
      id = std::max(id, p->id() + 1);
22!
388
    }
389
  }
390

391
  if (id_ == id)
968!
392
    return;
393

394
  // Check to make sure this ID doesn't already exist
395
  if (model::plot_map.find(id) != model::plot_map.end()) {
968!
396
    throw std::runtime_error {
×
397
      fmt::format("Two or more plots use the same unique ID: {}", id)};
×
398
  }
399

400
  id_ = id;
968✔
401
}
402

403
// Checks if png or ppm is already present
404
bool file_extension_present(
948✔
405
  const std::string& filename, const std::string& extension)
406
{
407
  std::string file_extension_if_present =
948✔
408
    filename.substr(filename.find_last_of(".") + 1);
948✔
409
  if (file_extension_if_present == extension)
948✔
410
    return true;
55✔
411
  return false;
412
}
948✔
413

414
void Plot::set_output_path(pugi::xml_node plot_node)
869✔
415
{
416
  // Set output file path
417
  std::string filename;
869✔
418

419
  if (check_for_node(plot_node, "filename")) {
869✔
420
    filename = get_node_value(plot_node, "filename");
242✔
421
  } else {
422
    filename = fmt::format("plot_{}", id());
627✔
423
  }
424
  const std::string dir_if_present =
869✔
425
    filename.substr(0, filename.find_last_of("/") + 1);
869✔
426
  if (dir_if_present.size() > 0 && !dir_exists(dir_if_present)) {
869✔
427
    fatal_error(fmt::format("Directory '{}' does not exist!", dir_if_present));
9✔
428
  }
429
  // add appropriate file extension to name
430
  switch (type_) {
860!
431
  case PlotType::slice:
805✔
432
#ifdef USE_LIBPNG
433
    if (!file_extension_present(filename, "png"))
805!
434
      filename.append(".png");
805✔
435
#else
436
    if (!file_extension_present(filename, "ppm"))
437
      filename.append(".ppm");
438
#endif
439
    break;
440
  case PlotType::voxel:
55✔
441
    if (!file_extension_present(filename, "h5"))
55!
442
      filename.append(".h5");
55✔
443
    break;
444
  }
445

446
  path_plot_ = filename;
860✔
447

448
  // Copy plot pixel size
449
  vector<int> pxls = get_node_array<int>(plot_node, "pixels");
1,720✔
450
  if (PlotType::slice == type_) {
860✔
451
    if (pxls.size() == 2) {
805!
452
      pixels()[0] = pxls[0];
805✔
453
      pixels()[1] = pxls[1];
805✔
454
    } else {
455
      fatal_error(
×
456
        fmt::format("<pixels> must be length 2 in slice plot {}", id()));
×
457
    }
458
  } else if (PlotType::voxel == type_) {
55!
459
    if (pxls.size() == 3) {
55!
460
      pixels()[0] = pxls[0];
55✔
461
      pixels()[1] = pxls[1];
55✔
462
      pixels()[2] = pxls[2];
55✔
463
    } else {
464
      fatal_error(
×
465
        fmt::format("<pixels> must be length 3 in voxel plot {}", id()));
×
466
    }
467
  }
468
}
860✔
469

470
void PlottableInterface::set_bg_color(pugi::xml_node plot_node)
957✔
471
{
472
  // Copy plot background color
473
  if (check_for_node(plot_node, "background")) {
957✔
474
    vector<int> bg_rgb = get_node_array<int>(plot_node, "background");
44✔
475
    if (bg_rgb.size() == 3) {
44!
476
      not_found_ = bg_rgb;
44✔
477
    } else {
478
      fatal_error(fmt::format("Bad background RGB in plot {}", id()));
×
479
    }
480
  }
44✔
481
}
957✔
482

483
void Plot::set_basis(pugi::xml_node plot_node)
860✔
484
{
485
  // Copy plot basis
486
  if (PlotType::slice == type_) {
860✔
487
    std::string pl_basis = "xy";
805✔
488
    if (check_for_node(plot_node, "basis")) {
805!
489
      pl_basis = get_node_value(plot_node, "basis", true);
805✔
490
    }
491
    if ("xy" == pl_basis) {
805✔
492
      basis_ = PlotBasis::xy;
731✔
493
    } else if ("xz" == pl_basis) {
74✔
494
      basis_ = PlotBasis::xz;
22✔
495
    } else if ("yz" == pl_basis) {
52!
496
      basis_ = PlotBasis::yz;
52✔
497
    } else {
498
      fatal_error(
×
499
        fmt::format("Unsupported plot basis '{}' in plot {}", pl_basis, id()));
×
500
    }
501
  }
805✔
502
}
860✔
503

504
void Plot::set_origin(pugi::xml_node plot_node)
860✔
505
{
506
  // Copy plotting origin
507
  auto pl_origin = get_node_array<double>(plot_node, "origin");
860✔
508
  if (pl_origin.size() == 3) {
860!
509
    origin_ = pl_origin;
860✔
510
  } else {
511
    fatal_error(fmt::format("Origin must be length 3 in plot {}", id()));
×
512
  }
513
}
860✔
514

515
void Plot::set_width(pugi::xml_node plot_node)
860✔
516
{
517
  // Copy plotting width
518
  vector<double> pl_width = get_node_array<double>(plot_node, "width");
860✔
519
  if (PlotType::slice == type_) {
860✔
520
    if (pl_width.size() == 2) {
805!
521
      width_.x = pl_width[0];
805✔
522
      width_.y = pl_width[1];
805✔
523
      switch (basis_) {
805!
524
      case PlotBasis::xy:
731✔
525
        u_span_ = {width_.x, 0.0, 0.0};
731✔
526
        v_span_ = {0.0, width_.y, 0.0};
731✔
527
        break;
731✔
528
      case PlotBasis::xz:
22✔
529
        u_span_ = {width_.x, 0.0, 0.0};
22✔
530
        v_span_ = {0.0, 0.0, width_.y};
22✔
531
        break;
22✔
532
      case PlotBasis::yz:
52✔
533
        u_span_ = {0.0, width_.x, 0.0};
52✔
534
        v_span_ = {0.0, 0.0, width_.y};
52✔
535
        break;
52✔
NEW
536
      default:
×
NEW
537
        UNREACHABLE();
×
538
      }
539
    } else {
540
      fatal_error(
×
541
        fmt::format("<width> must be length 2 in slice plot {}", id()));
×
542
    }
543
  } else if (PlotType::voxel == type_) {
55!
544
    if (pl_width.size() == 3) {
55!
545
      pl_width = get_node_array<double>(plot_node, "width");
110✔
546
      width_ = pl_width;
55✔
547
    } else {
548
      fatal_error(
×
549
        fmt::format("<width> must be length 3 in voxel plot {}", id()));
×
550
    }
551
  }
552
}
860✔
553

554
void PlottableInterface::set_universe(pugi::xml_node plot_node)
957✔
555
{
556
  // Copy plot universe level
557
  if (check_for_node(plot_node, "level")) {
957!
558
    level_ = std::stoi(get_node_value(plot_node, "level"));
×
559
    if (level_ < 0) {
×
560
      fatal_error(fmt::format("Bad universe level in plot {}", id()));
×
561
    }
562
  } else {
563
    level_ = PLOT_LEVEL_LOWEST;
957✔
564
  }
565
}
957✔
566

567
void PlottableInterface::set_color_by(pugi::xml_node plot_node)
957✔
568
{
569
  // Copy plot color type
570
  std::string pl_color_by = "cell";
957✔
571
  if (check_for_node(plot_node, "color_by")) {
957✔
572
    pl_color_by = get_node_value(plot_node, "color_by", true);
924✔
573
  }
574
  if ("cell" == pl_color_by) {
957✔
575
    color_by_ = PlotColorBy::cells;
287✔
576
  } else if ("material" == pl_color_by) {
670!
577
    color_by_ = PlotColorBy::mats;
670✔
578
  } else {
579
    fatal_error(fmt::format(
×
580
      "Unsupported plot color type '{}' in plot {}", pl_color_by, id()));
×
581
  }
582
}
957✔
583

584
void PlottableInterface::set_default_colors()
968✔
585
{
586
  // Copy plot color type and initialize all colors randomly
587
  if (PlotColorBy::cells == color_by_) {
968✔
588
    colors_.resize(model::cells.size());
287✔
589
  } else if (PlotColorBy::mats == color_by_) {
681!
590
    colors_.resize(model::materials.size());
681✔
591
  }
592

593
  for (auto& c : colors_) {
4,343✔
594
    c = random_color();
3,375✔
595
    // make sure we don't interfere with some default colors
596
    while (c == RED || c == WHITE) {
3,375!
597
      c = random_color();
×
598
    }
599
  }
600
}
968✔
601

602
void PlottableInterface::set_user_colors(pugi::xml_node plot_node)
957✔
603
{
604
  for (auto cn : plot_node.children("color")) {
1,144✔
605
    // Make sure 3 values are specified for RGB
606
    vector<int> user_rgb = get_node_array<int>(cn, "rgb");
187✔
607
    if (user_rgb.size() != 3) {
187!
608
      fatal_error(fmt::format("Bad RGB in plot {}", id()));
×
609
    }
610
    // Ensure that there is an id for this color specification
611
    int col_id;
187✔
612
    if (check_for_node(cn, "id")) {
187!
613
      col_id = std::stoi(get_node_value(cn, "id"));
374✔
614
    } else {
615
      fatal_error(fmt::format(
×
616
        "Must specify id for color specification in plot {}", id()));
×
617
    }
618
    // Add RGB
619
    if (PlotColorBy::cells == color_by_) {
187✔
620
      if (model::cell_map.find(col_id) != model::cell_map.end()) {
88!
621
        col_id = model::cell_map[col_id];
88✔
622
        colors_[col_id] = user_rgb;
88✔
623
      } else {
624
        warning(fmt::format(
×
625
          "Could not find cell {} specified in plot {}", col_id, id()));
×
626
      }
627
    } else if (PlotColorBy::mats == color_by_) {
99!
628
      if (model::material_map.find(col_id) != model::material_map.end()) {
99!
629
        col_id = model::material_map[col_id];
99✔
630
        colors_[col_id] = user_rgb;
99✔
631
      } else {
632
        warning(fmt::format(
×
633
          "Could not find material {} specified in plot {}", col_id, id()));
×
634
      }
635
    }
636
  } // color node loop
187✔
637
}
957✔
638

639
void Plot::set_meshlines(pugi::xml_node plot_node)
860✔
640
{
641
  // Deal with meshlines
642
  pugi::xpath_node_set mesh_line_nodes = plot_node.select_nodes("meshlines");
860✔
643

644
  if (!mesh_line_nodes.empty()) {
860✔
645
    if (PlotType::voxel == type_) {
33!
646
      warning(fmt::format("Meshlines ignored in voxel plot {}", id()));
×
647
    }
648

649
    if (mesh_line_nodes.size() == 1) {
33!
650
      // Get first meshline node
651
      pugi::xml_node meshlines_node = mesh_line_nodes[0].node();
33✔
652

653
      // Check mesh type
654
      std::string meshtype;
33✔
655
      if (check_for_node(meshlines_node, "meshtype")) {
33!
656
        meshtype = get_node_value(meshlines_node, "meshtype");
33✔
657
      } else {
658
        fatal_error(fmt::format(
×
659
          "Must specify a meshtype for meshlines specification in plot {}",
660
          id()));
×
661
      }
662

663
      // Ensure that there is a linewidth for this meshlines specification
664
      std::string meshline_width;
33✔
665
      if (check_for_node(meshlines_node, "linewidth")) {
33!
666
        meshline_width = get_node_value(meshlines_node, "linewidth");
33✔
667
        meshlines_width_ = std::stoi(meshline_width);
33✔
668
      } else {
669
        fatal_error(fmt::format(
×
670
          "Must specify a linewidth for meshlines specification in plot {}",
671
          id()));
×
672
      }
673

674
      // Check for color
675
      if (check_for_node(meshlines_node, "color")) {
33!
676
        // Check and make sure 3 values are specified for RGB
677
        vector<int> ml_rgb = get_node_array<int>(meshlines_node, "color");
×
678
        if (ml_rgb.size() != 3) {
×
679
          fatal_error(
×
680
            fmt::format("Bad RGB for meshlines color in plot {}", id()));
×
681
        }
682
        meshlines_color_ = ml_rgb;
×
683
      }
684

685
      // Set mesh based on type
686
      if ("ufs" == meshtype) {
33!
687
        if (!simulation::ufs_mesh) {
×
688
          fatal_error(
×
689
            fmt::format("No UFS mesh for meshlines on plot {}", id()));
×
690
        } else {
691
          for (int i = 0; i < model::meshes.size(); ++i) {
×
692
            if (const auto* m =
×
693
                  dynamic_cast<const RegularMesh*>(model::meshes[i].get())) {
×
694
              if (m == simulation::ufs_mesh) {
×
695
                index_meshlines_mesh_ = i;
×
696
              }
697
            }
698
          }
699
          if (index_meshlines_mesh_ == -1)
×
700
            fatal_error("Could not find the UFS mesh for meshlines plot");
×
701
        }
702
      } else if ("entropy" == meshtype) {
33✔
703
        if (!simulation::entropy_mesh) {
22!
704
          fatal_error(
×
705
            fmt::format("No entropy mesh for meshlines on plot {}", id()));
×
706
        } else {
707
          for (int i = 0; i < model::meshes.size(); ++i) {
55✔
708
            if (const auto* m =
66✔
709
                  dynamic_cast<const RegularMesh*>(model::meshes[i].get())) {
55!
710
              if (m == simulation::entropy_mesh) {
22!
711
                index_meshlines_mesh_ = i;
22✔
712
              }
713
            }
714
          }
715
          if (index_meshlines_mesh_ == -1)
22!
716
            fatal_error("Could not find the entropy mesh for meshlines plot");
×
717
        }
718
      } else if ("tally" == meshtype) {
11!
719
        // Ensure that there is a mesh id if the type is tally
720
        int tally_mesh_id;
11✔
721
        if (check_for_node(meshlines_node, "id")) {
11!
722
          tally_mesh_id = std::stoi(get_node_value(meshlines_node, "id"));
22✔
723
        } else {
724
          std::stringstream err_msg;
×
725
          fatal_error(fmt::format("Must specify a mesh id for meshlines tally "
×
726
                                  "mesh specification in plot {}",
727
            id()));
×
728
        }
×
729
        // find the tally index
730
        int idx;
11✔
731
        int err = openmc_get_mesh_index(tally_mesh_id, &idx);
11✔
732
        if (err != 0) {
11!
733
          fatal_error(fmt::format("Could not find mesh {} specified in "
×
734
                                  "meshlines for plot {}",
735
            tally_mesh_id, id()));
×
736
        }
737
        index_meshlines_mesh_ = idx;
11✔
738
      } else {
739
        fatal_error(fmt::format("Invalid type for meshlines on plot {}", id()));
×
740
      }
741
    } else {
33✔
742
      fatal_error(fmt::format("Mutliple meshlines specified in plot {}", id()));
×
743
    }
744
  }
745
}
860✔
746

747
void PlottableInterface::set_mask(pugi::xml_node plot_node)
957✔
748
{
749
  // Deal with masks
750
  pugi::xpath_node_set mask_nodes = plot_node.select_nodes("mask");
957✔
751

752
  if (!mask_nodes.empty()) {
957✔
753
    if (mask_nodes.size() == 1) {
33!
754
      // Get pointer to mask
755
      pugi::xml_node mask_node = mask_nodes[0].node();
33✔
756

757
      // Determine how many components there are and allocate
758
      vector<int> iarray = get_node_array<int>(mask_node, "components");
33✔
759
      if (iarray.size() == 0) {
33!
760
        fatal_error(
×
761
          fmt::format("Missing <components> in mask of plot {}", id()));
×
762
      }
763

764
      // First we need to change the user-specified identifiers to indices
765
      // in the cell and material arrays
766
      for (auto& col_id : iarray) {
99✔
767
        if (PlotColorBy::cells == color_by_) {
66!
768
          if (model::cell_map.find(col_id) != model::cell_map.end()) {
66!
769
            col_id = model::cell_map[col_id];
66✔
770
          } else {
771
            fatal_error(fmt::format("Could not find cell {} specified in the "
×
772
                                    "mask in plot {}",
773
              col_id, id()));
×
774
          }
775
        } else if (PlotColorBy::mats == color_by_) {
×
776
          if (model::material_map.find(col_id) != model::material_map.end()) {
×
777
            col_id = model::material_map[col_id];
×
778
          } else {
779
            fatal_error(fmt::format("Could not find material {} specified in "
×
780
                                    "the mask in plot {}",
781
              col_id, id()));
×
782
          }
783
        }
784
      }
785

786
      // Alter colors based on mask information
787
      for (int j = 0; j < colors_.size(); j++) {
132✔
788
        if (contains(iarray, j)) {
99✔
789
          if (check_for_node(mask_node, "background")) {
66!
790
            vector<int> bg_rgb = get_node_array<int>(mask_node, "background");
66✔
791
            colors_[j] = bg_rgb;
66✔
792
          } else {
66✔
793
            colors_[j] = WHITE;
×
794
          }
795
        }
796
      }
797

798
    } else {
33✔
799
      fatal_error(fmt::format("Mutliple masks specified in plot {}", id()));
×
800
    }
801
  }
802
}
957✔
803

804
void PlottableInterface::set_overlap_color(pugi::xml_node plot_node)
957✔
805
{
806
  color_overlaps_ = false;
957✔
807
  if (check_for_node(plot_node, "show_overlaps")) {
957✔
808
    color_overlaps_ = get_node_value_bool(plot_node, "show_overlaps");
22✔
809
    // check for custom overlap color
810
    if (check_for_node(plot_node, "overlap_color")) {
22✔
811
      if (!color_overlaps_) {
11!
812
        warning(fmt::format(
×
813
          "Overlap color specified in plot {} but overlaps won't be shown.",
814
          id()));
×
815
      }
816
      vector<int> olap_clr = get_node_array<int>(plot_node, "overlap_color");
11✔
817
      if (olap_clr.size() == 3) {
11!
818
        overlap_color_ = olap_clr;
11✔
819
      } else {
820
        fatal_error(fmt::format("Bad overlap RGB in plot {}", id()));
×
821
      }
822
    }
11✔
823
  }
824

825
  // make sure we allocate the vector for counting overlap checks if
826
  // they're going to be plotted
827
  if (color_overlaps_ && settings::run_mode == RunMode::PLOTTING) {
957!
828
    settings::check_overlaps = true;
22✔
829
    model::overlap_check_count.resize(model::cells.size(), 0);
22✔
830
  }
831
}
957✔
832

833
PlottableInterface::PlottableInterface(pugi::xml_node plot_node)
957✔
834
{
835
  set_id(plot_node);
957✔
836
  set_bg_color(plot_node);
957✔
837
  set_universe(plot_node);
957✔
838
  set_color_by(plot_node);
957✔
839
  set_default_colors();
957✔
840
  set_user_colors(plot_node);
957✔
841
  set_mask(plot_node);
957✔
842
  set_overlap_color(plot_node);
957✔
843
}
957✔
844

845
Plot::Plot(pugi::xml_node plot_node, PlotType type)
869✔
846
  : PlottableInterface(plot_node), type_(type), index_meshlines_mesh_ {-1}
869✔
847
{
848
  set_output_path(plot_node);
869✔
849
  set_basis(plot_node);
860✔
850
  set_origin(plot_node);
860✔
851
  set_width(plot_node);
860✔
852
  set_meshlines(plot_node);
860✔
853
  slice_level_ = level_; // Copy level employed in SlicePlotBase::get_map
860✔
854
  slice_color_overlaps_ = color_overlaps_;
860✔
855
}
860✔
856

857
//==============================================================================
858
// OUTPUT_PPM writes out a previously generated image to a PPM file
859
//==============================================================================
860

861
void output_ppm(const std::string& filename, const ImageData& data)
×
862
{
863
  // Open PPM file for writing
864
  std::string fname = filename;
×
865
  fname = strtrim(fname);
×
866
  std::ofstream of;
×
867

868
  of.open(fname);
×
869

870
  // Write header
871
  of << "P6\n";
×
872
  of << data.shape(0) << " " << data.shape(1) << "\n";
×
873
  of << "255\n";
×
874
  of.close();
×
875

876
  of.open(fname, std::ios::binary | std::ios::app);
×
877
  // Write color for each pixel
878
  for (int y = 0; y < data.shape(1); y++) {
×
879
    for (int x = 0; x < data.shape(0); x++) {
×
880
      RGBColor rgb = data(x, y);
×
881
      of << rgb.red << rgb.green << rgb.blue;
×
882
    }
883
  }
884
  of << "\n";
×
885
}
×
886

887
//==============================================================================
888
// OUTPUT_PNG writes out a previously generated image to a PNG file
889
//==============================================================================
890

891
#ifdef USE_LIBPNG
892
void output_png(const std::string& filename, const ImageData& data)
231✔
893
{
894
  // Open PNG file for writing
895
  std::string fname = filename;
231✔
896
  fname = strtrim(fname);
231✔
897
  auto fp = std::fopen(fname.c_str(), "wb");
231✔
898

899
  // Initialize write and info structures
900
  auto png_ptr =
231✔
901
    png_create_write_struct(PNG_LIBPNG_VER_STRING, nullptr, nullptr, nullptr);
231✔
902
  auto info_ptr = png_create_info_struct(png_ptr);
231✔
903

904
  // Setup exception handling
905
  if (setjmp(png_jmpbuf(png_ptr)))
231!
906
    fatal_error("Error during png creation");
×
907

908
  png_init_io(png_ptr, fp);
231✔
909

910
  // Write header (8 bit colour depth)
911
  int width = data.shape(0);
231!
912
  int height = data.shape(1);
231!
913
  png_set_IHDR(png_ptr, info_ptr, width, height, 8, PNG_COLOR_TYPE_RGB,
231✔
914
    PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE);
915
  png_write_info(png_ptr, info_ptr);
231✔
916

917
  // Allocate memory for one row (3 bytes per pixel - RGB)
918
  std::vector<png_byte> row(3 * width);
231✔
919

920
  // Write color for each pixel
921
  for (int y = 0; y < height; y++) {
47,751✔
922
    for (int x = 0; x < width; x++) {
11,159,720✔
923
      RGBColor rgb = data(x, y);
11,112,200✔
924
      row[3 * x] = rgb.red;
11,112,200✔
925
      row[3 * x + 1] = rgb.green;
11,112,200✔
926
      row[3 * x + 2] = rgb.blue;
11,112,200✔
927
    }
928
    png_write_row(png_ptr, row.data());
47,520✔
929
  }
930

931
  // End write
932
  png_write_end(png_ptr, nullptr);
231✔
933

934
  // Clean up data structures
935
  std::fclose(fp);
231✔
936
  png_free_data(png_ptr, info_ptr, PNG_FREE_ALL, -1);
231✔
937
  png_destroy_write_struct(&png_ptr, &info_ptr);
231✔
938
}
231✔
939
#endif
940

941
//==============================================================================
942
// DRAW_MESH_LINES draws mesh line boundaries on an image
943
//==============================================================================
944

945
void Plot::draw_mesh_lines(ImageData& data) const
33✔
946
{
947
  RGBColor rgb;
33!
948
  rgb = meshlines_color_;
33✔
949

950
  int ax1, ax2;
33✔
951
  switch (basis_) {
33!
952
  case PlotBasis::xy:
953
    ax1 = 0;
954
    ax2 = 1;
955
    break;
956
  case PlotBasis::xz:
957
    ax1 = 0;
958
    ax2 = 2;
959
    break;
960
  case PlotBasis::yz:
961
    ax1 = 1;
962
    ax2 = 2;
963
    break;
964
  default:
×
965
    UNREACHABLE();
×
966
  }
967

968
  Position ll_plot {origin_};
33✔
969
  Position ur_plot {origin_};
33✔
970

971
  ll_plot[ax1] -= width_[0] / 2.;
33✔
972
  ll_plot[ax2] -= width_[1] / 2.;
33✔
973
  ur_plot[ax1] += width_[0] / 2.;
33✔
974
  ur_plot[ax2] += width_[1] / 2.;
33✔
975

976
  Position width = ur_plot - ll_plot;
33✔
977

978
  // Find the (axis-aligned) lines of the mesh that intersect this plot.
979
  auto axis_lines =
33✔
980
    model::meshes[index_meshlines_mesh_]->plot(ll_plot, ur_plot);
33✔
981

982
  // Find the bounds along the second axis (accounting for low-D meshes).
983
  int ax2_min, ax2_max;
33✔
984
  if (axis_lines.second.size() > 0) {
33!
985
    double frac = (axis_lines.second.back() - ll_plot[ax2]) / width[ax2];
33✔
986
    ax2_min = (1.0 - frac) * pixels()[1];
33✔
987
    if (ax2_min < 0)
33✔
988
      ax2_min = 0;
989
    frac = (axis_lines.second.front() - ll_plot[ax2]) / width[ax2];
33✔
990
    ax2_max = (1.0 - frac) * pixels()[1];
33!
991
    if (ax2_max > pixels()[1])
33!
992
      ax2_max = pixels()[1];
×
993
  } else {
994
    ax2_min = 0;
×
995
    ax2_max = pixels()[1];
×
996
  }
997

998
  // Iterate across the first axis and draw lines.
999
  for (auto ax1_val : axis_lines.first) {
187✔
1000
    double frac = (ax1_val - ll_plot[ax1]) / width[ax1];
154✔
1001
    int ax1_ind = frac * pixels()[0];
154✔
1002
    for (int ax2_ind = ax2_min; ax2_ind < ax2_max; ++ax2_ind) {
24,948✔
1003
      for (int plus = 0; plus <= meshlines_width_; plus++) {
49,588✔
1004
        if (ax1_ind + plus >= 0 && ax1_ind + plus < pixels()[0])
24,794!
1005
          data(ax1_ind + plus, ax2_ind) = rgb;
24,794✔
1006
        if (ax1_ind - plus >= 0 && ax1_ind - plus < pixels()[0])
24,794!
1007
          data(ax1_ind - plus, ax2_ind) = rgb;
24,794✔
1008
      }
1009
    }
1010
  }
1011

1012
  // Find the bounds along the first axis.
1013
  int ax1_min, ax1_max;
33✔
1014
  if (axis_lines.first.size() > 0) {
33!
1015
    double frac = (axis_lines.first.front() - ll_plot[ax1]) / width[ax1];
33✔
1016
    ax1_min = frac * pixels()[0];
33✔
1017
    if (ax1_min < 0)
33✔
1018
      ax1_min = 0;
1019
    frac = (axis_lines.first.back() - ll_plot[ax1]) / width[ax1];
33✔
1020
    ax1_max = frac * pixels()[0];
33!
1021
    if (ax1_max > pixels()[0])
33!
1022
      ax1_max = pixels()[0];
×
1023
  } else {
1024
    ax1_min = 0;
×
1025
    ax1_max = pixels()[0];
×
1026
  }
1027

1028
  // Iterate across the second axis and draw lines.
1029
  for (auto ax2_val : axis_lines.second) {
209✔
1030
    double frac = (ax2_val - ll_plot[ax2]) / width[ax2];
176✔
1031
    int ax2_ind = (1.0 - frac) * pixels()[1];
176✔
1032
    for (int ax1_ind = ax1_min; ax1_ind < ax1_max; ++ax1_ind) {
28,336✔
1033
      for (int plus = 0; plus <= meshlines_width_; plus++) {
56,320✔
1034
        if (ax2_ind + plus >= 0 && ax2_ind + plus < pixels()[1])
28,160!
1035
          data(ax1_ind, ax2_ind + plus) = rgb;
28,160✔
1036
        if (ax2_ind - plus >= 0 && ax2_ind - plus < pixels()[1])
28,160!
1037
          data(ax1_ind, ax2_ind - plus) = rgb;
28,160✔
1038
      }
1039
    }
1040
  }
1041
}
33✔
1042

1043
/* outputs a binary file that can be input into silomesh for 3D geometry
1044
 * visualization.  It works the same way as create_image by dragging a particle
1045
 * across the geometry for the specified number of voxels. The first 3 int's in
1046
 * the binary are the number of x, y, and z voxels.  The next 3 double's are
1047
 * the widths of the voxels in the x, y, and z directions. The next 3 double's
1048
 * are the x, y, and z coordinates of the lower left point. Finally the binary
1049
 * is filled with entries of four int's each. Each 'row' in the binary contains
1050
 * four int's: 3 for x,y,z position and 1 for cell or material id.  For 1
1051
 * million voxels this produces a file of approximately 15MB.
1052
 */
1053
void Plot::create_voxel() const
55✔
1054
{
1055
  // compute voxel widths in each direction
1056
  array<double, 3> vox;
55✔
1057
  vox[0] = width_[0] / static_cast<double>(pixels()[0]);
55✔
1058
  vox[1] = width_[1] / static_cast<double>(pixels()[1]);
55✔
1059
  vox[2] = width_[2] / static_cast<double>(pixels()[2]);
55✔
1060

1061
  // initial particle position
1062
  Position ll = origin_ - width_ / 2.;
55✔
1063

1064
  // Open binary plot file for writing
1065
  std::ofstream of;
55✔
1066
  std::string fname = std::string(path_plot_);
55✔
1067
  fname = strtrim(fname);
55✔
1068
  hid_t file_id = file_open(fname, 'w');
55✔
1069

1070
  // write header info
1071
  write_attribute(file_id, "filetype", "voxel");
55✔
1072
  write_attribute(file_id, "version", VERSION_VOXEL);
55✔
1073
  write_attribute(file_id, "openmc_version", VERSION);
55✔
1074

1075
#ifdef GIT_SHA1
1076
  write_attribute(file_id, "git_sha1", GIT_SHA1);
1077
#endif
1078

1079
  // Write current date and time
1080
  write_attribute(file_id, "date_and_time", time_stamp().c_str());
110✔
1081
  array<int, 3> h5_pixels;
55✔
1082
  std::copy(pixels().begin(), pixels().end(), h5_pixels.begin());
55✔
1083
  write_attribute(file_id, "num_voxels", h5_pixels);
55✔
1084
  write_attribute(file_id, "voxel_width", vox);
55✔
1085
  write_attribute(file_id, "lower_left", ll);
55✔
1086

1087
  // Create dataset for voxel data -- note that the dimensions are reversed
1088
  // since we want the order in the file to be z, y, x
1089
  hsize_t dims[3];
55✔
1090
  dims[0] = pixels()[2];
55✔
1091
  dims[1] = pixels()[1];
55✔
1092
  dims[2] = pixels()[0];
55✔
1093
  hid_t dspace, dset, memspace;
55✔
1094
  voxel_init(file_id, &(dims[0]), &dspace, &dset, &memspace);
55✔
1095

1096
  SlicePlotBase pltbase;
55✔
1097
  pltbase.width_ = width_;
55✔
1098
  pltbase.origin_ = origin_;
55✔
1099
  pltbase.basis_ = PlotBasis::xy;
55✔
1100
  pltbase.u_span_ = {width_.x, 0.0, 0.0};
55✔
1101
  pltbase.v_span_ = {0.0, width_.y, 0.0};
55✔
1102
  pltbase.pixels() = pixels();
55✔
1103
  pltbase.slice_color_overlaps_ = color_overlaps_;
55✔
1104

1105
  ProgressBar pb;
55✔
1106
  for (int z = 0; z < pixels()[2]; z++) {
4,785✔
1107
    // update z coordinate
1108
    pltbase.origin_.z = ll.z + z * vox[2];
4,730✔
1109

1110
    // generate ids using plotbase
1111
    IdData ids = pltbase.get_map<IdData>();
4,730✔
1112

1113
    // select only cell/material ID data and flip the y-axis
1114
    int idx = color_by_ == PlotColorBy::cells ? 0 : 2;
4,730!
1115
    // Extract 2D slice at index idx from 3D data
1116
    size_t rows = ids.data_.shape(0);
4,730!
1117
    size_t cols = ids.data_.shape(1);
4,730!
1118
    tensor::Tensor<int32_t> data_slice({rows, cols});
4,730✔
1119
    for (size_t r = 0; r < rows; ++r)
912,230✔
1120
      for (size_t c = 0; c < cols; ++c)
179,382,500✔
1121
        data_slice(r, c) = ids.data_(r, c, idx);
178,475,000✔
1122
    tensor::Tensor<int32_t> data_flipped = data_slice.flip(0);
4,730✔
1123

1124
    // Write to HDF5 dataset
1125
    voxel_write_slice(z, dspace, dset, memspace, data_flipped.data());
4,730✔
1126

1127
    // update progress bar
1128
    pb.set_value(
4,730✔
1129
      100. * static_cast<double>(z + 1) / static_cast<double>((pixels()[2])));
4,730✔
1130
  }
14,190✔
1131

1132
  voxel_finalize(dspace, dset, memspace);
55✔
1133
  file_close(file_id);
55✔
1134
}
55✔
1135

1136
void voxel_init(hid_t file_id, const hsize_t* dims, hid_t* dspace, hid_t* dset,
55✔
1137
  hid_t* memspace)
1138
{
1139
  // Create dataspace/dataset for voxel data
1140
  *dspace = H5Screate_simple(3, dims, nullptr);
55✔
1141
  *dset = H5Dcreate(file_id, "data", H5T_NATIVE_INT, *dspace, H5P_DEFAULT,
55✔
1142
    H5P_DEFAULT, H5P_DEFAULT);
1143

1144
  // Create dataspace for a slice of the voxel
1145
  hsize_t dims_slice[2] {dims[1], dims[2]};
55✔
1146
  *memspace = H5Screate_simple(2, dims_slice, nullptr);
55✔
1147

1148
  // Select hyperslab in dataspace
1149
  hsize_t start[3] {0, 0, 0};
55✔
1150
  hsize_t count[3] {1, dims[1], dims[2]};
55✔
1151
  H5Sselect_hyperslab(*dspace, H5S_SELECT_SET, start, nullptr, count, nullptr);
55✔
1152
}
55✔
1153

1154
void voxel_write_slice(
4,730✔
1155
  int x, hid_t dspace, hid_t dset, hid_t memspace, void* buf)
1156
{
1157
  hssize_t offset[3] {x, 0, 0};
4,730✔
1158
  H5Soffset_simple(dspace, offset);
4,730✔
1159
  H5Dwrite(dset, H5T_NATIVE_INT, memspace, dspace, H5P_DEFAULT, buf);
4,730✔
1160
}
4,730✔
1161

1162
void voxel_finalize(hid_t dspace, hid_t dset, hid_t memspace)
55✔
1163
{
1164
  H5Dclose(dset);
55✔
1165
  H5Sclose(dspace);
55✔
1166
  H5Sclose(memspace);
55✔
1167
}
55✔
1168

1169
RGBColor random_color(void)
3,375✔
1170
{
1171
  return {int(prn(&model::plotter_seed) * 255),
3,375✔
1172
    int(prn(&model::plotter_seed) * 255), int(prn(&model::plotter_seed) * 255)};
3,375✔
1173
}
1174

1175
RayTracePlot::RayTracePlot(pugi::xml_node node) : PlottableInterface(node)
88✔
1176
{
1177
  set_look_at(node);
88✔
1178
  set_camera_position(node);
88✔
1179
  set_field_of_view(node);
88✔
1180
  set_pixels(node);
88✔
1181
  set_orthographic_width(node);
88✔
1182
  set_output_path(node);
88✔
1183

1184
  if (check_for_node(node, "orthographic_width") &&
99!
1185
      check_for_node(node, "field_of_view"))
11✔
1186
    fatal_error("orthographic_width and field_of_view are mutually exclusive "
×
1187
                "parameters.");
1188
}
88✔
1189

1190
void RayTracePlot::update_view()
110✔
1191
{
1192
  // Get centerline vector for camera-to-model. We create vectors around this
1193
  // that form a pixel array, and then trace rays along that.
1194
  auto up = up_ / up_.norm();
110✔
1195
  Direction looking_direction = look_at_ - camera_position_;
110✔
1196
  looking_direction /= looking_direction.norm();
110✔
1197
  if (std::abs(std::abs(looking_direction.dot(up)) - 1.0) < 1e-9)
110!
1198
    fatal_error("Up vector cannot align with vector between camera position "
×
1199
                "and look_at!");
1200
  Direction cam_yaxis = looking_direction.cross(up);
110✔
1201
  cam_yaxis /= cam_yaxis.norm();
110✔
1202
  Direction cam_zaxis = cam_yaxis.cross(looking_direction);
110✔
1203
  cam_zaxis /= cam_zaxis.norm();
110✔
1204

1205
  // Cache the camera-to-model matrix
1206
  camera_to_model_ = {looking_direction.x, cam_yaxis.x, cam_zaxis.x,
110✔
1207
    looking_direction.y, cam_yaxis.y, cam_zaxis.y, looking_direction.z,
110✔
1208
    cam_yaxis.z, cam_zaxis.z};
110✔
1209
}
110✔
1210

1211
WireframeRayTracePlot::WireframeRayTracePlot(pugi::xml_node node)
55✔
1212
  : RayTracePlot(node)
55✔
1213
{
1214
  set_opacities(node);
55✔
1215
  set_wireframe_thickness(node);
55✔
1216
  set_wireframe_ids(node);
55✔
1217
  set_wireframe_color(node);
55✔
1218
  update_view();
55✔
1219
}
55✔
1220

1221
void WireframeRayTracePlot::set_wireframe_color(pugi::xml_node plot_node)
55✔
1222
{
1223
  // Copy plot wireframe color
1224
  if (check_for_node(plot_node, "wireframe_color")) {
55!
1225
    vector<int> w_rgb = get_node_array<int>(plot_node, "wireframe_color");
×
1226
    if (w_rgb.size() == 3) {
×
1227
      wireframe_color_ = w_rgb;
×
1228
    } else {
1229
      fatal_error(fmt::format("Bad wireframe RGB in plot {}", id()));
×
1230
    }
1231
  }
×
1232
}
55✔
1233

1234
void RayTracePlot::set_output_path(pugi::xml_node node)
88✔
1235
{
1236
  // Set output file path
1237
  std::string filename;
88✔
1238

1239
  if (check_for_node(node, "filename")) {
88✔
1240
    filename = get_node_value(node, "filename");
77✔
1241
  } else {
1242
    filename = fmt::format("plot_{}", id());
11✔
1243
  }
1244

1245
#ifdef USE_LIBPNG
1246
  if (!file_extension_present(filename, "png"))
88✔
1247
    filename.append(".png");
33✔
1248
#else
1249
  if (!file_extension_present(filename, "ppm"))
1250
    filename.append(".ppm");
1251
#endif
1252
  path_plot_ = filename;
176✔
1253
}
88✔
1254

1255
bool WireframeRayTracePlot::trackstack_equivalent(
3,041,159✔
1256
  const std::vector<TrackSegment>& track1,
1257
  const std::vector<TrackSegment>& track2) const
1258
{
1259
  if (wireframe_ids_.empty()) {
3,041,159✔
1260
    // Draw wireframe for all surfaces/cells/materials
1261
    if (track1.size() != track2.size())
2,545,070✔
1262
      return false;
1263
    for (int i = 0; i < track1.size(); ++i) {
6,707,954✔
1264
      if (track1[i].id != track2[i].id ||
4,236,771✔
1265
          track1[i].surface_index != track2[i].surface_index) {
4,236,639✔
1266
        return false;
1267
      }
1268
    }
1269
    return true;
1270
  } else {
1271
    // This runs in O(nm) where n is the intersection stack size
1272
    // and m is the number of IDs we are wireframing. A simpler
1273
    // algorithm can likely be found.
1274
    for (const int id : wireframe_ids_) {
986,194✔
1275
      int t1_i = 0;
496,089✔
1276
      int t2_i = 0;
496,089✔
1277

1278
      // Advance to first instance of the ID
1279
      while (t1_i < track1.size() && t2_i < track2.size()) {
562,430✔
1280
        while (t1_i < track1.size() && track1[t1_i].id != id)
392,832✔
1281
          t1_i++;
229,053✔
1282
        while (t2_i < track2.size() && track2[t2_i].id != id)
393,668✔
1283
          t2_i++;
229,889✔
1284

1285
        // This one is really important!
1286
        if ((t1_i == track1.size() && t2_i != track2.size()) ||
163,779✔
1287
            (t1_i != track1.size() && t2_i == track2.size()))
162,096✔
1288
          return false;
3,718✔
1289
        if (t1_i == track1.size() && t2_i == track2.size())
160,061!
1290
          break;
1291
        // Check if surface different
1292
        if (track1[t1_i].surface_index != track2[t2_i].surface_index)
68,607✔
1293
          return false;
1294

1295
        // Pretty sure this should not be used:
1296
        // if (t2_i != track2.size() - 1 &&
1297
        //     t1_i != track1.size() - 1 &&
1298
        //     track1[t1_i+1].id != track2[t2_i+1].id) return false;
1299
        if (t2_i != 0 && t1_i != 0 &&
67,122✔
1300
            track1[t1_i - 1].surface_index != track2[t2_i - 1].surface_index)
53,944✔
1301
          return false;
1302

1303
        // Check if neighboring cells are different
1304
        // if (track1[t1_i ? t1_i - 1 : 0].id != track2[t2_i ? t2_i - 1 : 0].id)
1305
        // return false; if (track1[t1_i < track1.size() - 1 ? t1_i + 1 : t1_i
1306
        // ].id !=
1307
        //    track2[t2_i < track2.size() - 1 ? t2_i + 1 : t2_i].id) return
1308
        //    false;
1309
        t1_i++, t2_i++;
66,341✔
1310
      }
1311
    }
1312
    return true;
1313
  }
1314
}
1315

1316
std::pair<Position, Direction> RayTracePlot::get_pixel_ray(
3,521,056✔
1317
  int horiz, int vert) const
1318
{
1319
  // Compute field of view in radians
1320
  constexpr double DEGREE_TO_RADIAN = M_PI / 180.0;
3,521,056✔
1321
  double horiz_fov_radians = horizontal_field_of_view_ * DEGREE_TO_RADIAN;
3,521,056✔
1322
  double p0 = static_cast<double>(pixels()[0]);
3,521,056✔
1323
  double p1 = static_cast<double>(pixels()[1]);
3,521,056✔
1324
  double vert_fov_radians = horiz_fov_radians * p1 / p0;
3,521,056✔
1325

1326
  // focal_plane_dist can be changed to alter the perspective distortion
1327
  // effect. This is in units of cm. This seems to look good most of the
1328
  // time. TODO let this variable be set through XML.
1329
  constexpr double focal_plane_dist = 10.0;
3,521,056✔
1330
  const double dx = 2.0 * focal_plane_dist * std::tan(0.5 * horiz_fov_radians);
3,521,056✔
1331
  const double dy = p1 / p0 * dx;
3,521,056✔
1332

1333
  std::pair<Position, Direction> result;
3,521,056✔
1334

1335
  // Generate the starting position/direction of the ray
1336
  if (orthographic_width_ == C_NONE) { // perspective projection
3,521,056✔
1337
    Direction camera_local_vec;
3,081,056✔
1338
    camera_local_vec.x = focal_plane_dist;
3,081,056✔
1339
    camera_local_vec.y = -0.5 * dx + horiz * dx / p0;
3,081,056✔
1340
    camera_local_vec.z = 0.5 * dy - vert * dy / p1;
3,081,056✔
1341
    camera_local_vec /= camera_local_vec.norm();
3,081,056✔
1342

1343
    result.first = camera_position_;
3,081,056✔
1344
    result.second = camera_local_vec.rotate(camera_to_model_);
3,081,056✔
1345
  } else { // orthographic projection
1346

1347
    double x_pix_coord = (static_cast<double>(horiz) - p0 / 2.0) / p0;
440,000✔
1348
    double y_pix_coord = (static_cast<double>(vert) - p1 / 2.0) / p1;
440,000✔
1349

1350
    result.first = camera_position_ +
440,000✔
1351
                   camera_y_axis() * x_pix_coord * orthographic_width_ +
440,000✔
1352
                   camera_z_axis() * y_pix_coord * orthographic_width_;
440,000✔
1353
    result.second = camera_x_axis();
440,000✔
1354
  }
1355

1356
  return result;
3,521,056✔
1357
}
1358

1359
ImageData WireframeRayTracePlot::create_image() const
55✔
1360
{
1361
  size_t width = pixels()[0];
55✔
1362
  size_t height = pixels()[1];
55✔
1363
  ImageData data({width, height}, not_found_);
55✔
1364

1365
  // This array marks where the initial wireframe was drawn. We convolve it with
1366
  // a filter that gets adjusted with the wireframe thickness in order to
1367
  // thicken the lines.
1368
  tensor::Tensor<int> wireframe_initial(
55✔
1369
    {static_cast<size_t>(width), static_cast<size_t>(height)}, 0);
55✔
1370

1371
  /* Holds all of the track segments for the current rendered line of pixels.
1372
   * old_segments holds a copy of this_line_segments from the previous line.
1373
   * By holding both we can check if the cell/material intersection stack
1374
   * differs from the left or upper neighbor. This allows a robustly drawn
1375
   * wireframe. If only checking the left pixel (which requires substantially
1376
   * less memory), the wireframe tends to be spotty and be disconnected for
1377
   * surface edges oriented horizontally in the rendering.
1378
   *
1379
   * Note that a vector of vectors is required rather than a 2-tensor,
1380
   * since the stack size varies within each column.
1381
   */
1382
  const int n_threads = num_threads();
55✔
1383
  std::vector<std::vector<std::vector<TrackSegment>>> this_line_segments(
55✔
1384
    n_threads);
55✔
1385
  for (int t = 0; t < n_threads; ++t) {
140✔
1386
    this_line_segments[t].resize(pixels()[0]);
85✔
1387
  }
1388

1389
  // The last thread writes to this, and the first thread reads from it.
1390
  std::vector<std::vector<TrackSegment>> old_segments(pixels()[0]);
55✔
1391

1392
#pragma omp parallel
30✔
1393
  {
25✔
1394
    const int n_threads = num_threads();
25✔
1395
    const int tid = thread_num();
25✔
1396

1397
    int vert = tid;
25✔
1398
    for (int iter = 0; iter <= pixels()[1] / n_threads; iter++) {
5,050✔
1399

1400
      // Save bottom line of current work chunk to compare against later. This
1401
      // used to be inside the below if block, but it causes a spurious line to
1402
      // be drawn at the bottom of the image. Not sure why, but moving it here
1403
      // fixes things.
1404
      if (tid == n_threads - 1)
5,025✔
1405
        old_segments = this_line_segments[n_threads - 1];
5,025✔
1406

1407
      if (vert < pixels()[1]) {
5,025✔
1408

1409
        for (int horiz = 0; horiz < pixels()[0]; ++horiz) {
1,005,000✔
1410

1411
          // RayTracePlot implements camera ray generation
1412
          std::pair<Position, Direction> ru = get_pixel_ray(horiz, vert);
1,000,000✔
1413

1414
          this_line_segments[tid][horiz].clear();
1,000,000✔
1415
          ProjectionRay ray(
1,000,000✔
1416
            ru.first, ru.second, *this, this_line_segments[tid][horiz]);
1,000,000✔
1417

1418
          ray.trace();
1,000,000✔
1419

1420
          // Now color the pixel based on what we have intersected...
1421
          // Loops backwards over intersections.
1422
          Position current_color(
1,000,000✔
1423
            not_found_.red, not_found_.green, not_found_.blue);
1,000,000✔
1424
          const auto& segments = this_line_segments[tid][horiz];
1,000,000✔
1425

1426
          // There must be at least two cell intersections to color, front and
1427
          // back of the cell. Maybe an infinitely thick cell could be present
1428
          // with no back, but why would you want to color that? It's easier to
1429
          // just skip that edge case and not even color it.
1430
          if (segments.size() <= 1)
1,000,000✔
1431
            continue;
616,655✔
1432

1433
          for (int i = segments.size() - 2; i >= 0; --i) {
1,072,335✔
1434
            int colormap_idx = segments[i].id;
688,990✔
1435
            RGBColor seg_color = colors_[colormap_idx];
688,990✔
1436
            Position seg_color_vec(
688,990✔
1437
              seg_color.red, seg_color.green, seg_color.blue);
688,990✔
1438
            double mixing =
688,990✔
1439
              std::exp(-xs_[colormap_idx] *
1,377,980✔
1440
                       (segments[i + 1].length - segments[i].length));
688,990✔
1441
            current_color =
688,990✔
1442
              current_color * mixing + (1.0 - mixing) * seg_color_vec;
688,990✔
1443
          }
1444

1445
          // save result converting from double-precision color coordinates to
1446
          // byte-sized
1447
          RGBColor result;
383,345✔
1448
          result.red = static_cast<uint8_t>(current_color.x);
383,345✔
1449
          result.green = static_cast<uint8_t>(current_color.y);
383,345✔
1450
          result.blue = static_cast<uint8_t>(current_color.z);
383,345✔
1451
          data(horiz, vert) = result;
383,345✔
1452

1453
          // Check to draw wireframe in horizontal direction. No inter-thread
1454
          // comm.
1455
          if (horiz > 0) {
383,345✔
1456
            if (!trackstack_equivalent(this_line_segments[tid][horiz],
382,345✔
1457
                  this_line_segments[tid][horiz - 1])) {
382,345✔
1458
              wireframe_initial(horiz, vert) = 1;
15,710✔
1459
            }
1460
          }
1461
        }
1,000,000✔
1462
      } // end "if" vert in correct range
1463

1464
      // We require a barrier before comparing vertical neighbors' intersection
1465
      // stacks. i.e. all threads must be done with their line.
1466
#pragma omp barrier
1467

1468
      // Now that the horizontal line has finished rendering, we can fill in
1469
      // wireframe entries that require comparison among all the threads. Hence
1470
      // the omp barrier being used. It has to be OUTSIDE any if blocks!
1471
      if (vert < pixels()[1]) {
5,025✔
1472
        // Loop over horizontal pixels, checking intersection stack of upper
1473
        // neighbor
1474

1475
        const std::vector<std::vector<TrackSegment>>* top_cmp = nullptr;
1476
        if (tid == 0)
1477
          top_cmp = &old_segments;
1478
        else
1479
          top_cmp = &this_line_segments[tid - 1];
1480

1481
        for (int horiz = 0; horiz < pixels()[0]; ++horiz) {
1,005,000✔
1482
          if (!trackstack_equivalent(
1,000,000✔
1483
                this_line_segments[tid][horiz], (*top_cmp)[horiz])) {
1,000,000✔
1484
            wireframe_initial(horiz, vert) = 1;
20,595✔
1485
          }
1486
        }
1487
      }
1488

1489
      // We need another barrier to ensure threads don't proceed to modify their
1490
      // intersection stacks on that horizontal line while others are
1491
      // potentially still working on the above.
1492
#pragma omp barrier
1493
      vert += n_threads;
5,025✔
1494
    }
1495
  } // end omp parallel
1496

1497
  // Now thicken the wireframe lines and apply them to our image
1498
  for (int vert = 0; vert < pixels()[1]; ++vert) {
11,055✔
1499
    for (int horiz = 0; horiz < pixels()[0]; ++horiz) {
2,211,000✔
1500
      if (wireframe_initial(horiz, vert)) {
2,200,000✔
1501
        if (wireframe_thickness_ == 1)
70,983✔
1502
          data(horiz, vert) = wireframe_color_;
30,195✔
1503
        for (int i = -wireframe_thickness_ / 2; i < wireframe_thickness_ / 2;
195,723✔
1504
             ++i)
1505
          for (int j = -wireframe_thickness_ / 2; j < wireframe_thickness_ / 2;
546,876✔
1506
               ++j)
1507
            if (i * i + j * j < wireframe_thickness_ * wireframe_thickness_) {
422,136!
1508

1509
              // Check if wireframe pixel is out of bounds
1510
              int w_i = std::max(std::min(horiz + i, pixels()[0] - 1), 0);
422,136!
1511
              int w_j = std::max(std::min(vert + j, pixels()[1] - 1), 0);
422,268✔
1512
              data(w_i, w_j) = wireframe_color_;
422,136✔
1513
            }
1514
      }
1515
    }
1516
  }
1517

1518
  return data;
110✔
1519
}
110✔
1520

1521
void WireframeRayTracePlot::create_output() const
55✔
1522
{
1523
  ImageData data = create_image();
55✔
1524
  write_image(data);
55✔
1525
}
55✔
1526

1527
void RayTracePlot::print_info() const
88✔
1528
{
1529
  fmt::print("Camera position: {} {} {}\n", camera_position_.x,
176✔
1530
    camera_position_.y, camera_position_.z);
88✔
1531
  fmt::print("Look at: {} {} {}\n", look_at_.x, look_at_.y, look_at_.z);
88✔
1532
  fmt::print(
176✔
1533
    "Horizontal field of view: {} degrees\n", horizontal_field_of_view_);
88✔
1534
  fmt::print("Pixels: {} {}\n", pixels()[0], pixels()[1]);
88✔
1535
}
88✔
1536

1537
void WireframeRayTracePlot::print_info() const
55✔
1538
{
1539
  fmt::print("Plot Type: Wireframe ray-traced\n");
55✔
1540
  RayTracePlot::print_info();
55✔
1541
}
55✔
1542

1543
void WireframeRayTracePlot::set_opacities(pugi::xml_node node)
55✔
1544
{
1545
  xs_.resize(colors_.size(), 1e6); // set to large value for opaque by default
55✔
1546

1547
  for (auto cn : node.children("color")) {
121✔
1548
    // Make sure 3 values are specified for RGB
1549
    double user_xs = std::stod(get_node_value(cn, "xs"));
132✔
1550
    int col_id = std::stoi(get_node_value(cn, "id"));
132✔
1551

1552
    // Add RGB
1553
    if (PlotColorBy::cells == color_by_) {
66!
1554
      if (model::cell_map.find(col_id) != model::cell_map.end()) {
66!
1555
        col_id = model::cell_map[col_id];
66✔
1556
        xs_[col_id] = user_xs;
66✔
1557
      } else {
1558
        warning(fmt::format(
×
1559
          "Could not find cell {} specified in plot {}", col_id, id()));
×
1560
      }
1561
    } else if (PlotColorBy::mats == color_by_) {
×
1562
      if (model::material_map.find(col_id) != model::material_map.end()) {
×
1563
        col_id = model::material_map[col_id];
×
1564
        xs_[col_id] = user_xs;
×
1565
      } else {
1566
        warning(fmt::format(
×
1567
          "Could not find material {} specified in plot {}", col_id, id()));
×
1568
      }
1569
    }
1570
  }
1571
}
55✔
1572

1573
void RayTracePlot::set_orthographic_width(pugi::xml_node node)
88✔
1574
{
1575
  if (check_for_node(node, "orthographic_width")) {
88✔
1576
    double orthographic_width =
11✔
1577
      std::stod(get_node_value(node, "orthographic_width", true));
11✔
1578
    if (orthographic_width < 0.0)
11!
1579
      fatal_error("Requires positive orthographic_width");
×
1580
    orthographic_width_ = orthographic_width;
11✔
1581
  }
1582
}
88✔
1583

1584
void WireframeRayTracePlot::set_wireframe_thickness(pugi::xml_node node)
55✔
1585
{
1586
  if (check_for_node(node, "wireframe_thickness")) {
55✔
1587
    int wireframe_thickness =
22✔
1588
      std::stoi(get_node_value(node, "wireframe_thickness", true));
22✔
1589
    if (wireframe_thickness < 0)
22!
1590
      fatal_error("Requires non-negative wireframe thickness");
×
1591
    wireframe_thickness_ = wireframe_thickness;
22✔
1592
  }
1593
}
55✔
1594

1595
void WireframeRayTracePlot::set_wireframe_ids(pugi::xml_node node)
55✔
1596
{
1597
  if (check_for_node(node, "wireframe_ids")) {
55✔
1598
    wireframe_ids_ = get_node_array<int>(node, "wireframe_ids");
11✔
1599
    // It is read in as actual ID values, but we have to convert to indices in
1600
    // mat/cell array
1601
    for (auto& x : wireframe_ids_)
22✔
1602
      x = color_by_ == PlotColorBy::mats ? model::material_map[x]
22!
1603
                                         : model::cell_map[x];
×
1604
  }
1605
  // We make sure the list is sorted in order to later use
1606
  // std::binary_search.
1607
  std::sort(wireframe_ids_.begin(), wireframe_ids_.end());
55✔
1608
}
55✔
1609

1610
void RayTracePlot::set_pixels(pugi::xml_node node)
88✔
1611
{
1612
  vector<int> pxls = get_node_array<int>(node, "pixels");
88✔
1613
  if (pxls.size() != 2)
88!
1614
    fatal_error(
×
1615
      fmt::format("<pixels> must be length 2 in projection plot {}", id()));
×
1616
  pixels()[0] = pxls[0];
88✔
1617
  pixels()[1] = pxls[1];
88✔
1618
}
88✔
1619

1620
void RayTracePlot::set_camera_position(pugi::xml_node node)
88✔
1621
{
1622
  vector<double> camera_pos = get_node_array<double>(node, "camera_position");
88✔
1623
  if (camera_pos.size() != 3) {
88!
1624
    fatal_error(fmt::format(
×
1625
      "camera_position element must have three floating point values"));
1626
  }
1627
  camera_position_.x = camera_pos[0];
88✔
1628
  camera_position_.y = camera_pos[1];
88✔
1629
  camera_position_.z = camera_pos[2];
88✔
1630
}
88✔
1631

1632
void RayTracePlot::set_look_at(pugi::xml_node node)
88✔
1633
{
1634
  vector<double> look_at = get_node_array<double>(node, "look_at");
88✔
1635
  if (look_at.size() != 3) {
88!
1636
    fatal_error("look_at element must have three floating point values");
×
1637
  }
1638
  look_at_.x = look_at[0];
88✔
1639
  look_at_.y = look_at[1];
88✔
1640
  look_at_.z = look_at[2];
88✔
1641
}
88✔
1642

1643
void RayTracePlot::set_field_of_view(pugi::xml_node node)
88✔
1644
{
1645
  // Defaults to 70 degree horizontal field of view (see .h file)
1646
  if (check_for_node(node, "horizontal_field_of_view")) {
88!
1647
    double fov =
×
1648
      std::stod(get_node_value(node, "horizontal_field_of_view", true));
×
1649
    if (fov < 180.0 && fov > 0.0) {
×
1650
      horizontal_field_of_view_ = fov;
×
1651
    } else {
1652
      fatal_error(fmt::format("Horizontal field of view for plot {} "
×
1653
                              "out-of-range. Must be in (0, 180) degrees.",
1654
        id()));
×
1655
    }
1656
  }
1657
}
88✔
1658

1659
SolidRayTracePlot::SolidRayTracePlot(pugi::xml_node node) : RayTracePlot(node)
33✔
1660
{
1661
  set_opaque_ids(node);
33✔
1662
  set_diffuse_fraction(node);
33✔
1663
  set_light_position(node);
33✔
1664
  update_view();
33✔
1665
}
33✔
1666

1667
void SolidRayTracePlot::print_info() const
33✔
1668
{
1669
  fmt::print("Plot Type: Solid ray-traced\n");
33✔
1670
  RayTracePlot::print_info();
33✔
1671
}
33✔
1672

1673
ImageData SolidRayTracePlot::create_image() const
55✔
1674
{
1675
  size_t width = pixels()[0];
55✔
1676
  size_t height = pixels()[1];
55✔
1677
  ImageData data({width, height}, not_found_);
55✔
1678

1679
#pragma omp parallel for schedule(dynamic) collapse(2)
30✔
1680
  for (int horiz = 0; horiz < pixels()[0]; ++horiz) {
3,105✔
1681
    for (int vert = 0; vert < pixels()[1]; ++vert) {
603,560✔
1682
      // RayTracePlot implements camera ray generation
1683
      std::pair<Position, Direction> ru = get_pixel_ray(horiz, vert);
600,480✔
1684
      PhongRay ray(ru.first, ru.second, *this);
600,480✔
1685
      ray.trace();
600,480✔
1686
      data(horiz, vert) = ray.result_color();
600,480✔
1687
    }
600,480✔
1688
  }
1689

1690
  return data;
55✔
1691
}
1692

1693
void SolidRayTracePlot::create_output() const
33✔
1694
{
1695
  ImageData data = create_image();
33✔
1696
  write_image(data);
33✔
1697
}
33✔
1698

1699
void SolidRayTracePlot::set_opaque_ids(pugi::xml_node node)
33✔
1700
{
1701
  if (check_for_node(node, "opaque_ids")) {
33!
1702
    auto opaque_ids_tmp = get_node_array<int>(node, "opaque_ids");
33✔
1703

1704
    // It is read in as actual ID values, but we have to convert to indices in
1705
    // mat/cell array
1706
    for (auto& x : opaque_ids_tmp)
99✔
1707
      x = color_by_ == PlotColorBy::mats ? model::material_map[x]
132!
1708
                                         : model::cell_map[x];
×
1709

1710
    opaque_ids_.insert(opaque_ids_tmp.begin(), opaque_ids_tmp.end());
33✔
1711
  }
33✔
1712
}
33✔
1713

1714
void SolidRayTracePlot::set_light_position(pugi::xml_node node)
33✔
1715
{
1716
  if (check_for_node(node, "light_position")) {
33✔
1717
    auto light_pos_tmp = get_node_array<double>(node, "light_position");
11✔
1718

1719
    if (light_pos_tmp.size() != 3)
11!
1720
      fatal_error("Light position must be given as 3D coordinates");
×
1721

1722
    light_location_.x = light_pos_tmp[0];
11✔
1723
    light_location_.y = light_pos_tmp[1];
11✔
1724
    light_location_.z = light_pos_tmp[2];
11✔
1725
  } else {
11✔
1726
    light_location_ = camera_position();
22✔
1727
  }
1728
}
33✔
1729

1730
void SolidRayTracePlot::set_diffuse_fraction(pugi::xml_node node)
33✔
1731
{
1732
  if (check_for_node(node, "diffuse_fraction")) {
33✔
1733
    diffuse_fraction_ = std::stod(get_node_value(node, "diffuse_fraction"));
11✔
1734
    if (diffuse_fraction_ < 0.0 || diffuse_fraction_ > 1.0) {
11!
1735
      fatal_error("Must have 0 <= diffuse fraction <= 1");
×
1736
    }
1737
  }
1738
}
33✔
1739

1740
void Ray::compute_distance()
3,014,638✔
1741
{
1742
  boundary() = distance_to_boundary(*this);
3,014,638✔
1743
}
3,014,638✔
1744

1745
void Ray::trace()
3,521,056✔
1746
{
1747
  // To trace the ray from its origin all the way through the model, we have
1748
  // to proceed in two phases. In the first, the ray may or may not be found
1749
  // inside the model. If the ray is already in the model, phase one can be
1750
  // skipped. Otherwise, the ray has to be advanced to the boundary of the
1751
  // model where all the cells are defined. Importantly, this is assuming that
1752
  // the model is convex, which is a very reasonable assumption for any
1753
  // radiation transport model.
1754
  //
1755
  // After phase one is done, we can starting tracing from cell to cell within
1756
  // the model. This step can use neighbor lists to accelerate the ray tracing.
1757

1758
  bool inside_cell;
3,521,056✔
1759
  // Check for location if the particle is already known
1760
  if (lowest_coord().cell() == C_NONE) {
3,521,056!
1761
    // The geometry position of the particle is either unknown or outside of the
1762
    // edge of the model.
1763
    if (lowest_coord().universe() == C_NONE) {
3,521,056!
1764
      // Attempt to initialize the particle. We may have to
1765
      // enter a loop to move it up to the edge of the model.
1766
      inside_cell = exhaustive_find_cell(*this, settings::verbosity >= 10);
3,521,056✔
1767
    } else {
1768
      // It has been already calculated that the current position is outside of
1769
      // the edge of the model.
1770
      inside_cell = false;
1771
    }
1772
  } else {
1773
    // Availability of the cell means that the particle is located inside the
1774
    // edge.
1775
    inside_cell = true;
1776
  }
1777

1778
  // Advance to the boundary of the model
1779
  while (!inside_cell) {
15,618,438!
1780
    advance_to_boundary_from_void();
15,618,438✔
1781
    inside_cell = exhaustive_find_cell(*this, settings::verbosity >= 10);
15,618,438✔
1782

1783
    // If true this means no surface was intersected. See cell.cpp and search
1784
    // for numeric_limits to see where we return it.
1785
    if (surface() == std::numeric_limits<int>::max()) {
15,618,438!
1786
      warning(fmt::format("Lost a ray, r = {}, u = {}", r(), u()));
×
1787
      return;
×
1788
    }
1789

1790
    // Exit this loop and enter into cell-to-cell ray tracing (which uses
1791
    // neighbor lists)
1792
    if (inside_cell)
15,618,438✔
1793
      break;
1794

1795
    // if there is no intersection with the model, we're done
1796
    if (boundary().surface() == SURFACE_NONE)
14,068,604✔
1797
      return;
1798

1799
    event_counter_++;
12,097,382✔
1800
    if (event_counter_ > MAX_INTERSECTIONS) {
12,097,382!
1801
      warning("Likely infinite loop in ray traced plot");
×
1802
      return;
×
1803
    }
1804
  }
1805

1806
  // Call the specialized logic for this type of ray. This is for the
1807
  // intersection for the first intersection if we had one.
1808
  if (boundary().surface() != SURFACE_NONE) {
1,549,834!
1809
    // set the geometry state's surface attribute to be used for
1810
    // surface normal computation
1811
    surface() = boundary().surface();
1,549,834✔
1812
    on_intersection();
1,549,834✔
1813
    if (stop_)
1,549,834!
1814
      return;
1815
  }
1816

1817
  // reset surface attribute to zero after the first intersection so that it
1818
  // doesn't perturb surface crossing logic from here on out
1819
  surface() = 0;
1,549,834✔
1820

1821
  // This is the ray tracing loop within the model. It exits after exiting
1822
  // the model, which is equivalent to assuming that the model is convex.
1823
  // It would be nice to factor out the on_intersection at the end of this
1824
  // loop and then do "while (inside_cell)", but we can't guarantee it's
1825
  // on a surface in that case. There might be some other way to set it
1826
  // up that is perhaps a little more elegant, but this is what works just
1827
  // fine.
1828
  while (true) {
2,308,570✔
1829

1830
    compute_distance();
2,308,570✔
1831

1832
    // There are no more intersections to process
1833
    // if we hit the edge of the model, so stop
1834
    // the particle in that case. Also, just exit
1835
    // if a negative distance was somehow computed.
1836
    if (boundary().distance() == INFTY || boundary().distance() == INFINITY ||
2,308,570!
1837
        boundary().distance() < 0) {
2,308,570!
1838
      return;
1839
    }
1840

1841
    // See below comment where call_on_intersection is checked in an
1842
    // if statement for an explanation of this.
1843
    bool call_on_intersection {true};
2,308,570✔
1844
    if (boundary().distance() < 10 * TINY_BIT) {
2,308,570✔
1845
      call_on_intersection = false;
593,285✔
1846
    }
1847

1848
    // DAGMC surfaces expect us to go a little bit further than the advance
1849
    // distance to properly check cell inclusion.
1850
    boundary().distance() += TINY_BIT;
2,308,570✔
1851

1852
    // Advance particle, prepare for next intersection
1853
    for (int lev = 0; lev < n_coord(); ++lev) {
4,617,140✔
1854
      coord(lev).r() += boundary().distance() * coord(lev).u();
2,308,570✔
1855
    }
1856
    surface() = boundary().surface();
2,308,570✔
1857
    // Initialize last cells from the current cell, because the cell() variable
1858
    // does not contain the data for the case of a single-segment ray
1859
    for (int j = 0; j < n_coord(); ++j) {
4,617,140✔
1860
      cell_last(j) = coord(j).cell();
2,308,570✔
1861
    }
1862
    n_coord_last() = n_coord();
2,308,570✔
1863
    n_coord() = boundary().coord_level();
2,308,570!
1864
    if (boundary().lattice_translation()[0] != 0 ||
2,308,570!
1865
        boundary().lattice_translation()[1] != 0 ||
2,308,570!
1866
        boundary().lattice_translation()[2] != 0) {
2,308,570!
1867
      cross_lattice(*this, boundary(), settings::verbosity >= 10);
×
1868
    }
1869

1870
    // Record how far the ray has traveled
1871
    traversal_distance_ += boundary().distance();
2,308,570✔
1872
    inside_cell = neighbor_list_find_cell(*this, settings::verbosity >= 10);
2,308,570✔
1873

1874
    // Call the specialized logic for this type of ray. Note that we do not
1875
    // call this if the advance distance is very small. Unfortunately, it seems
1876
    // darn near impossible to get the particle advanced to the model boundary
1877
    // and through it without sometimes accidentally calling on_intersection
1878
    // twice. This incorrectly shades the region as occluded when it might not
1879
    // actually be. By screening out intersection distances smaller than a
1880
    // threshold 10x larger than the scoot distance used to advance up to the
1881
    // model boundary, we can avoid that situation.
1882
    if (call_on_intersection) {
2,308,570✔
1883
      on_intersection();
1,715,285✔
1884
      if (stop_)
1,715,285✔
1885
        return;
1886
    }
1887

1888
    if (!inside_cell)
2,273,007✔
1889
      return;
1890

1891
    event_counter_++;
758,736✔
1892
    if (event_counter_ > MAX_INTERSECTIONS) {
758,736!
1893
      warning("Likely infinite loop in ray traced plot");
×
1894
      return;
×
1895
    }
1896
  }
1897
}
1898

1899
void ProjectionRay::on_intersection()
2,359,148✔
1900
{
1901
  // This records a tuple with the following info
1902
  //
1903
  // 1) ID (material or cell depending on color_by_)
1904
  // 2) Distance traveled by the ray through that ID
1905
  // 3) Index of the intersected surface (starting from 1)
1906

1907
  line_segments_.emplace_back(
2,359,148✔
1908
    plot_.color_by_ == PlottableInterface::PlotColorBy::mats
2,359,148✔
1909
      ? material()
545,919✔
1910
      : lowest_coord().cell(),
1,813,229✔
1911
    traversal_distance_, boundary().surface_index());
2,359,148✔
1912
}
2,359,148✔
1913

1914
void PhongRay::on_intersection()
905,971✔
1915
{
1916
  // Check if we hit an opaque material or cell
1917
  int hit_id = plot_.color_by_ == PlottableInterface::PlotColorBy::mats
905,971✔
1918
                 ? material()
905,971!
1919
                 : lowest_coord().cell();
×
1920

1921
  // If we are reflected and have advanced beyond the camera,
1922
  // the ray is done. This is checked here because we should
1923
  // kill the ray even if the material is not opaque.
1924
  if (reflected_ && (r() - plot_.camera_position()).dot(u()) >= 0.0) {
905,971!
1925
    stop();
×
1926
    return;
164,340✔
1927
  }
1928

1929
  // Anything that's not opaque has zero impact on the plot.
1930
  if (plot_.opaque_ids_.find(hit_id) == plot_.opaque_ids_.end())
905,971✔
1931
    return;
1932

1933
  if (!reflected_) {
741,631✔
1934
    // reflect the particle and set the color to be colored by
1935
    // the normal or the diffuse lighting contribution
1936
    reflected_ = true;
706,068✔
1937
    result_color_ = plot_.colors_[hit_id];
706,068✔
1938
    Direction to_light = plot_.light_location_ - r();
706,068✔
1939
    to_light /= to_light.norm();
706,068✔
1940

1941
    // TODO
1942
    // Not sure what can cause a surface token to be invalid here, although it
1943
    // sometimes happens for a few pixels. It's very very rare, so proceed by
1944
    // coloring the pixel with the overlap color. It seems to happen only for a
1945
    // few pixels on the outer boundary of a hex lattice.
1946
    //
1947
    // We cannot detect it in the outer loop, and it only matters here, so
1948
    // that's why the error handling is a little different than for a lost
1949
    // ray.
1950
    if (surface() == 0) {
706,068!
1951
      result_color_ = plot_.overlap_color_;
×
1952
      stop();
×
1953
      return;
×
1954
    }
1955

1956
    // Get surface pointer
1957
    const auto& surf = model::surfaces.at(surface_index());
706,068✔
1958

1959
    Direction normal = surf->normal(r_local());
706,068✔
1960
    normal /= normal.norm();
706,068✔
1961

1962
    // Need to apply translations to find the normal vector in
1963
    // the base level universe's coordinate system.
1964
    for (int lev = n_coord() - 2; lev >= 0; --lev) {
706,068!
1965
      if (coord(lev + 1).rotated()) {
×
1966
        const Cell& c {*model::cells[coord(lev).cell()]};
×
1967
        normal = normal.inverse_rotate(c.rotation_);
×
1968
      }
1969
    }
1970

1971
    // use the normal opposed to the ray direction
1972
    if (normal.dot(u()) > 0.0) {
706,068✔
1973
      normal *= -1.0;
63,789✔
1974
    }
1975

1976
    // Facing away from the light means no lighting
1977
    double dotprod = normal.dot(to_light);
706,068✔
1978
    dotprod = std::max(0.0, dotprod);
706,068✔
1979

1980
    double modulation =
706,068✔
1981
      plot_.diffuse_fraction_ + (1.0 - plot_.diffuse_fraction_) * dotprod;
706,068✔
1982
    result_color_ *= modulation;
706,068✔
1983

1984
    // Now point the particle to the camera. We now begin
1985
    // checking to see if it's occluded by another surface
1986
    u() = to_light;
706,068✔
1987

1988
    orig_hit_id_ = hit_id;
706,068✔
1989

1990
    // OpenMC native CSG and DAGMC surfaces have some slight differences
1991
    // in how they interpret particles that are sitting on a surface.
1992
    // I don't know exactly why, but this makes everything work beautifully.
1993
    if (surf->geom_type() == GeometryType::DAG) {
706,068!
1994
      surface() = 0;
×
1995
    } else {
1996
      surface() = -surface(); // go to other side
706,068✔
1997
    }
1998

1999
    // Must fully restart coordinate search. Why? Not sure.
2000
    clear();
706,068✔
2001

2002
    // Note this could likely be faster if we cached the previous
2003
    // cell we were in before the reflection. This is the easiest
2004
    // way to fully initialize all the sub-universe coordinates and
2005
    // directions though.
2006
    bool found = exhaustive_find_cell(*this);
706,068✔
2007
    if (!found) {
706,068!
2008
      fatal_error("Lost particle after reflection.");
×
2009
    }
2010

2011
    // Must recalculate distance to boundary due to the
2012
    // direction change
2013
    compute_distance();
706,068✔
2014

2015
  } else {
2016
    // If it's not facing the light, we color with the diffuse contribution, so
2017
    // next we check if we're going to occlude the last reflected surface. if
2018
    // so, color by the diffuse contribution instead
2019

2020
    if (orig_hit_id_ == -1)
35,563!
2021
      fatal_error("somehow a ray got reflected but not original ID set?");
×
2022

2023
    result_color_ = plot_.colors_[orig_hit_id_];
35,563✔
2024
    result_color_ *= plot_.diffuse_fraction_;
35,563✔
2025
    stop();
741,631✔
2026
  }
2027
}
2028

2029
extern "C" int openmc_id_map(const void* plot, int32_t* data_out)
278✔
2030
{
2031

2032
  auto plt = reinterpret_cast<const SlicePlotBase*>(plot);
278✔
2033
  if (!plt) {
278!
2034
    set_errmsg("Invalid slice pointer passed to openmc_id_map");
×
2035
    return OPENMC_E_INVALID_ARGUMENT;
×
2036
  }
2037

2038
  if (plt->slice_color_overlaps_ && model::overlap_check_count.size() == 0) {
278!
2039
    model::overlap_check_count.resize(model::cells.size());
22✔
2040
  }
2041

2042
  auto ids = plt->get_map<IdData>();
278✔
2043

2044
  // write id data to array
2045
  std::copy(ids.data_.begin(), ids.data_.end(), data_out);
278✔
2046

2047
  return 0;
278✔
2048
}
278✔
2049

2050
extern "C" int openmc_property_map(const void* plot, double* data_out)
11✔
2051
{
2052

2053
  auto plt = reinterpret_cast<const SlicePlotBase*>(plot);
11✔
2054
  if (!plt) {
11!
NEW
2055
    set_errmsg("Invalid slice pointer passed to openmc_property_map");
×
2056
    return OPENMC_E_INVALID_ARGUMENT;
×
2057
  }
2058

2059
  if (plt->slice_color_overlaps_ && model::overlap_check_count.size() == 0) {
11!
2060
    model::overlap_check_count.resize(model::cells.size());
×
2061
  }
2062

2063
  auto props = plt->get_map<PropertyData>();
11✔
2064

2065
  // write id data to array
2066
  std::copy(props.data_.begin(), props.data_.end(), data_out);
11✔
2067

2068
  return 0;
11✔
2069
}
11✔
2070

2071
extern "C" int openmc_slice_plot(const double origin[3], const double u_span[3],
110✔
2072
  const double v_span[3], const size_t pixels[2], bool color_overlaps,
2073
  int level, int32_t filter_index, int32_t* geom_data, double* property_data)
2074
{
2075
  // Validate span vectors
2076
  Position u_span_pos {u_span[0], u_span[1], u_span[2]};
110✔
2077
  Position v_span_pos {v_span[0], v_span[1], v_span[2]};
110✔
2078
  double u_norm = u_span_pos.norm();
110✔
2079
  double v_norm = v_span_pos.norm();
110✔
2080
  if (u_norm == 0.0 || v_norm == 0.0) {
110!
NEW
2081
    set_errmsg("Slice span vectors must be non-zero.");
×
NEW
2082
    return OPENMC_E_INVALID_ARGUMENT;
×
2083
  }
2084

2085
  constexpr double ORTHO_REL_TOL = 1e-10;
110✔
2086
  double dot = u_span_pos.dot(v_span_pos);
110!
2087
  if (std::abs(dot) > ORTHO_REL_TOL * u_norm * v_norm) {
110!
NEW
2088
    set_errmsg("Slice span vectors must be orthogonal.");
×
NEW
2089
    return OPENMC_E_INVALID_ARGUMENT;
×
2090
  }
2091

2092
  // Validate filter index if provided
2093
  if (filter_index >= 0) {
110✔
2094
    if (int err = verify_filter(filter_index))
22!
2095
      return err;
2096
  }
2097

2098
  // Initialize overlap check vector if needed
2099
  if (color_overlaps && model::overlap_check_count.size() == 0) {
110!
2100
    model::overlap_check_count.resize(model::cells.size());
22✔
2101
  }
2102

2103
  try {
110✔
2104
    // Create a temporary SlicePlotBase object to reuse get_map logic
2105
    SlicePlotBase plot_params;
110✔
2106
    plot_params.origin_ = Position {origin[0], origin[1], origin[2]};
110✔
2107
    plot_params.u_span_ = u_span_pos;
110✔
2108
    plot_params.v_span_ = v_span_pos;
110✔
2109
    plot_params.width_ = Position {u_norm, v_norm, 0.0};
110✔
2110
    plot_params.basis_ = SlicePlotBase::PlotBasis::xy;
110✔
2111
    plot_params.pixels_[0] = pixels[0];
110✔
2112
    plot_params.pixels_[1] = pixels[1];
110✔
2113
    plot_params.slice_color_overlaps_ = color_overlaps;
110✔
2114
    plot_params.slice_level_ = level;
110✔
2115

2116
    // Use get_map<RasterData> to generate data
2117
    auto data = plot_params.get_map<RasterData>(filter_index);
110✔
2118

2119
    // Copy geometry data
2120
    std::copy(data.id_data_.begin(), data.id_data_.end(), geom_data);
110✔
2121

2122
    // Copy property data if requested
2123
    if (property_data != nullptr) {
110✔
2124
      std::copy(
55✔
2125
        data.property_data_.begin(), data.property_data_.end(), property_data);
2126
    }
2127
  } catch (const std::exception& e) {
110!
NEW
2128
    set_errmsg(e.what());
×
NEW
2129
    return OPENMC_E_UNASSIGNED;
×
NEW
2130
  }
×
2131

2132
  return 0;
110✔
2133
}
2134

2135
extern "C" int openmc_get_plot_index(int32_t id, int32_t* index)
22✔
2136
{
2137
  auto it = model::plot_map.find(id);
22!
2138
  if (it == model::plot_map.end()) {
22!
2139
    set_errmsg("No plot exists with ID=" + std::to_string(id) + ".");
×
2140
    return OPENMC_E_INVALID_ID;
×
2141
  }
2142

2143
  *index = it->second;
22✔
2144
  return 0;
22✔
2145
}
2146

2147
extern "C" int openmc_plot_get_id(int32_t index, int32_t* id)
55✔
2148
{
2149
  if (index < 0 || index >= model::plots.size()) {
55!
2150
    set_errmsg("Index in plots array is out of bounds.");
×
2151
    return OPENMC_E_OUT_OF_BOUNDS;
×
2152
  }
2153

2154
  *id = model::plots[index]->id();
55✔
2155
  return 0;
55✔
2156
}
2157

2158
extern "C" int openmc_plot_set_id(int32_t index, int32_t id)
×
2159
{
2160
  if (index < 0 || index >= model::plots.size()) {
×
2161
    set_errmsg("Index in plots array is out of bounds.");
×
2162
    return OPENMC_E_OUT_OF_BOUNDS;
×
2163
  }
2164

2165
  if (id < 0 && id != C_NONE) {
×
2166
    set_errmsg("Invalid plot ID.");
×
2167
    return OPENMC_E_INVALID_ARGUMENT;
×
2168
  }
2169

2170
  auto* plot = model::plots[index].get();
×
2171
  int32_t old_id = plot->id();
×
2172
  if (id == old_id)
×
2173
    return 0;
2174

2175
  model::plot_map.erase(old_id);
×
2176
  try {
×
2177
    plot->set_id(id);
×
2178
  } catch (const std::runtime_error& e) {
×
2179
    model::plot_map[old_id] = index;
×
2180
    set_errmsg(e.what());
×
2181
    return OPENMC_E_INVALID_ID;
×
2182
  }
×
2183
  model::plot_map[plot->id()] = index;
×
2184
  return 0;
×
2185
}
2186

2187
extern "C" size_t openmc_plots_size()
22✔
2188
{
2189
  return model::plots.size();
22✔
2190
}
2191

2192
int map_phong_domain_id(
55✔
2193
  const SolidRayTracePlot* plot, int32_t id, int32_t* index_out)
2194
{
2195
  if (!plot || !index_out) {
55!
2196
    set_errmsg("Invalid plot pointer passed to map_phong_domain_id");
×
2197
    return OPENMC_E_INVALID_ARGUMENT;
×
2198
  }
2199

2200
  if (plot->color_by_ == PlottableInterface::PlotColorBy::mats) {
55!
2201
    auto it = model::material_map.find(id);
55!
2202
    if (it == model::material_map.end()) {
55!
2203
      set_errmsg("Invalid material ID for SolidRayTracePlot");
×
2204
      return OPENMC_E_INVALID_ID;
×
2205
    }
2206
    *index_out = it->second;
55✔
2207
    return 0;
55✔
2208
  }
2209

2210
  if (plot->color_by_ == PlottableInterface::PlotColorBy::cells) {
×
2211
    auto it = model::cell_map.find(id);
×
2212
    if (it == model::cell_map.end()) {
×
2213
      set_errmsg("Invalid cell ID for SolidRayTracePlot");
×
2214
      return OPENMC_E_INVALID_ID;
×
2215
    }
2216
    *index_out = it->second;
×
2217
    return 0;
×
2218
  }
2219

2220
  set_errmsg("Unsupported color_by for SolidRayTracePlot");
×
2221
  return OPENMC_E_INVALID_TYPE;
×
2222
}
2223

2224
int get_solidraytrace_plot_by_index(int32_t index, SolidRayTracePlot** plot)
308✔
2225
{
2226
  if (!plot) {
308!
2227
    set_errmsg("Null output pointer passed to get_solidraytrace_plot_by_index");
×
2228
    return OPENMC_E_INVALID_ARGUMENT;
×
2229
  }
2230

2231
  if (index < 0 || index >= model::plots.size()) {
308!
2232
    set_errmsg("Index in plots array is out of bounds.");
×
2233
    return OPENMC_E_OUT_OF_BOUNDS;
×
2234
  }
2235

2236
  auto* plottable = model::plots[index].get();
308!
2237
  auto* solid_plot = dynamic_cast<SolidRayTracePlot*>(plottable);
308!
2238
  if (!solid_plot) {
308!
2239
    set_errmsg("Plot at index=" + std::to_string(index) +
×
2240
               " is not a solid raytrace plot.");
2241
    return OPENMC_E_INVALID_TYPE;
×
2242
  }
2243

2244
  *plot = solid_plot;
308✔
2245
  return 0;
308✔
2246
}
2247

2248
extern "C" int openmc_solidraytrace_plot_create(int32_t* index)
11✔
2249
{
2250
  if (!index) {
11!
2251
    set_errmsg(
×
2252
      "Null output pointer passed to openmc_solidraytrace_plot_create");
2253
    return OPENMC_E_INVALID_ARGUMENT;
×
2254
  }
2255

2256
  try {
11✔
2257
    auto new_plot = std::make_unique<SolidRayTracePlot>();
11✔
2258
    new_plot->set_id();
11✔
2259
    int32_t new_plot_id = new_plot->id();
11✔
2260
#ifdef USE_LIBPNG
2261
    new_plot->path_plot() = fmt::format("plot_{}.png", new_plot_id);
11✔
2262
#else
2263
    new_plot->path_plot() = fmt::format("plot_{}.ppm", new_plot_id);
2264
#endif
2265
    int32_t new_plot_index = model::plots.size();
11✔
2266
    model::plots.emplace_back(std::move(new_plot));
11✔
2267
    model::plot_map[new_plot_id] = new_plot_index;
11✔
2268
    *index = new_plot_index;
11✔
2269
  } catch (const std::exception& e) {
11!
2270
    set_errmsg(e.what());
×
2271
    return OPENMC_E_ALLOCATE;
×
2272
  }
×
2273

2274
  return 0;
11✔
2275
}
2276

2277
extern "C" int openmc_solidraytrace_plot_get_pixels(
33✔
2278
  int32_t index, int32_t* width, int32_t* height)
2279
{
2280
  if (!width || !height) {
33!
2281
    set_errmsg(
×
2282
      "Invalid arguments passed to openmc_solidraytrace_plot_get_pixels");
2283
    return OPENMC_E_INVALID_ARGUMENT;
×
2284
  }
2285

2286
  SolidRayTracePlot* plt = nullptr;
33✔
2287
  int err = get_solidraytrace_plot_by_index(index, &plt);
33✔
2288
  if (err)
33!
2289
    return err;
2290

2291
  *width = plt->pixels()[0];
33✔
2292
  *height = plt->pixels()[1];
33✔
2293
  return 0;
33✔
2294
}
2295

2296
extern "C" int openmc_solidraytrace_plot_set_pixels(
11✔
2297
  int32_t index, int32_t width, int32_t height)
2298
{
2299
  if (width <= 0 || height <= 0) {
11!
2300
    set_errmsg(
×
2301
      "Invalid arguments passed to openmc_solidraytrace_plot_set_pixels");
2302
    return OPENMC_E_INVALID_ARGUMENT;
×
2303
  }
2304

2305
  SolidRayTracePlot* plt = nullptr;
11✔
2306
  int err = get_solidraytrace_plot_by_index(index, &plt);
11✔
2307
  if (err)
11!
2308
    return err;
2309

2310
  plt->pixels()[0] = width;
11✔
2311
  plt->pixels()[1] = height;
11✔
2312
  return 0;
11✔
2313
}
2314

2315
extern "C" int openmc_solidraytrace_plot_get_color_by(
11✔
2316
  int32_t index, int32_t* color_by)
2317
{
2318
  if (!color_by) {
11!
2319
    set_errmsg(
×
2320
      "Invalid arguments passed to openmc_solidraytrace_plot_get_color_by");
2321
    return OPENMC_E_INVALID_ARGUMENT;
×
2322
  }
2323

2324
  SolidRayTracePlot* plt = nullptr;
11✔
2325
  int err = get_solidraytrace_plot_by_index(index, &plt);
11✔
2326
  if (err)
11!
2327
    return err;
2328

2329
  if (plt->color_by_ == PlottableInterface::PlotColorBy::mats) {
11!
2330
    *color_by = 0;
11✔
2331
  } else if (plt->color_by_ == PlottableInterface::PlotColorBy::cells) {
×
2332
    *color_by = 1;
×
2333
  } else {
2334
    set_errmsg("Unsupported color_by for SolidRayTracePlot");
×
2335
    return OPENMC_E_INVALID_TYPE;
×
2336
  }
2337

2338
  return 0;
2339
}
2340

2341
extern "C" int openmc_solidraytrace_plot_set_color_by(
11✔
2342
  int32_t index, int32_t color_by)
2343
{
2344
  SolidRayTracePlot* plt = nullptr;
11✔
2345
  int err = get_solidraytrace_plot_by_index(index, &plt);
11✔
2346
  if (err)
11!
2347
    return err;
2348

2349
  if (color_by == 0) {
11!
2350
    plt->color_by_ = PlottableInterface::PlotColorBy::mats;
11✔
2351
  } else if (color_by == 1) {
×
2352
    plt->color_by_ = PlottableInterface::PlotColorBy::cells;
×
2353
  } else {
2354
    set_errmsg("Invalid color_by value for SolidRayTracePlot");
×
2355
    return OPENMC_E_INVALID_ARGUMENT;
×
2356
  }
2357

2358
  return 0;
2359
}
2360

2361
extern "C" int openmc_solidraytrace_plot_set_default_colors(int32_t index)
11✔
2362
{
2363
  SolidRayTracePlot* plt = nullptr;
11✔
2364
  int err = get_solidraytrace_plot_by_index(index, &plt);
11✔
2365
  if (err)
11!
2366
    return err;
2367

2368
  plt->set_default_colors();
11✔
2369
  return 0;
2370
}
2371

2372
extern "C" int openmc_solidraytrace_plot_set_all_opaque(int32_t index)
×
2373
{
2374
  SolidRayTracePlot* plt = nullptr;
×
2375
  int err = get_solidraytrace_plot_by_index(index, &plt);
×
2376
  if (err)
×
2377
    return err;
2378

2379
  plt->opaque_ids().clear();
×
2380
  if (plt->color_by_ == PlottableInterface::PlotColorBy::mats) {
×
2381
    for (int32_t i = 0; i < model::materials.size(); ++i) {
×
2382
      plt->opaque_ids().insert(i);
×
2383
    }
2384
    return 0;
×
2385
  }
2386

2387
  if (plt->color_by_ == PlottableInterface::PlotColorBy::cells) {
×
2388
    for (int32_t i = 0; i < model::cells.size(); ++i) {
×
2389
      plt->opaque_ids().insert(i);
×
2390
    }
2391
    return 0;
×
2392
  }
2393

2394
  set_errmsg("Unsupported color_by for SolidRayTracePlot");
×
2395
  return OPENMC_E_INVALID_TYPE;
×
2396
}
2397

2398
extern "C" int openmc_solidraytrace_plot_set_opaque(
22✔
2399
  int32_t index, int32_t id, bool visible)
2400
{
2401
  SolidRayTracePlot* plt = nullptr;
22✔
2402
  int err = get_solidraytrace_plot_by_index(index, &plt);
22✔
2403
  if (err)
22!
2404
    return err;
2405

2406
  int32_t domain_index = -1;
22✔
2407
  err = map_phong_domain_id(plt, id, &domain_index);
22✔
2408
  if (err)
22!
2409
    return err;
2410

2411
  if (visible) {
22✔
2412
    plt->opaque_ids().insert(domain_index);
11✔
2413
  } else {
2414
    plt->opaque_ids().erase(domain_index);
11✔
2415
  }
2416

2417
  return 0;
2418
}
2419

2420
extern "C" int openmc_solidraytrace_plot_set_color(
22✔
2421
  int32_t index, int32_t id, uint8_t r, uint8_t g, uint8_t b)
2422
{
2423
  SolidRayTracePlot* plt = nullptr;
22✔
2424
  int err = get_solidraytrace_plot_by_index(index, &plt);
22✔
2425
  if (err)
22!
2426
    return err;
2427

2428
  int32_t domain_index = -1;
22✔
2429
  err = map_phong_domain_id(plt, id, &domain_index);
22✔
2430
  if (err)
22!
2431
    return err;
2432

2433
  if (domain_index < 0 ||
22!
2434
      static_cast<size_t>(domain_index) >= plt->colors_.size()) {
22!
2435
    set_errmsg("Color index out of range for SolidRayTracePlot");
×
2436
    return OPENMC_E_OUT_OF_BOUNDS;
×
2437
  }
2438

2439
  plt->colors_[domain_index] = RGBColor(r, g, b);
22✔
2440
  return 0;
22✔
2441
}
2442

2443
extern "C" int openmc_solidraytrace_plot_get_camera_position(
11✔
2444
  int32_t index, double* x, double* y, double* z)
2445
{
2446
  if (!x || !y || !z) {
11!
2447
    set_errmsg("Invalid arguments passed to "
×
2448
               "openmc_solidraytrace_plot_get_camera_position");
2449
    return OPENMC_E_INVALID_ARGUMENT;
×
2450
  }
2451

2452
  SolidRayTracePlot* plt = nullptr;
11✔
2453
  int err = get_solidraytrace_plot_by_index(index, &plt);
11✔
2454
  if (err)
11!
2455
    return err;
2456

2457
  const auto& camera_position = plt->camera_position();
11✔
2458
  *x = camera_position.x;
11✔
2459
  *y = camera_position.y;
11✔
2460
  *z = camera_position.z;
11✔
2461
  return 0;
11✔
2462
}
2463

2464
extern "C" int openmc_solidraytrace_plot_set_camera_position(
11✔
2465
  int32_t index, double x, double y, double z)
2466
{
2467
  SolidRayTracePlot* plt = nullptr;
11✔
2468
  int err = get_solidraytrace_plot_by_index(index, &plt);
11✔
2469
  if (err)
11!
2470
    return err;
2471

2472
  plt->camera_position() = {x, y, z};
11✔
2473
  return 0;
11✔
2474
}
2475

2476
extern "C" int openmc_solidraytrace_plot_get_look_at(
11✔
2477
  int32_t index, double* x, double* y, double* z)
2478
{
2479
  if (!x || !y || !z) {
11!
2480
    set_errmsg(
×
2481
      "Invalid arguments passed to openmc_solidraytrace_plot_get_look_at");
2482
    return OPENMC_E_INVALID_ARGUMENT;
×
2483
  }
2484

2485
  SolidRayTracePlot* plt = nullptr;
11✔
2486
  int err = get_solidraytrace_plot_by_index(index, &plt);
11✔
2487
  if (err)
11!
2488
    return err;
2489

2490
  const auto& look_at = plt->look_at();
11✔
2491
  *x = look_at.x;
11✔
2492
  *y = look_at.y;
11✔
2493
  *z = look_at.z;
11✔
2494
  return 0;
11✔
2495
}
2496

2497
extern "C" int openmc_solidraytrace_plot_set_look_at(
11✔
2498
  int32_t index, double x, double y, double z)
2499
{
2500
  SolidRayTracePlot* plt = nullptr;
11✔
2501
  int err = get_solidraytrace_plot_by_index(index, &plt);
11✔
2502
  if (err)
11!
2503
    return err;
2504

2505
  plt->look_at() = {x, y, z};
11✔
2506
  return 0;
11✔
2507
}
2508

2509
extern "C" int openmc_solidraytrace_plot_get_up(
11✔
2510
  int32_t index, double* x, double* y, double* z)
2511
{
2512
  if (!x || !y || !z) {
11!
2513
    set_errmsg("Invalid arguments passed to openmc_solidraytrace_plot_get_up");
×
2514
    return OPENMC_E_INVALID_ARGUMENT;
×
2515
  }
2516

2517
  SolidRayTracePlot* plt = nullptr;
11✔
2518
  int err = get_solidraytrace_plot_by_index(index, &plt);
11✔
2519
  if (err)
11!
2520
    return err;
2521

2522
  const auto& up = plt->up();
11✔
2523
  *x = up.x;
11✔
2524
  *y = up.y;
11✔
2525
  *z = up.z;
11✔
2526
  return 0;
11✔
2527
}
2528

2529
extern "C" int openmc_solidraytrace_plot_set_up(
11✔
2530
  int32_t index, double x, double y, double z)
2531
{
2532
  SolidRayTracePlot* plt = nullptr;
11✔
2533
  int err = get_solidraytrace_plot_by_index(index, &plt);
11✔
2534
  if (err)
11!
2535
    return err;
2536

2537
  plt->up() = {x, y, z};
11✔
2538
  return 0;
11✔
2539
}
2540

2541
extern "C" int openmc_solidraytrace_plot_get_light_position(
11✔
2542
  int32_t index, double* x, double* y, double* z)
2543
{
2544
  if (!x || !y || !z) {
11!
2545
    set_errmsg("Invalid arguments passed to "
×
2546
               "openmc_solidraytrace_plot_get_light_position");
2547
    return OPENMC_E_INVALID_ARGUMENT;
×
2548
  }
2549

2550
  SolidRayTracePlot* plt = nullptr;
11✔
2551
  int err = get_solidraytrace_plot_by_index(index, &plt);
11✔
2552
  if (err)
11!
2553
    return err;
2554

2555
  const auto& light_position = plt->light_location();
11✔
2556
  *x = light_position.x;
11✔
2557
  *y = light_position.y;
11✔
2558
  *z = light_position.z;
11✔
2559
  return 0;
11✔
2560
}
2561

2562
extern "C" int openmc_solidraytrace_plot_set_light_position(
11✔
2563
  int32_t index, double x, double y, double z)
2564
{
2565
  SolidRayTracePlot* plt = nullptr;
11✔
2566
  int err = get_solidraytrace_plot_by_index(index, &plt);
11✔
2567
  if (err)
11!
2568
    return err;
2569

2570
  plt->light_location() = {x, y, z};
11✔
2571
  return 0;
11✔
2572
}
2573

2574
extern "C" int openmc_solidraytrace_plot_get_fov(int32_t index, double* fov)
11✔
2575
{
2576
  if (!fov) {
11!
2577
    set_errmsg("Invalid arguments passed to openmc_solidraytrace_plot_get_fov");
×
2578
    return OPENMC_E_INVALID_ARGUMENT;
×
2579
  }
2580

2581
  SolidRayTracePlot* plt = nullptr;
11✔
2582
  int err = get_solidraytrace_plot_by_index(index, &plt);
11✔
2583
  if (err)
11!
2584
    return err;
2585

2586
  *fov = plt->horizontal_field_of_view();
11✔
2587
  return 0;
11✔
2588
}
2589

2590
extern "C" int openmc_solidraytrace_plot_set_fov(int32_t index, double fov)
11✔
2591
{
2592
  SolidRayTracePlot* plt = nullptr;
11✔
2593
  int err = get_solidraytrace_plot_by_index(index, &plt);
11✔
2594
  if (err)
11!
2595
    return err;
2596

2597
  plt->horizontal_field_of_view() = fov;
11✔
2598
  return 0;
11✔
2599
}
2600

2601
extern "C" int openmc_solidraytrace_plot_update_view(int32_t index)
22✔
2602
{
2603
  SolidRayTracePlot* plt = nullptr;
22✔
2604
  int err = get_solidraytrace_plot_by_index(index, &plt);
22✔
2605
  if (err)
22!
2606
    return err;
2607

2608
  plt->update_view();
22✔
2609
  return 0;
2610
}
2611

2612
extern "C" int openmc_solidraytrace_plot_create_image(
22✔
2613
  int32_t index, uint8_t* data_out, int32_t width, int32_t height)
2614
{
2615
  if (!data_out || width <= 0 || height <= 0) {
22!
2616
    set_errmsg(
×
2617
      "Invalid arguments passed to openmc_solidraytrace_plot_create_image");
2618
    return OPENMC_E_INVALID_ARGUMENT;
×
2619
  }
2620

2621
  SolidRayTracePlot* plt = nullptr;
22✔
2622
  int err = get_solidraytrace_plot_by_index(index, &plt);
22✔
2623
  if (err)
22!
2624
    return err;
2625

2626
  if (plt->pixels()[0] != width || plt->pixels()[1] != height) {
22!
2627
    set_errmsg(
×
2628
      "Requested image size does not match SolidRayTracePlot pixel settings");
2629
    return OPENMC_E_INVALID_SIZE;
×
2630
  }
2631

2632
  ImageData data = plt->create_image();
22✔
2633
  if (static_cast<int32_t>(data.shape()[0]) != width ||
22!
2634
      static_cast<int32_t>(data.shape()[1]) != height) {
22!
2635
    set_errmsg("Unexpected image size from SolidRayTracePlot create_image");
×
2636
    return OPENMC_E_INVALID_SIZE;
×
2637
  }
2638

2639
  for (int32_t y = 0; y < height; ++y) {
154✔
2640
    for (int32_t x = 0; x < width; ++x) {
1,188✔
2641
      const auto& color = data(x, y);
1,056✔
2642
      size_t idx = (static_cast<size_t>(y) * width + x) * 3;
1,056✔
2643
      data_out[idx + 0] = color.red;
1,056✔
2644
      data_out[idx + 1] = color.green;
1,056✔
2645
      data_out[idx + 2] = color.blue;
1,056✔
2646
    }
2647
  }
2648

2649
  return 0;
2650
}
22✔
2651

2652
extern "C" int openmc_solidraytrace_plot_get_color(
11✔
2653
  int32_t index, int32_t id, uint8_t* r, uint8_t* g, uint8_t* b)
2654
{
2655
  if (!r || !g || !b) {
11!
2656
    set_errmsg(
×
2657
      "Invalid arguments passed to openmc_solidraytrace_plot_get_color");
2658
    return OPENMC_E_INVALID_ARGUMENT;
×
2659
  }
2660

2661
  SolidRayTracePlot* plt = nullptr;
11✔
2662
  int err = get_solidraytrace_plot_by_index(index, &plt);
11✔
2663
  if (err)
11!
2664
    return err;
2665

2666
  int32_t domain_index = -1;
11✔
2667
  err = map_phong_domain_id(plt, id, &domain_index);
11✔
2668
  if (err)
11!
2669
    return err;
2670

2671
  if (domain_index < 0 ||
11!
2672
      static_cast<size_t>(domain_index) >= plt->colors_.size()) {
11!
2673
    set_errmsg("Color index out of range for SolidRayTracePlot");
×
2674
    return OPENMC_E_OUT_OF_BOUNDS;
×
2675
  }
2676

2677
  const auto& color = plt->colors_[domain_index];
11✔
2678
  *r = color.red;
11✔
2679
  *g = color.green;
11✔
2680
  *b = color.blue;
11✔
2681
  return 0;
11✔
2682
}
2683

2684
extern "C" int openmc_solidraytrace_plot_get_diffuse_fraction(
11✔
2685
  int32_t index, double* diffuse_fraction)
2686
{
2687
  if (!diffuse_fraction) {
11!
2688
    set_errmsg("Invalid arguments passed to "
×
2689
               "openmc_solidraytrace_plot_get_diffuse_fraction");
2690
    return OPENMC_E_INVALID_ARGUMENT;
×
2691
  }
2692

2693
  SolidRayTracePlot* plt = nullptr;
11✔
2694
  int err = get_solidraytrace_plot_by_index(index, &plt);
11✔
2695
  if (err)
11!
2696
    return err;
2697

2698
  *diffuse_fraction = plt->diffuse_fraction();
11✔
2699
  return 0;
11✔
2700
}
2701

2702
extern "C" int openmc_solidraytrace_plot_set_diffuse_fraction(
11✔
2703
  int32_t index, double diffuse_fraction)
2704
{
2705
  SolidRayTracePlot* plt = nullptr;
11✔
2706
  int err = get_solidraytrace_plot_by_index(index, &plt);
11✔
2707
  if (err)
11!
2708
    return err;
2709

2710
  if (diffuse_fraction < 0.0 || diffuse_fraction > 1.0) {
11!
2711
    set_errmsg("Diffuse fraction must be between 0 and 1");
×
2712
    return OPENMC_E_INVALID_ARGUMENT;
×
2713
  }
2714

2715
  plt->diffuse_fraction() = diffuse_fraction;
11✔
2716
  return 0;
11✔
2717
}
2718

2719
} // namespace openmc
STATUS · Troubleshooting · Open an Issue · Sales · Support · CAREERS · ENTERPRISE · START FREE · SCHEDULE DEMO
ANNOUNCEMENTS · TWITTER · TOS & SLA · Supported CI Services · What's a CI service? · Automated Testing

© 2026 Coveralls, Inc