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

openmc-dev / openmc / 18554907707

16 Oct 2025 08:18AM UTC coverage: 81.605% (-0.3%) from 81.917%
18554907707

Pull #3598

github

web-flow
Merge a2c05a8c8 into b94b49611
Pull Request #3598: load mesh objects from weight_windows.h5 file

16605 of 23339 branches covered (71.15%)

Branch coverage included in aggregate %.

12 of 119 new or added lines in 3 files covered. (10.08%)

231 existing lines in 2 files now uncovered.

53725 of 62844 relevant lines covered (85.49%)

42595066.85 hits per line

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

67.75
/src/mesh.cpp
1
#include "openmc/mesh.h"
2
#include <algorithm> // for copy, equal, min, min_element
3
#include <cassert>
4
#define _USE_MATH_DEFINES // to make M_PI declared in Intel and MSVC compilers
5
#include <cmath>          // for ceil
6
#include <cstddef>        // for size_t
7
#include <string>
8

9
#ifdef _MSC_VER
10
#include <intrin.h> // for _InterlockedCompareExchange
11
#endif
12

13
#ifdef OPENMC_MPI
14
#include "mpi.h"
15
#endif
16

17
#include "xtensor/xadapt.hpp"
18
#include "xtensor/xbuilder.hpp"
19
#include "xtensor/xeval.hpp"
20
#include "xtensor/xmath.hpp"
21
#include "xtensor/xsort.hpp"
22
#include "xtensor/xtensor.hpp"
23
#include "xtensor/xview.hpp"
24
#include <fmt/core.h> // for fmt
25

26
#include "openmc/capi.h"
27
#include "openmc/constants.h"
28
#include "openmc/container_util.h"
29
#include "openmc/error.h"
30
#include "openmc/file_utils.h"
31
#include "openmc/geometry.h"
32
#include "openmc/hdf5_interface.h"
33
#include "openmc/material.h"
34
#include "openmc/memory.h"
35
#include "openmc/message_passing.h"
36
#include "openmc/openmp_interface.h"
37
#include "openmc/output.h"
38
#include "openmc/particle_data.h"
39
#include "openmc/plot.h"
40
#include "openmc/random_dist.h"
41
#include "openmc/search.h"
42
#include "openmc/settings.h"
43
#include "openmc/string_utils.h"
44
#include "openmc/tallies/filter.h"
45
#include "openmc/tallies/tally.h"
46
#include "openmc/timer.h"
47
#include "openmc/volume_calc.h"
48
#include "openmc/xml_interface.h"
49

50
#ifdef OPENMC_LIBMESH_ENABLED
51
#include "libmesh/mesh_modification.h"
52
#include "libmesh/mesh_tools.h"
53
#include "libmesh/numeric_vector.h"
54
#endif
55

56
#ifdef OPENMC_DAGMC_ENABLED
57
#include "moab/FileOptions.hpp"
58
#endif
59

60
namespace openmc {
61

62
//==============================================================================
63
// Global variables
64
//==============================================================================
65

66
#ifdef OPENMC_LIBMESH_ENABLED
67
const bool LIBMESH_ENABLED = true;
68
#else
69
const bool LIBMESH_ENABLED = false;
70
#endif
71

72
// Value used to indicate an empty slot in the hash table. We use -2 because
73
// the value -1 is used to indicate a void material.
74
constexpr int32_t EMPTY = -2;
75

76
namespace model {
77

78
std::unordered_map<int32_t, int32_t> mesh_map;
79
vector<unique_ptr<Mesh>> meshes;
80

81
} // namespace model
82

83
#ifdef OPENMC_LIBMESH_ENABLED
84
namespace settings {
85
unique_ptr<libMesh::LibMeshInit> libmesh_init;
86
const libMesh::Parallel::Communicator* libmesh_comm {nullptr};
87
} // namespace settings
88
#endif
89

90
//==============================================================================
91
// Helper functions
92
//==============================================================================
93

94
//! Update an intersection point if the given candidate is closer.
95
//
96
//! The first 6 arguments are coordinates for the starting point of a particle
97
//! and its intersection with a mesh surface.  If the distance between these
98
//! two points is shorter than the given `min_distance`, then the `r` argument
99
//! will be updated to match the intersection point, and `min_distance` will
100
//! also be updated.
101

102
inline bool check_intersection_point(double x1, double x0, double y1, double y0,
103
  double z1, double z0, Position& r, double& min_distance)
104
{
105
  double dist =
106
    std::pow(x1 - x0, 2) + std::pow(y1 - y0, 2) + std::pow(z1 - z0, 2);
107
  if (dist < min_distance) {
108
    r.x = x1;
109
    r.y = y1;
110
    r.z = z1;
111
    min_distance = dist;
112
    return true;
113
  }
114
  return false;
115
}
116

117
//! Atomic compare-and-swap for signed 32-bit integer
118
//
119
//! \param[in,out] ptr Pointer to value to update
120
//! \param[in,out] expected Value to compare to
121
//! \param[in] desired If comparison is successful, value to update to
122
//! \return True if the comparison was successful and the value was updated
123
inline bool atomic_cas_int32(int32_t* ptr, int32_t& expected, int32_t desired)
1,396✔
124
{
125
#if defined(__GNUC__) || defined(__clang__)
126
  // For gcc/clang, use the __atomic_compare_exchange_n intrinsic
127
  return __atomic_compare_exchange_n(
1,396✔
128
    ptr, &expected, desired, false, __ATOMIC_SEQ_CST, __ATOMIC_SEQ_CST);
1,396✔
129

130
#elif defined(_MSC_VER)
131
  // For MSVC, use the _InterlockedCompareExchange intrinsic
132
  int32_t old_val =
133
    _InterlockedCompareExchange(reinterpret_cast<volatile long*>(ptr),
134
      static_cast<long>(desired), static_cast<long>(expected));
135
  return (old_val == expected);
136

137
#else
138
#error "No compare-and-swap implementation available for this compiler."
139
#endif
140
}
141

142
namespace detail {
143

144
//==============================================================================
145
// MaterialVolumes implementation
146
//==============================================================================
147

148
void MaterialVolumes::add_volume(
2,653,387✔
149
  int index_elem, int index_material, double volume)
150
{
151
  // This method handles adding elements to the materials hash table,
152
  // implementing open addressing with linear probing. Consistency across
153
  // multiple threads is handled by with an atomic compare-and-swap operation.
154
  // Ideally, we would use #pragma omp atomic compare, but it was introduced in
155
  // OpenMP 5.1 and is not widely supported yet.
156

157
  // Loop for linear probing
158
  for (int attempt = 0; attempt < table_size_; ++attempt) {
2,659,073!
159
    // Determine slot to check, making sure it is positive
160
    int slot = (index_material + attempt) % table_size_;
2,659,073✔
161
    if (slot < 0)
2,659,073✔
162
      slot += table_size_;
208,978✔
163
    int32_t* slot_ptr = &this->materials(index_elem, slot);
2,659,073✔
164

165
    // Non-atomic read of current material
166
    int32_t current_val = *slot_ptr;
2,659,073✔
167

168
    // Found the desired material; accumulate volume
169
    if (current_val == index_material) {
2,659,073✔
170
#pragma omp atomic
1,447,002✔
171
      this->volumes(index_elem, slot) += volume;
2,651,991✔
172
      return;
2,651,991✔
173
    }
174

175
    // Slot appears to be empty; attempt to claim
176
    if (current_val == EMPTY) {
7,082✔
177
      // Attempt compare-and-swap from EMPTY to index_material
178
      int32_t expected_val = EMPTY;
1,396✔
179
      bool claimed_slot =
180
        atomic_cas_int32(slot_ptr, expected_val, index_material);
1,396✔
181

182
      // If we claimed the slot or another thread claimed it but the same
183
      // material was inserted, proceed to accumulate
184
      if (claimed_slot || (expected_val == index_material)) {
1,396!
185
#pragma omp atomic
768✔
186
        this->volumes(index_elem, slot) += volume;
1,396✔
187
        return;
1,396✔
188
      }
189
    }
190
  }
191

192
  // If table is full, set a flag that can be checked later
193
  table_full_ = true;
×
194
}
195

196
void MaterialVolumes::add_volume_unsafe(
×
197
  int index_elem, int index_material, double volume)
198
{
199
  // Linear probe
200
  for (int attempt = 0; attempt < table_size_; ++attempt) {
×
201
    // Determine slot to check, making sure it is positive
202
    int slot = (index_material + attempt) % table_size_;
×
203
    if (slot < 0)
×
204
      slot += table_size_;
×
205

206
    // Read current material
207
    int32_t current_val = this->materials(index_elem, slot);
×
208

209
    // Found the desired material; accumulate volume
210
    if (current_val == index_material) {
×
211
      this->volumes(index_elem, slot) += volume;
×
212
      return;
×
213
    }
214

215
    // Claim empty slot
216
    if (current_val == EMPTY) {
×
217
      this->materials(index_elem, slot) = index_material;
×
218
      this->volumes(index_elem, slot) += volume;
×
219
      return;
×
220
    }
221
  }
222

223
  // If table is full, set a flag that can be checked later
224
  table_full_ = true;
×
225
}
226

227
} // namespace detail
228

229
//==============================================================================
230
// Mesh implementation
231
//==============================================================================
232

233
Mesh::Mesh(pugi::xml_node node)
2,711✔
234
{
235
  // Read mesh id
236
  id_ = std::stoi(get_node_value(node, "id"));
2,711✔
237
  if (check_for_node(node, "name"))
2,711✔
238
    name_ = get_node_value(node, "name");
16✔
239
}
2,711✔
240

NEW
241
Mesh::Mesh(hid_t group)
×
242
{
243
  // Read mesh ID
NEW
244
  read_attribute(group, "id", id_);
×
245

246
  // Read mesh name
NEW
247
  if (object_exists(group, "name")) {
×
NEW
248
    read_dataset(group, "name", name_);
×
249
  }
NEW
250
}
×
251

252
void Mesh::set_id(int32_t id)
1✔
253
{
254
  assert(id >= 0 || id == C_NONE);
1!
255

256
  // Clear entry in mesh map in case one was already assigned
257
  if (id_ != C_NONE) {
1!
258
    model::mesh_map.erase(id_);
×
259
    id_ = C_NONE;
×
260
  }
261

262
  // Ensure no other mesh has the same ID
263
  if (model::mesh_map.find(id) != model::mesh_map.end()) {
1!
264
    throw std::runtime_error {
×
265
      fmt::format("Two meshes have the same ID: {}", id)};
×
266
  }
267

268
  // If no ID is specified, auto-assign the next ID in the sequence
269
  if (id == C_NONE) {
1!
270
    id = 0;
1✔
271
    for (const auto& m : model::meshes) {
3✔
272
      id = std::max(id, m->id_);
2✔
273
    }
274
    ++id;
1✔
275
  }
276

277
  // Update ID and entry in the mesh map
278
  id_ = id;
1✔
279
  model::mesh_map[id] = model::meshes.size() - 1;
1✔
280
}
1✔
281

282
vector<double> Mesh::volumes() const
252✔
283
{
284
  vector<double> volumes(n_bins());
252✔
285
  for (int i = 0; i < n_bins(); i++) {
1,125,127✔
286
    volumes[i] = this->volume(i);
1,124,875✔
287
  }
288
  return volumes;
252✔
UNCOV
289
}
×
290

291
void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
165✔
292
  int32_t* materials, double* volumes) const
293
{
294
  if (mpi::master) {
165!
295
    header("MESH MATERIAL VOLUMES CALCULATION", 7);
165✔
296
  }
297
  write_message(7, "Number of mesh elements = {}", n_bins());
165✔
298
  write_message(7, "Number of rays (x) = {}", nx);
165✔
299
  write_message(7, "Number of rays (y) = {}", ny);
165✔
300
  write_message(7, "Number of rays (z) = {}", nz);
165✔
301
  int64_t n_total = static_cast<int64_t>(nx) * ny +
165✔
302
                    static_cast<int64_t>(ny) * nz +
165✔
303
                    static_cast<int64_t>(nx) * nz;
165✔
304
  write_message(7, "Total number of rays = {}", n_total);
165✔
305
  write_message(7, "Table size per mesh element = {}", table_size);
165✔
306

307
  Timer timer;
165✔
308
  timer.start();
165✔
309

310
  // Create object for keeping track of materials/volumes
311
  detail::MaterialVolumes result(materials, volumes, table_size);
165✔
312

313
  // Determine bounding box
314
  auto bbox = this->bounding_box();
165✔
315

316
  std::array<int, 3> n_rays = {nx, ny, nz};
165✔
317

318
  // Determine effective width of rays
319
  Position width((nx > 0) ? (bbox.xmax - bbox.xmin) / nx : 0.0,
297✔
320
    (ny > 0) ? (bbox.ymax - bbox.ymin) / ny : 0.0,
319✔
321
    (nz > 0) ? (bbox.zmax - bbox.zmin) / nz : 0.0);
165✔
322

323
  // Set flag for mesh being contained within model
324
  bool out_of_model = false;
165✔
325

326
#pragma omp parallel
90✔
327
  {
328
    // Preallocate vector for mesh indices and length fractions and particle
329
    std::vector<int> bins;
75✔
330
    std::vector<double> length_fractions;
75✔
331
    Particle p;
75✔
332

333
    SourceSite site;
75✔
334
    site.E = 1.0;
75✔
335
    site.particle = ParticleType::neutron;
75✔
336

337
    for (int axis = 0; axis < 3; ++axis) {
300✔
338
      // Set starting position and direction
339
      site.r = {0.0, 0.0, 0.0};
225✔
340
      site.r[axis] = bbox.min()[axis];
225✔
341
      site.u = {0.0, 0.0, 0.0};
225✔
342
      site.u[axis] = 1.0;
225✔
343

344
      // Determine width of rays and number of rays in other directions
345
      int ax1 = (axis + 1) % 3;
225✔
346
      int ax2 = (axis + 2) % 3;
225✔
347
      double min1 = bbox.min()[ax1];
225✔
348
      double min2 = bbox.min()[ax2];
225✔
349
      double d1 = width[ax1];
225✔
350
      double d2 = width[ax2];
225✔
351
      int n1 = n_rays[ax1];
225✔
352
      int n2 = n_rays[ax2];
225✔
353
      if (n1 == 0 || n2 == 0) {
225✔
354
        continue;
60✔
355
      }
356

357
      // Divide rays in first direction over MPI processes by computing starting
358
      // and ending indices
359
      int min_work = n1 / mpi::n_procs;
165✔
360
      int remainder = n1 % mpi::n_procs;
165✔
361
      int n1_local = (mpi::rank < remainder) ? min_work + 1 : min_work;
165!
362
      int i1_start = mpi::rank * min_work + std::min(mpi::rank, remainder);
165✔
363
      int i1_end = i1_start + n1_local;
165✔
364

365
      // Loop over rays on face of bounding box
366
#pragma omp for collapse(2)
367
      for (int i1 = i1_start; i1 < i1_end; ++i1) {
8,900✔
368
        for (int i2 = 0; i2 < n2; ++i2) {
521,010✔
369
          site.r[ax1] = min1 + (i1 + 0.5) * d1;
512,275✔
370
          site.r[ax2] = min2 + (i2 + 0.5) * d2;
512,275✔
371

372
          p.from_source(&site);
512,275✔
373

374
          // Determine particle's location
375
          if (!exhaustive_find_cell(p)) {
512,275✔
376
            out_of_model = true;
39,930✔
377
            continue;
39,930✔
378
          }
379

380
          // Set birth cell attribute
381
          if (p.cell_born() == C_NONE)
472,345!
382
            p.cell_born() = p.lowest_coord().cell();
472,345✔
383

384
          // Initialize last cells from current cell
385
          for (int j = 0; j < p.n_coord(); ++j) {
944,690✔
386
            p.cell_last(j) = p.coord(j).cell();
472,345✔
387
          }
388
          p.n_coord_last() = p.n_coord();
472,345✔
389

390
          while (true) {
391
            // Ray trace from r_start to r_end
392
            Position r0 = p.r();
987,435✔
393
            double max_distance = bbox.max()[axis] - r0[axis];
987,435✔
394

395
            // Find the distance to the nearest boundary
396
            BoundaryInfo boundary = distance_to_boundary(p);
987,435✔
397

398
            // Advance particle forward
399
            double distance = std::min(boundary.distance(), max_distance);
987,435✔
400
            p.move_distance(distance);
987,435✔
401

402
            // Determine what mesh elements were crossed by particle
403
            bins.clear();
987,435✔
404
            length_fractions.clear();
987,435✔
405
            this->bins_crossed(r0, p.r(), p.u(), bins, length_fractions);
987,435✔
406

407
            // Add volumes to any mesh elements that were crossed
408
            int i_material = p.material();
987,435✔
409
            if (i_material != C_NONE) {
987,435✔
410
              i_material = model::materials[i_material]->id();
881,745✔
411
            }
412
            for (int i_bin = 0; i_bin < bins.size(); i_bin++) {
2,193,520✔
413
              int mesh_index = bins[i_bin];
1,206,085✔
414
              double length = distance * length_fractions[i_bin];
1,206,085✔
415

416
              // Add volume to result
417
              result.add_volume(mesh_index, i_material, length * d1 * d2);
1,206,085✔
418
            }
419

420
            if (distance == max_distance)
987,435✔
421
              break;
472,345✔
422

423
            // cross next geometric surface
424
            for (int j = 0; j < p.n_coord(); ++j) {
1,030,180✔
425
              p.cell_last(j) = p.coord(j).cell();
515,090✔
426
            }
427
            p.n_coord_last() = p.n_coord();
515,090✔
428

429
            // Set surface that particle is on and adjust coordinate levels
430
            p.surface() = boundary.surface();
515,090✔
431
            p.n_coord() = boundary.coord_level();
515,090✔
432

433
            if (boundary.lattice_translation()[0] != 0 ||
515,090✔
434
                boundary.lattice_translation()[1] != 0 ||
1,030,180!
435
                boundary.lattice_translation()[2] != 0) {
515,090!
436
              // Particle crosses lattice boundary
437
              cross_lattice(p, boundary);
×
438
            } else {
439
              // Particle crosses surface
440
              const auto& surf {model::surfaces[p.surface_index()].get()};
515,090✔
441
              p.cross_surface(*surf);
515,090✔
442
            }
443
          }
515,090✔
444
        }
445
      }
446
    }
447
  }
75✔
448

449
  // Check for errors
450
  if (out_of_model) {
165✔
451
    throw std::runtime_error("Mesh not fully contained in geometry.");
11✔
452
  } else if (result.table_full()) {
154!
UNCOV
453
    throw std::runtime_error("Maximum number of materials for mesh material "
×
UNCOV
454
                             "volume calculation insufficient.");
×
455
  }
456

457
  // Compute time for raytracing
458
  double t_raytrace = timer.elapsed();
154✔
459

460
#ifdef OPENMC_MPI
461
  // Combine results from multiple MPI processes
462
  if (mpi::n_procs > 1) {
70!
463
    int total = this->n_bins() * table_size;
×
464
    if (mpi::master) {
×
465
      // Allocate temporary buffer for receiving data
466
      std::vector<int32_t> mats(total);
×
467
      std::vector<double> vols(total);
×
468

469
      for (int i = 1; i < mpi::n_procs; ++i) {
×
470
        // Receive material indices and volumes from process i
471
        MPI_Recv(mats.data(), total, MPI_INT32_T, i, i, mpi::intracomm,
×
472
          MPI_STATUS_IGNORE);
473
        MPI_Recv(vols.data(), total, MPI_DOUBLE, i, i, mpi::intracomm,
×
474
          MPI_STATUS_IGNORE);
475

476
        // Combine with existing results; we can call thread unsafe version of
477
        // add_volume because each thread is operating on a different element
478
#pragma omp for
479
        for (int index_elem = 0; index_elem < n_bins(); ++index_elem) {
×
480
          for (int k = 0; k < table_size; ++k) {
×
481
            int index = index_elem * table_size + k;
482
            if (mats[index] != EMPTY) {
×
483
              result.add_volume_unsafe(index_elem, mats[index], vols[index]);
×
484
            }
485
          }
486
        }
487
      }
488
    } else {
489
      // Send material indices and volumes to process 0
490
      MPI_Send(materials, total, MPI_INT32_T, 0, mpi::rank, mpi::intracomm);
×
491
      MPI_Send(volumes, total, MPI_DOUBLE, 0, mpi::rank, mpi::intracomm);
×
492
    }
493
  }
494

495
  // Report time for MPI communication
496
  double t_mpi = timer.elapsed() - t_raytrace;
70✔
497
#else
498
  double t_mpi = 0.0;
84✔
499
#endif
500

501
  // Normalize based on known volumes of elements
502
  for (int i = 0; i < this->n_bins(); ++i) {
1,023✔
503
    // Estimated total volume in element i
504
    double volume = 0.0;
869✔
505
    for (int j = 0; j < table_size; ++j) {
7,821✔
506
      volume += result.volumes(i, j);
6,952✔
507
    }
508
    // Renormalize volumes based on known volume of element i
509
    double norm = this->volume(i) / volume;
869✔
510
    for (int j = 0; j < table_size; ++j) {
7,821✔
511
      result.volumes(i, j) *= norm;
6,952✔
512
    }
513
  }
514

515
  // Get total time and normalization time
516
  timer.stop();
154✔
517
  double t_total = timer.elapsed();
154✔
518
  double t_norm = t_total - t_raytrace - t_mpi;
154✔
519

520
  // Show timing statistics
521
  if (settings::verbosity < 7 || !mpi::master)
154!
522
    return;
44✔
523
  header("Timing Statistics", 7);
110✔
524
  fmt::print(" Total time elapsed            = {:.4e} seconds\n", t_total);
110✔
525
  fmt::print("   Ray tracing                 = {:.4e} seconds\n", t_raytrace);
110✔
526
  fmt::print("   MPI communication           = {:.4e} seconds\n", t_mpi);
110✔
527
  fmt::print("   Normalization               = {:.4e} seconds\n", t_norm);
90✔
528
  fmt::print(" Calculation rate              = {:.4e} rays/seconds\n",
90✔
529
    n_total / t_raytrace);
110✔
530
  fmt::print(" Calculation rate (per thread) = {:.4e} rays/seconds\n",
90✔
531
    n_total / (t_raytrace * mpi::n_procs * num_threads()));
110✔
532
  std::fflush(stdout);
110✔
533
}
534

535
void Mesh::to_hdf5(hid_t group) const
2,650✔
536
{
537
  // Create group for mesh
538
  std::string group_name = fmt::format("mesh {}", id_);
4,800✔
539
  hid_t mesh_group = create_group(group, group_name.c_str());
2,650✔
540

541
  // Write mesh type
542
  write_dataset(mesh_group, "type", this->get_mesh_type());
2,650✔
543

544
  // Write mesh ID
545
  write_attribute(mesh_group, "id", id_);
2,650✔
546

547
  // Write mesh name
548
  write_dataset(mesh_group, "name", name_);
2,650✔
549

550
  // Write mesh data
551
  this->to_hdf5_inner(mesh_group);
2,650✔
552

553
  // Close group
554
  close_group(mesh_group);
2,650✔
555
}
2,650✔
556

557
//==============================================================================
558
// Structured Mesh implementation
559
//==============================================================================
560

561
std::string StructuredMesh::bin_label(int bin) const
5,127,555✔
562
{
563
  MeshIndex ijk = get_indices_from_bin(bin);
5,127,555✔
564

565
  if (n_dimension_ > 2) {
5,127,555✔
566
    return fmt::format("Mesh Index ({}, {}, {})", ijk[0], ijk[1], ijk[2]);
10,225,608✔
567
  } else if (n_dimension_ > 1) {
14,751✔
568
    return fmt::format("Mesh Index ({}, {})", ijk[0], ijk[1]);
28,952✔
569
  } else {
570
    return fmt::format("Mesh Index ({})", ijk[0]);
550✔
571
  }
572
}
573

574
xt::xtensor<int, 1> StructuredMesh::get_x_shape() const
2,300✔
575
{
576
  // because method is const, shape_ is const as well and can't be adapted
577
  auto tmp_shape = shape_;
2,300✔
578
  return xt::adapt(tmp_shape, {n_dimension_});
4,600✔
579
}
580

581
Position StructuredMesh::sample_element(
1,419,966✔
582
  const MeshIndex& ijk, uint64_t* seed) const
583
{
584
  // lookup the lower/upper bounds for the mesh element
585
  double x_min = negative_grid_boundary(ijk, 0);
1,419,966✔
586
  double x_max = positive_grid_boundary(ijk, 0);
1,419,966✔
587

588
  double y_min = (n_dimension_ >= 2) ? negative_grid_boundary(ijk, 1) : 0.0;
1,419,966!
589
  double y_max = (n_dimension_ >= 2) ? positive_grid_boundary(ijk, 1) : 0.0;
1,419,966!
590

591
  double z_min = (n_dimension_ == 3) ? negative_grid_boundary(ijk, 2) : 0.0;
1,419,966!
592
  double z_max = (n_dimension_ == 3) ? positive_grid_boundary(ijk, 2) : 0.0;
1,419,966!
593

594
  return {x_min + (x_max - x_min) * prn(seed),
1,419,966✔
595
    y_min + (y_max - y_min) * prn(seed), z_min + (z_max - z_min) * prn(seed)};
1,419,966✔
596
}
597

598
//==============================================================================
599
// Unstructured Mesh implementation
600
//==============================================================================
601

602
UnstructuredMesh::UnstructuredMesh(pugi::xml_node node) : Mesh(node)
46✔
603
{
604
  n_dimension_ = 3;
46✔
605

606
  // check the mesh type
607
  if (check_for_node(node, "type")) {
46!
608
    auto temp = get_node_value(node, "type", true, true);
46!
609
    if (temp != mesh_type) {
46!
UNCOV
610
      fatal_error(fmt::format("Invalid mesh type: {}", temp));
×
611
    }
612
  }
46✔
613

614
  // check if a length unit multiplier was specified
615
  if (check_for_node(node, "length_multiplier")) {
46!
UNCOV
616
    length_multiplier_ = std::stod(get_node_value(node, "length_multiplier"));
×
617
  }
618

619
  // get the filename of the unstructured mesh to load
620
  if (check_for_node(node, "filename")) {
46!
621
    filename_ = get_node_value(node, "filename");
46!
622
    if (!file_exists(filename_)) {
46!
UNCOV
623
      fatal_error("Mesh file '" + filename_ + "' does not exist!");
×
624
    }
625
  } else {
UNCOV
626
    fatal_error(fmt::format(
×
UNCOV
627
      "No filename supplied for unstructured mesh with ID: {}", id_));
×
628
  }
629

630
  if (check_for_node(node, "options")) {
46!
631
    options_ = get_node_value(node, "options");
16!
632
  }
633

634
  // check if mesh tally data should be written with
635
  // statepoint files
636
  if (check_for_node(node, "output")) {
46!
UNCOV
637
    output_ = get_node_value_bool(node, "output");
×
638
  }
639
}
46✔
640

NEW
UNCOV
641
UnstructuredMesh::UnstructuredMesh(hid_t group) : Mesh(group)
×
642
{
NEW
UNCOV
643
  n_dimension_ = 3;
×
644

645
  // check the mesh type
NEW
UNCOV
646
  if (object_exists(group, "type")) {
×
NEW
UNCOV
647
    std::string temp;
×
NEW
UNCOV
648
    read_dataset(group, "type", temp);
×
NEW
UNCOV
649
    if (temp != mesh_type) {
×
NEW
UNCOV
650
      fatal_error(fmt::format("Invalid mesh type: {}", temp));
×
651
    }
NEW
652
  }
×
653

654
  // check if a length unit multiplier was specified
NEW
UNCOV
655
  if (object_exists(group, "length_multiplier")) {
×
NEW
UNCOV
656
    read_dataset(group, "length_multiplier", length_multiplier_);
×
657
  }
658

659
  // get the filename of the unstructured mesh to load
NEW
UNCOV
660
  if (object_exists(group, "filename")) {
×
NEW
UNCOV
661
    read_dataset(group, "filename", filename_);
×
NEW
UNCOV
662
    if (!file_exists(filename_)) {
×
NEW
UNCOV
663
      fatal_error("Mesh file '" + filename_ + "' does not exist!");
×
664
    }
665
  } else {
NEW
UNCOV
666
    fatal_error(fmt::format(
×
NEW
UNCOV
667
      "No filename supplied for unstructured mesh with ID: {}", id_));
×
668
  }
669

NEW
UNCOV
670
  if (attribute_exists(group, "options")) {
×
NEW
UNCOV
671
    read_attribute(group, "options", options_);
×
672
  }
673

674
  // check if mesh tally data should be written with
675
  // statepoint files
NEW
UNCOV
676
  if (attribute_exists(group, "output")) {
×
NEW
UNCOV
677
    read_attribute(group, "output", output_);
×
678
  }
NEW
679
}
×
680

681
void UnstructuredMesh::determine_bounds()
24✔
682
{
683
  double xmin = INFTY;
24✔
684
  double ymin = INFTY;
24✔
685
  double zmin = INFTY;
24✔
686
  double xmax = -INFTY;
24✔
687
  double ymax = -INFTY;
24✔
688
  double zmax = -INFTY;
24✔
689
  int n = this->n_vertices();
24!
690
  for (int i = 0; i < n; ++i) {
55,936✔
691
    auto v = this->vertex(i);
55,912!
692
    xmin = std::min(v.x, xmin);
55,912✔
693
    ymin = std::min(v.y, ymin);
55,912✔
694
    zmin = std::min(v.z, zmin);
55,912✔
695
    xmax = std::max(v.x, xmax);
55,912✔
696
    ymax = std::max(v.y, ymax);
55,912✔
697
    zmax = std::max(v.z, zmax);
55,912✔
698
  }
699
  lower_left_ = {xmin, ymin, zmin};
24!
700
  upper_right_ = {xmax, ymax, zmax};
24!
701
}
24✔
702

703
Position UnstructuredMesh::sample_tet(
601,230✔
704
  std::array<Position, 4> coords, uint64_t* seed) const
705
{
706
  // Uniform distribution
707
  double s = prn(seed);
601,230✔
708
  double t = prn(seed);
601,230✔
709
  double u = prn(seed);
601,230✔
710

711
  // From PyNE implementation of moab tet sampling C. Rocchini & P. Cignoni
712
  // (2000) Generating Random Points in a Tetrahedron, Journal of Graphics
713
  // Tools, 5:4, 9-12, DOI: 10.1080/10867651.2000.10487528
714
  if (s + t > 1) {
601,230✔
715
    s = 1.0 - s;
300,284✔
716
    t = 1.0 - t;
300,284✔
717
  }
718
  if (s + t + u > 1) {
601,230✔
719
    if (t + u > 1) {
400,918✔
720
      double old_t = t;
200,864✔
721
      t = 1.0 - u;
200,864✔
722
      u = 1.0 - s - old_t;
200,864✔
723
    } else if (t + u <= 1) {
200,054!
724
      double old_s = s;
200,054✔
725
      s = 1.0 - t - u;
200,054✔
726
      u = old_s + t + u - 1;
200,054✔
727
    }
728
  }
729
  return s * (coords[1] - coords[0]) + t * (coords[2] - coords[0]) +
1,202,460✔
730
         u * (coords[3] - coords[0]) + coords[0];
1,803,690✔
731
}
732

733
const std::string UnstructuredMesh::mesh_type = "unstructured";
734

735
std::string UnstructuredMesh::get_mesh_type() const
31✔
736
{
737
  return mesh_type;
31✔
738
}
739

UNCOV
740
void UnstructuredMesh::surface_bins_crossed(
×
741
  Position r0, Position r1, const Direction& u, vector<int>& bins) const
742
{
UNCOV
743
  fatal_error("Unstructured mesh surface tallies are not implemented.");
×
744
}
745

746
std::string UnstructuredMesh::bin_label(int bin) const
205,712✔
747
{
748
  return fmt::format("Mesh Index ({})", bin);
205,712!
749
};
750

751
void UnstructuredMesh::to_hdf5_inner(hid_t mesh_group) const
31✔
752
{
753
  write_dataset(mesh_group, "filename", filename_);
31!
754
  write_dataset(mesh_group, "library", this->library());
31!
755
  if (!options_.empty()) {
31✔
756
    write_attribute(mesh_group, "options", options_);
8!
757
  }
758

759
  if (length_multiplier_ > 0.0)
31!
UNCOV
760
    write_dataset(mesh_group, "length_multiplier", length_multiplier_);
×
761

762
  // write vertex coordinates
763
  xt::xtensor<double, 2> vertices({static_cast<size_t>(this->n_vertices()), 3});
31!
764
  for (int i = 0; i < this->n_vertices(); i++) {
70,260!
765
    auto v = this->vertex(i);
70,229!
766
    xt::view(vertices, i, xt::all()) = xt::xarray<double>({v.x, v.y, v.z});
70,229!
767
  }
768
  write_dataset(mesh_group, "vertices", vertices);
31!
769

770
  int num_elem_skipped = 0;
31✔
771

772
  // write element types and connectivity
773
  vector<double> volumes;
31✔
774
  xt::xtensor<int, 2> connectivity({static_cast<size_t>(this->n_bins()), 8});
31!
775
  xt::xtensor<int, 2> elem_types({static_cast<size_t>(this->n_bins()), 1});
31!
776
  for (int i = 0; i < this->n_bins(); i++) {
349,743!
777
    auto conn = this->connectivity(i);
349,712!
778

779
    volumes.emplace_back(this->volume(i));
349,712!
780

781
    // write linear tet element
782
    if (conn.size() == 4) {
349,712✔
783
      xt::view(elem_types, i, xt::all()) =
695,424!
784
        static_cast<int>(ElementType::LINEAR_TET);
695,424!
785
      xt::view(connectivity, i, xt::all()) =
695,424!
786
        xt::xarray<int>({conn[0], conn[1], conn[2], conn[3], -1, -1, -1, -1});
1,043,136!
787
      // write linear hex element
788
    } else if (conn.size() == 8) {
2,000!
789
      xt::view(elem_types, i, xt::all()) =
4,000!
790
        static_cast<int>(ElementType::LINEAR_HEX);
4,000!
791
      xt::view(connectivity, i, xt::all()) = xt::xarray<int>({conn[0], conn[1],
8,000!
792
        conn[2], conn[3], conn[4], conn[5], conn[6], conn[7]});
6,000!
793
    } else {
UNCOV
794
      num_elem_skipped++;
×
UNCOV
795
      xt::view(elem_types, i, xt::all()) =
×
UNCOV
796
        static_cast<int>(ElementType::UNSUPPORTED);
×
UNCOV
797
      xt::view(connectivity, i, xt::all()) = -1;
×
798
    }
799
  }
349,712✔
800

801
  // warn users that some elements were skipped
802
  if (num_elem_skipped > 0) {
31!
UNCOV
803
    warning(fmt::format("The connectivity of {} elements "
×
804
                        "on mesh {} were not written "
805
                        "because they are not of type linear tet/hex.",
UNCOV
806
      num_elem_skipped, this->id_));
×
807
  }
808

809
  write_dataset(mesh_group, "volumes", volumes);
31!
810
  write_dataset(mesh_group, "connectivity", connectivity);
31!
811
  write_dataset(mesh_group, "element_types", elem_types);
31!
812
}
31✔
813

814
void UnstructuredMesh::set_length_multiplier(double length_multiplier)
23✔
815
{
816
  length_multiplier_ = length_multiplier;
23✔
817
}
23✔
818

819
ElementType UnstructuredMesh::element_type(int bin) const
120,000✔
820
{
821
  auto conn = connectivity(bin);
120,000!
822

823
  if (conn.size() == 4)
120,000!
824
    return ElementType::LINEAR_TET;
120,000✔
UNCOV
825
  else if (conn.size() == 8)
×
UNCOV
826
    return ElementType::LINEAR_HEX;
×
827
  else
UNCOV
828
    return ElementType::UNSUPPORTED;
×
829
}
120,000✔
830

831
StructuredMesh::MeshIndex StructuredMesh::get_indices(
1,118,160,656✔
832
  Position r, bool& in_mesh) const
833
{
834
  MeshIndex ijk;
835
  in_mesh = true;
1,118,160,656✔
836
  for (int i = 0; i < n_dimension_; ++i) {
2,147,483,647✔
837
    ijk[i] = get_index_in_direction(r[i], i);
2,147,483,647✔
838

839
    if (ijk[i] < 1 || ijk[i] > shape_[i])
2,147,483,647✔
840
      in_mesh = false;
102,385,917✔
841
  }
842
  return ijk;
1,118,160,656✔
843
}
844

845
int StructuredMesh::get_bin_from_indices(const MeshIndex& ijk) const
1,659,381,085✔
846
{
847
  switch (n_dimension_) {
1,659,381,085!
848
  case 1:
880,605✔
849
    return ijk[0] - 1;
880,605✔
850
  case 2:
70,207,324✔
851
    return (ijk[1] - 1) * shape_[0] + ijk[0] - 1;
70,207,324✔
852
  case 3:
1,588,293,156✔
853
    return ((ijk[2] - 1) * shape_[1] + (ijk[1] - 1)) * shape_[0] + ijk[0] - 1;
1,588,293,156✔
UNCOV
854
  default:
×
UNCOV
855
    throw std::runtime_error {"Invalid number of mesh dimensions"};
×
856
  }
857
}
858

859
StructuredMesh::MeshIndex StructuredMesh::get_indices_from_bin(int bin) const
7,763,685✔
860
{
861
  MeshIndex ijk;
862
  if (n_dimension_ == 1) {
7,763,685✔
863
    ijk[0] = bin + 1;
275✔
864
  } else if (n_dimension_ == 2) {
7,763,410✔
865
    ijk[0] = bin % shape_[0] + 1;
14,476✔
866
    ijk[1] = bin / shape_[0] + 1;
14,476✔
867
  } else if (n_dimension_ == 3) {
7,748,934!
868
    ijk[0] = bin % shape_[0] + 1;
7,748,934✔
869
    ijk[1] = (bin % (shape_[0] * shape_[1])) / shape_[0] + 1;
7,748,934✔
870
    ijk[2] = bin / (shape_[0] * shape_[1]) + 1;
7,748,934✔
871
  }
872
  return ijk;
7,763,685✔
873
}
874

875
int StructuredMesh::get_bin(Position r) const
247,671,287✔
876
{
877
  // Determine indices
878
  bool in_mesh;
879
  MeshIndex ijk = get_indices(r, in_mesh);
247,671,287✔
880
  if (!in_mesh)
247,671,287✔
881
    return -1;
20,990,292✔
882

883
  // Convert indices to bin
884
  return get_bin_from_indices(ijk);
226,680,995✔
885
}
886

887
int StructuredMesh::n_bins() const
1,139,306✔
888
{
889
  return std::accumulate(
1,139,306✔
890
    shape_.begin(), shape_.begin() + n_dimension_, 1, std::multiplies<>());
2,278,612✔
891
}
892

893
int StructuredMesh::n_surface_bins() const
380✔
894
{
895
  return 4 * n_dimension_ * n_bins();
380✔
896
}
897

UNCOV
898
xt::xtensor<double, 1> StructuredMesh::count_sites(
×
899
  const SourceSite* bank, int64_t length, bool* outside) const
900
{
901
  // Determine shape of array for counts
UNCOV
902
  std::size_t m = this->n_bins();
×
UNCOV
903
  vector<std::size_t> shape = {m};
×
904

905
  // Create array of zeros
UNCOV
906
  xt::xarray<double> cnt {shape, 0.0};
×
UNCOV
907
  bool outside_ = false;
×
908

UNCOV
909
  for (int64_t i = 0; i < length; i++) {
×
UNCOV
910
    const auto& site = bank[i];
×
911

912
    // determine scoring bin for entropy mesh
UNCOV
913
    int mesh_bin = get_bin(site.r);
×
914

915
    // if outside mesh, skip particle
UNCOV
916
    if (mesh_bin < 0) {
×
UNCOV
917
      outside_ = true;
×
UNCOV
918
      continue;
×
919
    }
920

921
    // Add to appropriate bin
UNCOV
922
    cnt(mesh_bin) += site.wgt;
×
923
  }
924

925
  // Create copy of count data. Since ownership will be acquired by xtensor,
926
  // std::allocator must be used to avoid Valgrind mismatched free() / delete
927
  // warnings.
UNCOV
928
  int total = cnt.size();
×
UNCOV
929
  double* cnt_reduced = std::allocator<double> {}.allocate(total);
×
930

931
#ifdef OPENMC_MPI
932
  // collect values from all processors
933
  MPI_Reduce(
×
934
    cnt.data(), cnt_reduced, total, MPI_DOUBLE, MPI_SUM, 0, mpi::intracomm);
935

936
  // Check if there were sites outside the mesh for any processor
937
  if (outside) {
×
938
    MPI_Reduce(&outside_, outside, 1, MPI_C_BOOL, MPI_LOR, 0, mpi::intracomm);
×
939
  }
940
#else
941
  std::copy(cnt.data(), cnt.data() + total, cnt_reduced);
×
942
  if (outside)
×
943
    *outside = outside_;
944
#endif
945

946
  // Adapt reduced values in array back into an xarray
UNCOV
947
  auto arr = xt::adapt(cnt_reduced, total, xt::acquire_ownership(), shape);
×
948
  xt::xarray<double> counts = arr;
×
949

UNCOV
950
  return counts;
×
951
}
×
952

953
// raytrace through the mesh. The template class T will do the tallying.
954
// A modern optimizing compiler can recognize the noop method of T and
955
// eliminate that call entirely.
956
template<class T>
957
void StructuredMesh::raytrace_mesh(
875,479,606✔
958
  Position r0, Position r1, const Direction& u, T tally) const
959
{
960
  // TODO: when c++-17 is available, use "if constexpr ()" to compile-time
961
  // enable/disable tally calls for now, T template type needs to provide both
962
  // surface and track methods, which might be empty. modern optimizing
963
  // compilers will (hopefully) eliminate the complete code (including
964
  // calculation of parameters) but for the future: be explicit
965

966
  // Compute the length of the entire track.
967
  double total_distance = (r1 - r0).norm();
875,479,606✔
968
  if (total_distance == 0.0 && settings::solver_type != SolverType::RANDOM_RAY)
875,479,606✔
969
    return;
11,632,807✔
970

971
  // keep a copy of the original global position to pass to get_indices,
972
  // which performs its own transformation to local coordinates
973
  Position global_r = r0;
863,846,799✔
974
  Position local_r = local_coords(r0);
863,846,799✔
975

976
  const int n = n_dimension_;
863,846,799✔
977

978
  // Flag if position is inside the mesh
979
  bool in_mesh;
980

981
  // Position is r = r0 + u * traveled_distance, start at r0
982
  double traveled_distance {0.0};
863,846,799✔
983

984
  // Calculate index of current cell. Offset the position a tiny bit in
985
  // direction of flight
986
  MeshIndex ijk = get_indices(global_r + TINY_BIT * u, in_mesh);
863,846,799✔
987

988
  // if track is very short, assume that it is completely inside one cell.
989
  // Only the current cell will score and no surfaces
990
  if (total_distance < 2 * TINY_BIT) {
863,846,799✔
991
    if (in_mesh) {
331,527✔
992
      tally.track(ijk, 1.0);
331,043✔
993
    }
994
    return;
331,527✔
995
  }
996

997
  // Calculate initial distances to next surfaces in all three dimensions
998
  std::array<MeshDistance, 3> distances;
1,727,030,544✔
999
  for (int k = 0; k < n; ++k) {
2,147,483,647✔
1000
    distances[k] = distance_to_grid_boundary(ijk, k, local_r, u, 0.0);
2,147,483,647✔
1001
  }
1002

1003
  // Loop until r = r1 is eventually reached
1004
  while (true) {
736,590,918✔
1005

1006
    if (in_mesh) {
1,600,106,190✔
1007

1008
      // find surface with minimal distance to current position
1009
      const auto k = std::min_element(distances.begin(), distances.end()) -
1,514,254,422✔
1010
                     distances.begin();
1,514,254,422✔
1011

1012
      // Tally track length delta since last step
1013
      tally.track(ijk,
1,514,254,422✔
1014
        (std::min(distances[k].distance, total_distance) - traveled_distance) /
1,514,254,422✔
1015
          total_distance);
1016

1017
      // update position and leave, if we have reached end position
1018
      traveled_distance = distances[k].distance;
1,514,254,422✔
1019
      if (traveled_distance >= total_distance)
1,514,254,422✔
1020
        return;
784,306,074✔
1021

1022
      // If we have not reached r1, we have hit a surface. Tally outward
1023
      // current
1024
      tally.surface(ijk, k, distances[k].max_surface, false);
729,948,348✔
1025

1026
      // Update cell and calculate distance to next surface in k-direction.
1027
      // The two other directions are still valid!
1028
      ijk[k] = distances[k].next_index;
729,948,348✔
1029
      distances[k] =
729,948,348✔
1030
        distance_to_grid_boundary(ijk, k, local_r, u, traveled_distance);
729,948,348✔
1031

1032
      // Check if we have left the interior of the mesh
1033
      in_mesh = ((ijk[k] >= 1) && (ijk[k] <= shape_[k]));
729,948,348✔
1034

1035
      // If we are still inside the mesh, tally inward current for the next
1036
      // cell
1037
      if (in_mesh)
729,948,348✔
1038
        tally.surface(ijk, k, !distances[k].max_surface, true);
716,792,597✔
1039

1040
    } else { // not inside mesh
1041

1042
      // For all directions outside the mesh, find the distance that we need
1043
      // to travel to reach the next surface. Use the largest distance, as
1044
      // only this will cross all outer surfaces.
1045
      int k_max {-1};
85,851,768✔
1046
      for (int k = 0; k < n; ++k) {
341,962,706✔
1047
        if ((ijk[k] < 1 || ijk[k] > shape_[k]) &&
349,953,923✔
1048
            (distances[k].distance > traveled_distance)) {
93,842,985✔
1049
          traveled_distance = distances[k].distance;
88,864,873✔
1050
          k_max = k;
88,864,873✔
1051
        }
1052
      }
1053
      // Assure some distance is traveled
1054
      if (k_max == -1) {
85,851,768✔
1055
        traveled_distance += TINY_BIT;
110✔
1056
      }
1057

1058
      // If r1 is not inside the mesh, exit here
1059
      if (traveled_distance >= total_distance)
85,851,768✔
1060
        return;
79,209,198✔
1061

1062
      // Calculate the new cell index and update all distances to next
1063
      // surfaces.
1064
      ijk = get_indices(global_r + (traveled_distance + TINY_BIT) * u, in_mesh);
6,642,570✔
1065
      for (int k = 0; k < n; ++k) {
26,361,742✔
1066
        distances[k] =
19,719,172✔
1067
          distance_to_grid_boundary(ijk, k, local_r, u, traveled_distance);
19,719,172✔
1068
      }
1069

1070
      // If inside the mesh, Tally inward current
1071
      if (in_mesh && k_max >= 0)
6,642,570!
1072
        tally.surface(ijk, k_max, !distances[k_max].max_surface, true);
6,223,255✔
1073
    }
1074
  }
1075
}
1076

1077
void StructuredMesh::bins_crossed(Position r0, Position r1, const Direction& u,
763,352,041✔
1078
  vector<int>& bins, vector<double>& lengths) const
1079
{
1080

1081
  // Helper tally class.
1082
  // stores a pointer to the mesh class and references to bins and lengths
1083
  // parameters. Performs the actual tally through the track method.
1084
  struct TrackAggregator {
1085
    TrackAggregator(
763,352,041✔
1086
      const StructuredMesh* _mesh, vector<int>& _bins, vector<double>& _lengths)
1087
      : mesh(_mesh), bins(_bins), lengths(_lengths)
763,352,041✔
1088
    {}
763,352,041✔
1089
    void surface(const MeshIndex& ijk, int k, bool max, bool inward) const {}
1,394,805,011✔
1090
    void track(const MeshIndex& ijk, double l) const
1,374,540,901✔
1091
    {
1092
      bins.push_back(mesh->get_bin_from_indices(ijk));
1,374,540,901✔
1093
      lengths.push_back(l);
1,374,540,901✔
1094
    }
1,374,540,901✔
1095

1096
    const StructuredMesh* mesh;
1097
    vector<int>& bins;
1098
    vector<double>& lengths;
1099
  };
1100

1101
  // Perform the mesh raytrace with the helper class.
1102
  raytrace_mesh(r0, r1, u, TrackAggregator(this, bins, lengths));
763,352,041✔
1103
}
763,352,041✔
1104

1105
void StructuredMesh::surface_bins_crossed(
112,127,565✔
1106
  Position r0, Position r1, const Direction& u, vector<int>& bins) const
1107
{
1108

1109
  // Helper tally class.
1110
  // stores a pointer to the mesh class and a reference to the bins parameter.
1111
  // Performs the actual tally through the surface method.
1112
  struct SurfaceAggregator {
1113
    SurfaceAggregator(const StructuredMesh* _mesh, vector<int>& _bins)
112,127,565✔
1114
      : mesh(_mesh), bins(_bins)
112,127,565✔
1115
    {}
112,127,565✔
1116
    void surface(const MeshIndex& ijk, int k, bool max, bool inward) const
58,159,189✔
1117
    {
1118
      int i_bin =
1119
        4 * mesh->n_dimension_ * mesh->get_bin_from_indices(ijk) + 4 * k;
58,159,189✔
1120
      if (max)
58,159,189✔
1121
        i_bin += 2;
29,051,440✔
1122
      if (inward)
58,159,189✔
1123
        i_bin += 1;
28,582,290✔
1124
      bins.push_back(i_bin);
58,159,189✔
1125
    }
58,159,189✔
1126
    void track(const MeshIndex& idx, double l) const {}
140,044,564✔
1127

1128
    const StructuredMesh* mesh;
1129
    vector<int>& bins;
1130
  };
1131

1132
  // Perform the mesh raytrace with the helper class.
1133
  raytrace_mesh(r0, r1, u, SurfaceAggregator(this, bins));
112,127,565✔
1134
}
112,127,565✔
1135

1136
//==============================================================================
1137
// RegularMesh implementation
1138
//==============================================================================
1139

1140
RegularMesh::RegularMesh(pugi::xml_node node) : StructuredMesh {node}
1,837✔
1141
{
1142
  // Determine number of dimensions for mesh
1143
  if (!check_for_node(node, "dimension")) {
1,837!
UNCOV
1144
    fatal_error("Must specify <dimension> on a regular mesh.");
×
1145
  }
1146

1147
  xt::xtensor<int, 1> shape = get_node_xarray<int>(node, "dimension");
1,837✔
1148
  int n = n_dimension_ = shape.size();
1,837✔
1149
  if (n != 1 && n != 2 && n != 3) {
1,837!
UNCOV
1150
    fatal_error("Mesh must be one, two, or three dimensions.");
×
1151
  }
1152
  std::copy(shape.begin(), shape.end(), shape_.begin());
1,837✔
1153

1154
  // Check that dimensions are all greater than zero
1155
  if (xt::any(shape <= 0)) {
1,837!
UNCOV
1156
    fatal_error("All entries on the <dimension> element for a tally "
×
1157
                "mesh must be positive.");
1158
  }
1159

1160
  // Check for lower-left coordinates
1161
  if (check_for_node(node, "lower_left")) {
1,837!
1162
    // Read mesh lower-left corner location
1163
    lower_left_ = get_node_xarray<double>(node, "lower_left");
1,837✔
1164
  } else {
UNCOV
1165
    fatal_error("Must specify <lower_left> on a mesh.");
×
1166
  }
1167

1168
  // Make sure lower_left and dimension match
1169
  if (shape.size() != lower_left_.size()) {
1,837!
UNCOV
1170
    fatal_error("Number of entries on <lower_left> must be the same "
×
1171
                "as the number of entries on <dimension>.");
1172
  }
1173

1174
  if (check_for_node(node, "width")) {
1,837✔
1175
    // Make sure one of upper-right or width were specified
1176
    if (check_for_node(node, "upper_right")) {
49!
UNCOV
1177
      fatal_error("Cannot specify both <upper_right> and <width> on a mesh.");
×
1178
    }
1179

1180
    width_ = get_node_xarray<double>(node, "width");
49✔
1181

1182
    // Check to ensure width has same dimensions
1183
    auto n = width_.size();
49✔
1184
    if (n != lower_left_.size()) {
49!
1185
      fatal_error("Number of entries on <width> must be the same as "
×
1186
                  "the number of entries on <lower_left>.");
1187
    }
1188

1189
    // Check for negative widths
1190
    if (xt::any(width_ < 0.0)) {
49!
1191
      fatal_error("Cannot have a negative <width> on a tally mesh.");
×
1192
    }
1193

1194
    // Set width and upper right coordinate
1195
    upper_right_ = xt::eval(lower_left_ + shape * width_);
49✔
1196

1197
  } else if (check_for_node(node, "upper_right")) {
1,788!
1198
    upper_right_ = get_node_xarray<double>(node, "upper_right");
1,788✔
1199

1200
    // Check to ensure width has same dimensions
1201
    auto n = upper_right_.size();
1,788✔
1202
    if (n != lower_left_.size()) {
1,788!
1203
      fatal_error("Number of entries on <upper_right> must be the "
×
1204
                  "same as the number of entries on <lower_left>.");
1205
    }
1206

1207
    // Check that upper-right is above lower-left
1208
    if (xt::any(upper_right_ < lower_left_)) {
1,788!
1209
      fatal_error("The <upper_right> coordinates must be greater than "
×
1210
                  "the <lower_left> coordinates on a tally mesh.");
1211
    }
1212

1213
    // Set width
1214
    width_ = xt::eval((upper_right_ - lower_left_) / shape);
1,788✔
1215
  } else {
1216
    fatal_error("Must specify either <upper_right> or <width> on a mesh.");
×
1217
  }
1218

1219
  // Set material volumes
1220
  volume_frac_ = 1.0 / xt::prod(shape)();
1,837✔
1221

1222
  element_volume_ = 1.0;
1,837✔
1223
  for (int i = 0; i < n_dimension_; i++) {
6,974✔
1224
    element_volume_ *= width_[i];
5,137✔
1225
  }
1226
}
1,837✔
1227

NEW
1228
RegularMesh::RegularMesh(hid_t group) : StructuredMesh {group}
×
1229
{
1230
  // Determine number of dimensions for mesh
NEW
1231
  if (!object_exists(group, "dimension")) {
×
NEW
1232
    fatal_error("Must specify <dimension> on a regular mesh.");
×
1233
  }
1234

NEW
1235
  xt::xtensor<int, 1> shape;
×
NEW
1236
  read_dataset(group, "dimension", shape);
×
NEW
1237
  int n = n_dimension_ = shape.size();
×
NEW
1238
  if (n != 1 && n != 2 && n != 3) {
×
NEW
1239
    fatal_error("Mesh must be one, two, or three dimensions.");
×
1240
  }
NEW
1241
  std::copy(shape.begin(), shape.end(), shape_.begin());
×
1242

1243
  // Check that dimensions are all greater than zero
NEW
1244
  if (xt::any(shape <= 0)) {
×
NEW
1245
    fatal_error("All entries on the <dimension> element for a tally "
×
1246
                "mesh must be positive.");
1247
  }
1248

1249
  // Check for lower-left coordinates
NEW
UNCOV
1250
  if (object_exists(group, "lower_left")) {
×
1251
    // Read mesh lower-left corner location
NEW
1252
    read_dataset(group, "lower_left", lower_left_);
×
1253
  } else {
NEW
UNCOV
1254
    fatal_error("Must specify <lower_left> on a mesh.");
×
1255
  }
1256

1257
  // Make sure lower_left and dimension match
NEW
1258
  if (shape.size() != lower_left_.size()) {
×
NEW
UNCOV
1259
    fatal_error("Number of entries on <lower_left> must be the same "
×
1260
                "as the number of entries on <dimension>.");
1261
  }
1262

NEW
UNCOV
1263
  if (object_exists(group, "upper_right")) {
×
1264

NEW
UNCOV
1265
    read_dataset(group, "upper_right", upper_right_);
×
1266

1267
    // Check to ensure width has same dimensions
NEW
UNCOV
1268
    auto n = upper_right_.size();
×
NEW
UNCOV
1269
    if (n != lower_left_.size()) {
×
NEW
UNCOV
1270
      fatal_error("Number of entries on <upper_right> must be the "
×
1271
                  "same as the number of entries on <lower_left>.");
1272
    }
1273

1274
    // Check that upper-right is above lower-left
NEW
UNCOV
1275
    if (xt::any(upper_right_ < lower_left_)) {
×
NEW
UNCOV
1276
      fatal_error("The <upper_right> coordinates must be greater than "
×
1277
                  "the <lower_left> coordinates on a tally mesh.");
1278
    }
1279

1280
    // Set width
NEW
UNCOV
1281
    width_ = xt::eval((upper_right_ - lower_left_) / shape);
×
1282
  } else {
NEW
1283
    fatal_error("Must specify either <upper_right> or <width> on a mesh.");
×
1284
  }
1285

1286
  // Set material volumes
NEW
1287
  volume_frac_ = 1.0 / xt::prod(shape)();
×
1288

NEW
1289
  element_volume_ = 1.0;
×
NEW
UNCOV
1290
  for (int i = 0; i < n_dimension_; i++) {
×
NEW
1291
    element_volume_ *= width_[i];
×
1292
  }
NEW
1293
}
×
1294

1295
int RegularMesh::get_index_in_direction(double r, int i) const
2,147,483,647✔
1296
{
1297
  return std::ceil((r - lower_left_[i]) / width_[i]);
2,147,483,647✔
1298
}
1299

1300
const std::string RegularMesh::mesh_type = "regular";
1301

1302
std::string RegularMesh::get_mesh_type() const
2,960✔
1303
{
1304
  return mesh_type;
2,960✔
1305
}
1306

1307
double RegularMesh::positive_grid_boundary(const MeshIndex& ijk, int i) const
1,406,200,992✔
1308
{
1309
  return lower_left_[i] + ijk[i] * width_[i];
1,406,200,992✔
1310
}
1311

1312
double RegularMesh::negative_grid_boundary(const MeshIndex& ijk, int i) const
1,344,486,367✔
1313
{
1314
  return lower_left_[i] + (ijk[i] - 1) * width_[i];
1,344,486,367✔
1315
}
1316

1317
StructuredMesh::MeshDistance RegularMesh::distance_to_grid_boundary(
2,147,483,647✔
1318
  const MeshIndex& ijk, int i, const Position& r0, const Direction& u,
1319
  double l) const
1320
{
1321
  MeshDistance d;
2,147,483,647✔
1322
  d.next_index = ijk[i];
2,147,483,647✔
1323
  if (std::abs(u[i]) < FP_PRECISION)
2,147,483,647✔
1324
    return d;
1,458,182✔
1325

1326
  d.max_surface = (u[i] > 0);
2,147,483,647✔
1327
  if (d.max_surface && (ijk[i] <= shape_[i])) {
2,147,483,647✔
1328
    d.next_index++;
1,401,941,094✔
1329
    d.distance = (positive_grid_boundary(ijk, i) - r0[i]) / u[i];
1,401,941,094✔
1330
  } else if (!d.max_surface && (ijk[i] >= 1)) {
1,361,387,142✔
1331
    d.next_index--;
1,340,226,469✔
1332
    d.distance = (negative_grid_boundary(ijk, i) - r0[i]) / u[i];
1,340,226,469✔
1333
  }
1334

1335
  return d;
2,147,483,647✔
1336
}
1337

1338
std::pair<vector<double>, vector<double>> RegularMesh::plot(
22✔
1339
  Position plot_ll, Position plot_ur) const
1340
{
1341
  // Figure out which axes lie in the plane of the plot.
1342
  array<int, 2> axes {-1, -1};
22✔
1343
  if (plot_ur.z == plot_ll.z) {
22!
1344
    axes[0] = 0;
22✔
1345
    if (n_dimension_ > 1)
22!
1346
      axes[1] = 1;
22✔
UNCOV
1347
  } else if (plot_ur.y == plot_ll.y) {
×
UNCOV
1348
    axes[0] = 0;
×
UNCOV
1349
    if (n_dimension_ > 2)
×
UNCOV
1350
      axes[1] = 2;
×
UNCOV
1351
  } else if (plot_ur.x == plot_ll.x) {
×
UNCOV
1352
    if (n_dimension_ > 1)
×
UNCOV
1353
      axes[0] = 1;
×
UNCOV
1354
    if (n_dimension_ > 2)
×
UNCOV
1355
      axes[1] = 2;
×
1356
  } else {
UNCOV
1357
    fatal_error("Can only plot mesh lines on an axis-aligned plot");
×
1358
  }
1359

1360
  // Get the coordinates of the mesh lines along both of the axes.
1361
  array<vector<double>, 2> axis_lines;
22✔
1362
  for (int i_ax = 0; i_ax < 2; ++i_ax) {
66✔
1363
    int axis = axes[i_ax];
44✔
1364
    if (axis == -1)
44!
UNCOV
1365
      continue;
×
1366
    auto& lines {axis_lines[i_ax]};
44✔
1367

1368
    double coord = lower_left_[axis];
44✔
1369
    for (int i = 0; i < shape_[axis] + 1; ++i) {
286✔
1370
      if (coord >= plot_ll[axis] && coord <= plot_ur[axis])
242!
1371
        lines.push_back(coord);
242✔
1372
      coord += width_[axis];
242✔
1373
    }
1374
  }
1375

1376
  return {axis_lines[0], axis_lines[1]};
44✔
1377
}
22✔
1378

1379
void RegularMesh::to_hdf5_inner(hid_t mesh_group) const
1,860✔
1380
{
1381
  write_dataset(mesh_group, "dimension", get_x_shape());
1,860✔
1382
  write_dataset(mesh_group, "lower_left", lower_left_);
1,860✔
1383
  write_dataset(mesh_group, "upper_right", upper_right_);
1,860✔
1384
  write_dataset(mesh_group, "width", width_);
1,860✔
1385
}
1,860✔
1386

1387
xt::xtensor<double, 1> RegularMesh::count_sites(
7,839✔
1388
  const SourceSite* bank, int64_t length, bool* outside) const
1389
{
1390
  // Determine shape of array for counts
1391
  std::size_t m = this->n_bins();
7,839✔
1392
  vector<std::size_t> shape = {m};
7,839✔
1393

1394
  // Create array of zeros
1395
  xt::xarray<double> cnt {shape, 0.0};
7,839✔
1396
  bool outside_ = false;
7,839✔
1397

1398
  for (int64_t i = 0; i < length; i++) {
7,675,290✔
1399
    const auto& site = bank[i];
7,667,451✔
1400

1401
    // determine scoring bin for entropy mesh
1402
    int mesh_bin = get_bin(site.r);
7,667,451✔
1403

1404
    // if outside mesh, skip particle
1405
    if (mesh_bin < 0) {
7,667,451!
UNCOV
1406
      outside_ = true;
×
UNCOV
1407
      continue;
×
1408
    }
1409

1410
    // Add to appropriate bin
1411
    cnt(mesh_bin) += site.wgt;
7,667,451✔
1412
  }
1413

1414
  // Create copy of count data. Since ownership will be acquired by xtensor,
1415
  // std::allocator must be used to avoid Valgrind mismatched free() / delete
1416
  // warnings.
1417
  int total = cnt.size();
7,839✔
1418
  double* cnt_reduced = std::allocator<double> {}.allocate(total);
7,839✔
1419

1420
#ifdef OPENMC_MPI
1421
  // collect values from all processors
1422
  MPI_Reduce(
3,615✔
1423
    cnt.data(), cnt_reduced, total, MPI_DOUBLE, MPI_SUM, 0, mpi::intracomm);
3,615✔
1424

1425
  // Check if there were sites outside the mesh for any processor
1426
  if (outside) {
3,615!
1427
    MPI_Reduce(&outside_, outside, 1, MPI_C_BOOL, MPI_LOR, 0, mpi::intracomm);
3,615✔
1428
  }
1429
#else
1430
  std::copy(cnt.data(), cnt.data() + total, cnt_reduced);
4,224✔
1431
  if (outside)
4,224!
1432
    *outside = outside_;
4,224✔
1433
#endif
1434

1435
  // Adapt reduced values in array back into an xarray
1436
  auto arr = xt::adapt(cnt_reduced, total, xt::acquire_ownership(), shape);
7,839✔
1437
  xt::xarray<double> counts = arr;
7,839✔
1438

1439
  return counts;
15,678✔
1440
}
7,839✔
1441

1442
double RegularMesh::volume(const MeshIndex& ijk) const
1,126,085✔
1443
{
1444
  return element_volume_;
1,126,085✔
1445
}
1446

1447
//==============================================================================
1448
// RectilinearMesh implementation
1449
//==============================================================================
1450

1451
RectilinearMesh::RectilinearMesh(pugi::xml_node node) : StructuredMesh {node}
114✔
1452
{
1453
  n_dimension_ = 3;
114✔
1454

1455
  grid_[0] = get_node_array<double>(node, "x_grid");
114✔
1456
  grid_[1] = get_node_array<double>(node, "y_grid");
114✔
1457
  grid_[2] = get_node_array<double>(node, "z_grid");
114✔
1458

1459
  if (int err = set_grid()) {
114!
UNCOV
1460
    fatal_error(openmc_err_msg);
×
1461
  }
1462
}
114✔
1463

NEW
UNCOV
1464
RectilinearMesh::RectilinearMesh(hid_t group) : StructuredMesh {group}
×
1465
{
NEW
UNCOV
1466
  n_dimension_ = 3;
×
1467

NEW
UNCOV
1468
  read_dataset(group, "x_grid", grid_[0]);
×
NEW
UNCOV
1469
  read_dataset(group, "y_grid", grid_[1]);
×
NEW
UNCOV
1470
  read_dataset(group, "z_grid", grid_[2]);
×
1471

NEW
UNCOV
1472
  if (int err = set_grid()) {
×
NEW
UNCOV
1473
    fatal_error(openmc_err_msg);
×
1474
  }
NEW
UNCOV
1475
}
×
1476

1477
const std::string RectilinearMesh::mesh_type = "rectilinear";
1478

1479
std::string RectilinearMesh::get_mesh_type() const
264✔
1480
{
1481
  return mesh_type;
264✔
1482
}
1483

1484
double RectilinearMesh::positive_grid_boundary(
26,505,963✔
1485
  const MeshIndex& ijk, int i) const
1486
{
1487
  return grid_[i][ijk[i]];
26,505,963✔
1488
}
1489

1490
double RectilinearMesh::negative_grid_boundary(
25,739,406✔
1491
  const MeshIndex& ijk, int i) const
1492
{
1493
  return grid_[i][ijk[i] - 1];
25,739,406✔
1494
}
1495

1496
StructuredMesh::MeshDistance RectilinearMesh::distance_to_grid_boundary(
53,602,087✔
1497
  const MeshIndex& ijk, int i, const Position& r0, const Direction& u,
1498
  double l) const
1499
{
1500
  MeshDistance d;
53,602,087✔
1501
  d.next_index = ijk[i];
53,602,087✔
1502
  if (std::abs(u[i]) < FP_PRECISION)
53,602,087✔
1503
    return d;
571,824✔
1504

1505
  d.max_surface = (u[i] > 0);
53,030,263✔
1506
  if (d.max_surface && (ijk[i] <= shape_[i])) {
53,030,263✔
1507
    d.next_index++;
26,505,963✔
1508
    d.distance = (positive_grid_boundary(ijk, i) - r0[i]) / u[i];
26,505,963✔
1509
  } else if (!d.max_surface && (ijk[i] > 0)) {
26,524,300✔
1510
    d.next_index--;
25,739,406✔
1511
    d.distance = (negative_grid_boundary(ijk, i) - r0[i]) / u[i];
25,739,406✔
1512
  }
1513
  return d;
53,030,263✔
1514
}
1515

1516
int RectilinearMesh::set_grid()
158✔
1517
{
1518
  shape_ = {static_cast<int>(grid_[0].size()) - 1,
158✔
1519
    static_cast<int>(grid_[1].size()) - 1,
158✔
1520
    static_cast<int>(grid_[2].size()) - 1};
158✔
1521

1522
  for (const auto& g : grid_) {
632✔
1523
    if (g.size() < 2) {
474!
UNCOV
1524
      set_errmsg("x-, y-, and z- grids for rectilinear meshes "
×
1525
                 "must each have at least 2 points");
UNCOV
1526
      return OPENMC_E_INVALID_ARGUMENT;
×
1527
    }
1528
    if (std::adjacent_find(g.begin(), g.end(), std::greater_equal<>()) !=
474✔
1529
        g.end()) {
948!
UNCOV
1530
      set_errmsg("Values in for x-, y-, and z- grids for "
×
1531
                 "rectilinear meshes must be sorted and unique.");
UNCOV
1532
      return OPENMC_E_INVALID_ARGUMENT;
×
1533
    }
1534
  }
1535

1536
  lower_left_ = {grid_[0].front(), grid_[1].front(), grid_[2].front()};
158✔
1537
  upper_right_ = {grid_[0].back(), grid_[1].back(), grid_[2].back()};
158✔
1538

1539
  return 0;
158✔
1540
}
1541

1542
int RectilinearMesh::get_index_in_direction(double r, int i) const
74,108,892✔
1543
{
1544
  return lower_bound_index(grid_[i].begin(), grid_[i].end(), r) + 1;
74,108,892✔
1545
}
1546

1547
std::pair<vector<double>, vector<double>> RectilinearMesh::plot(
11✔
1548
  Position plot_ll, Position plot_ur) const
1549
{
1550
  // Figure out which axes lie in the plane of the plot.
1551
  array<int, 2> axes {-1, -1};
11✔
1552
  if (plot_ur.z == plot_ll.z) {
11!
UNCOV
1553
    axes = {0, 1};
×
1554
  } else if (plot_ur.y == plot_ll.y) {
11!
1555
    axes = {0, 2};
11✔
1556
  } else if (plot_ur.x == plot_ll.x) {
×
UNCOV
1557
    axes = {1, 2};
×
1558
  } else {
UNCOV
1559
    fatal_error("Can only plot mesh lines on an axis-aligned plot");
×
1560
  }
1561

1562
  // Get the coordinates of the mesh lines along both of the axes.
1563
  array<vector<double>, 2> axis_lines;
11✔
1564
  for (int i_ax = 0; i_ax < 2; ++i_ax) {
33✔
1565
    int axis = axes[i_ax];
22✔
1566
    vector<double>& lines {axis_lines[i_ax]};
22✔
1567

1568
    for (auto coord : grid_[axis]) {
110✔
1569
      if (coord >= plot_ll[axis] && coord <= plot_ur[axis])
88!
1570
        lines.push_back(coord);
88✔
1571
    }
1572
  }
1573

1574
  return {axis_lines[0], axis_lines[1]};
22✔
1575
}
11✔
1576

1577
void RectilinearMesh::to_hdf5_inner(hid_t mesh_group) const
99✔
1578
{
1579
  write_dataset(mesh_group, "x_grid", grid_[0]);
99✔
1580
  write_dataset(mesh_group, "y_grid", grid_[1]);
99✔
1581
  write_dataset(mesh_group, "z_grid", grid_[2]);
99✔
1582
}
99✔
1583

1584
double RectilinearMesh::volume(const MeshIndex& ijk) const
132✔
1585
{
1586
  double vol {1.0};
132✔
1587

1588
  for (int i = 0; i < n_dimension_; i++) {
528✔
1589
    vol *= grid_[i][ijk[i]] - grid_[i][ijk[i] - 1];
396✔
1590
  }
1591
  return vol;
132✔
1592
}
1593

1594
//==============================================================================
1595
// CylindricalMesh implementation
1596
//==============================================================================
1597

1598
CylindricalMesh::CylindricalMesh(pugi::xml_node node)
390✔
1599
  : PeriodicStructuredMesh {node}
390✔
1600
{
1601
  n_dimension_ = 3;
390✔
1602
  grid_[0] = get_node_array<double>(node, "r_grid");
390✔
1603
  grid_[1] = get_node_array<double>(node, "phi_grid");
390✔
1604
  grid_[2] = get_node_array<double>(node, "z_grid");
390✔
1605
  origin_ = get_node_position(node, "origin");
390✔
1606

1607
  if (int err = set_grid()) {
390!
UNCOV
1608
    fatal_error(openmc_err_msg);
×
1609
  }
1610
}
390✔
1611

NEW
UNCOV
1612
CylindricalMesh::CylindricalMesh(hid_t group) : PeriodicStructuredMesh {group}
×
1613
{
NEW
UNCOV
1614
  n_dimension_ = 3;
×
NEW
UNCOV
1615
  read_dataset(group, "r_grid", grid_[0]);
×
NEW
UNCOV
1616
  read_dataset(group, "phi_grid", grid_[1]);
×
NEW
UNCOV
1617
  read_dataset(group, "z_grid", grid_[2]);
×
NEW
UNCOV
1618
  read_dataset(group, "origin", origin_);
×
1619

NEW
UNCOV
1620
  if (int err = set_grid()) {
×
NEW
UNCOV
1621
    fatal_error(openmc_err_msg);
×
1622
  }
NEW
UNCOV
1623
}
×
1624

1625
const std::string CylindricalMesh::mesh_type = "cylindrical";
1626

1627
std::string CylindricalMesh::get_mesh_type() const
473✔
1628
{
1629
  return mesh_type;
473✔
1630
}
1631

1632
StructuredMesh::MeshIndex CylindricalMesh::get_indices(
47,726,668✔
1633
  Position r, bool& in_mesh) const
1634
{
1635
  r = local_coords(r);
47,726,668✔
1636

1637
  Position mapped_r;
47,726,668✔
1638
  mapped_r[0] = std::hypot(r.x, r.y);
47,726,668✔
1639
  mapped_r[2] = r[2];
47,726,668✔
1640

1641
  if (mapped_r[0] < FP_PRECISION) {
47,726,668!
UNCOV
1642
    mapped_r[1] = 0.0;
×
1643
  } else {
1644
    mapped_r[1] = std::atan2(r.y, r.x);
47,726,668✔
1645
    if (mapped_r[1] < 0)
47,726,668✔
1646
      mapped_r[1] += 2 * M_PI;
23,872,431✔
1647
  }
1648

1649
  MeshIndex idx = StructuredMesh::get_indices(mapped_r, in_mesh);
47,726,668✔
1650

1651
  idx[1] = sanitize_phi(idx[1]);
47,726,668✔
1652

1653
  return idx;
47,726,668✔
1654
}
1655

1656
Position CylindricalMesh::sample_element(
88,110✔
1657
  const MeshIndex& ijk, uint64_t* seed) const
1658
{
1659
  double r_min = this->r(ijk[0] - 1);
88,110✔
1660
  double r_max = this->r(ijk[0]);
88,110✔
1661

1662
  double phi_min = this->phi(ijk[1] - 1);
88,110✔
1663
  double phi_max = this->phi(ijk[1]);
88,110✔
1664

1665
  double z_min = this->z(ijk[2] - 1);
88,110✔
1666
  double z_max = this->z(ijk[2]);
88,110✔
1667

1668
  double r_min_sq = r_min * r_min;
88,110✔
1669
  double r_max_sq = r_max * r_max;
88,110✔
1670
  double r = std::sqrt(uniform_distribution(r_min_sq, r_max_sq, seed));
88,110✔
1671
  double phi = uniform_distribution(phi_min, phi_max, seed);
88,110✔
1672
  double z = uniform_distribution(z_min, z_max, seed);
88,110✔
1673

1674
  double x = r * std::cos(phi);
88,110✔
1675
  double y = r * std::sin(phi);
88,110✔
1676

1677
  return origin_ + Position(x, y, z);
88,110✔
1678
}
1679

1680
double CylindricalMesh::find_r_crossing(
142,566,864✔
1681
  const Position& r, const Direction& u, double l, int shell) const
1682
{
1683

1684
  if ((shell < 0) || (shell > shape_[0]))
142,566,864!
1685
    return INFTY;
17,913,852✔
1686

1687
  // solve r.x^2 + r.y^2 == r0^2
1688
  // x^2 + 2*s*u*x + s^2*u^2 + s^2*v^2+2*s*v*y + y^2 -r0^2 = 0
1689
  // s^2 * (u^2 + v^2) + 2*s*(u*x+v*y) + x^2+y^2-r0^2 = 0
1690

1691
  const double r0 = grid_[0][shell];
124,653,012✔
1692
  if (r0 == 0.0)
124,653,012✔
1693
    return INFTY;
7,130,651✔
1694

1695
  const double denominator = u.x * u.x + u.y * u.y;
117,522,361✔
1696

1697
  // Direction of flight is in z-direction. Will never intersect r.
1698
  if (std::abs(denominator) < FP_PRECISION)
117,522,361✔
1699
    return INFTY;
58,960✔
1700

1701
  // inverse of dominator to help the compiler to speed things up
1702
  const double inv_denominator = 1.0 / denominator;
117,463,401✔
1703

1704
  const double p = (u.x * r.x + u.y * r.y) * inv_denominator;
117,463,401✔
1705
  double c = r.x * r.x + r.y * r.y - r0 * r0;
117,463,401✔
1706
  double D = p * p - c * inv_denominator;
117,463,401✔
1707

1708
  if (D < 0.0)
117,463,401✔
1709
    return INFTY;
9,733,570✔
1710

1711
  D = std::sqrt(D);
107,729,831✔
1712

1713
  // the solution -p - D is always smaller as -p + D : Check this one first
1714
  if (std::abs(c) <= RADIAL_MESH_TOL)
107,729,831✔
1715
    return INFTY;
6,611,374✔
1716

1717
  if (-p - D > l)
101,118,457✔
1718
    return -p - D;
20,205,570✔
1719
  if (-p + D > l)
80,912,887✔
1720
    return -p + D;
50,089,809✔
1721

1722
  return INFTY;
30,823,078✔
1723
}
1724

1725
double CylindricalMesh::find_phi_crossing(
74,445,558✔
1726
  const Position& r, const Direction& u, double l, int shell) const
1727
{
1728
  // Phi grid is [0, 2Ï€], thus there is no real surface to cross
1729
  if (full_phi_ && (shape_[1] == 1))
74,445,558✔
1730
    return INFTY;
30,474,840✔
1731

1732
  shell = sanitize_phi(shell);
43,970,718✔
1733

1734
  const double p0 = grid_[1][shell];
43,970,718✔
1735

1736
  // solve y(s)/x(s) = tan(p0) = sin(p0)/cos(p0)
1737
  // => x(s) * cos(p0) = y(s) * sin(p0)
1738
  // => (y + s * v) * cos(p0) = (x + s * u) * sin(p0)
1739
  // = s * (v * cos(p0) - u * sin(p0)) = - (y * cos(p0) - x * sin(p0))
1740

1741
  const double c0 = std::cos(p0);
43,970,718✔
1742
  const double s0 = std::sin(p0);
43,970,718✔
1743

1744
  const double denominator = (u.x * s0 - u.y * c0);
43,970,718✔
1745

1746
  // Check if direction of flight is not parallel to phi surface
1747
  if (std::abs(denominator) > FP_PRECISION) {
43,970,718✔
1748
    const double s = -(r.x * s0 - r.y * c0) / denominator;
43,709,974✔
1749
    // Check if solution is in positive direction of flight and crosses the
1750
    // correct phi surface (not -phi)
1751
    if ((s > l) && ((c0 * (r.x + s * u.x) + s0 * (r.y + s * u.y)) > 0.0))
43,709,974✔
1752
      return s;
20,219,859✔
1753
  }
1754

1755
  return INFTY;
23,750,859✔
1756
}
1757

1758
StructuredMesh::MeshDistance CylindricalMesh::find_z_crossing(
36,690,324✔
1759
  const Position& r, const Direction& u, double l, int shell) const
1760
{
1761
  MeshDistance d;
36,690,324✔
1762
  d.next_index = shell;
36,690,324✔
1763

1764
  // Direction of flight is within xy-plane. Will never intersect z.
1765
  if (std::abs(u.z) < FP_PRECISION)
36,690,324✔
1766
    return d;
1,118,216✔
1767

1768
  d.max_surface = (u.z > 0.0);
35,572,108✔
1769
  if (d.max_surface && (shell <= shape_[2])) {
35,572,108✔
1770
    d.next_index += 1;
16,873,241✔
1771
    d.distance = (grid_[2][shell] - r.z) / u.z;
16,873,241✔
1772
  } else if (!d.max_surface && (shell > 0)) {
18,698,867✔
1773
    d.next_index -= 1;
16,843,453✔
1774
    d.distance = (grid_[2][shell - 1] - r.z) / u.z;
16,843,453✔
1775
  }
1776
  return d;
35,572,108✔
1777
}
1778

1779
StructuredMesh::MeshDistance CylindricalMesh::distance_to_grid_boundary(
145,196,535✔
1780
  const MeshIndex& ijk, int i, const Position& r0, const Direction& u,
1781
  double l) const
1782
{
1783
  if (i == 0) {
145,196,535✔
1784

1785
    return std::min(
71,283,432✔
1786
      MeshDistance(ijk[i] + 1, true, find_r_crossing(r0, u, l, ijk[i])),
71,283,432✔
1787
      MeshDistance(ijk[i] - 1, false, find_r_crossing(r0, u, l, ijk[i] - 1)));
142,566,864✔
1788

1789
  } else if (i == 1) {
73,913,103✔
1790

1791
    return std::min(MeshDistance(sanitize_phi(ijk[i] + 1), true,
37,222,779✔
1792
                      find_phi_crossing(r0, u, l, ijk[i])),
37,222,779✔
1793
      MeshDistance(sanitize_phi(ijk[i] - 1), false,
37,222,779✔
1794
        find_phi_crossing(r0, u, l, ijk[i] - 1)));
74,445,558✔
1795

1796
  } else {
1797
    return find_z_crossing(r0, u, l, ijk[i]);
36,690,324✔
1798
  }
1799
}
1800

1801
int CylindricalMesh::set_grid()
412✔
1802
{
1803
  shape_ = {static_cast<int>(grid_[0].size()) - 1,
412✔
1804
    static_cast<int>(grid_[1].size()) - 1,
412✔
1805
    static_cast<int>(grid_[2].size()) - 1};
412✔
1806

1807
  for (const auto& g : grid_) {
1,648✔
1808
    if (g.size() < 2) {
1,236!
UNCOV
1809
      set_errmsg("r-, phi-, and z- grids for cylindrical meshes "
×
1810
                 "must each have at least 2 points");
UNCOV
1811
      return OPENMC_E_INVALID_ARGUMENT;
×
1812
    }
1813
    if (std::adjacent_find(g.begin(), g.end(), std::greater_equal<>()) !=
1,236✔
1814
        g.end()) {
2,472!
UNCOV
1815
      set_errmsg("Values in for r-, phi-, and z- grids for "
×
1816
                 "cylindrical meshes must be sorted and unique.");
UNCOV
1817
      return OPENMC_E_INVALID_ARGUMENT;
×
1818
    }
1819
  }
1820
  if (grid_[0].front() < 0.0) {
412!
UNCOV
1821
    set_errmsg("r-grid for "
×
1822
               "cylindrical meshes must start at r >= 0.");
UNCOV
1823
    return OPENMC_E_INVALID_ARGUMENT;
×
1824
  }
1825
  if (grid_[1].front() < 0.0) {
412!
UNCOV
1826
    set_errmsg("phi-grid for "
×
1827
               "cylindrical meshes must start at phi >= 0.");
UNCOV
1828
    return OPENMC_E_INVALID_ARGUMENT;
×
1829
  }
1830
  if (grid_[1].back() > 2.0 * PI) {
412!
UNCOV
1831
    set_errmsg("phi-grids for "
×
1832
               "cylindrical meshes must end with theta <= 2*pi.");
1833

UNCOV
1834
    return OPENMC_E_INVALID_ARGUMENT;
×
1835
  }
1836

1837
  full_phi_ = (grid_[1].front() == 0.0) && (grid_[1].back() == 2.0 * PI);
412!
1838

1839
  lower_left_ = {origin_[0] - grid_[0].back(), origin_[1] - grid_[0].back(),
824✔
1840
    origin_[2] + grid_[2].front()};
824✔
1841
  upper_right_ = {origin_[0] + grid_[0].back(), origin_[1] + grid_[0].back(),
824✔
1842
    origin_[2] + grid_[2].back()};
824✔
1843

1844
  return 0;
412✔
1845
}
1846

1847
int CylindricalMesh::get_index_in_direction(double r, int i) const
143,180,004✔
1848
{
1849
  return lower_bound_index(grid_[i].begin(), grid_[i].end(), r) + 1;
143,180,004✔
1850
}
1851

UNCOV
1852
std::pair<vector<double>, vector<double>> CylindricalMesh::plot(
×
1853
  Position plot_ll, Position plot_ur) const
1854
{
1855
  fatal_error("Plot of cylindrical Mesh not implemented");
×
1856

1857
  // Figure out which axes lie in the plane of the plot.
1858
  array<vector<double>, 2> axis_lines;
1859
  return {axis_lines[0], axis_lines[1]};
1860
}
1861

1862
void CylindricalMesh::to_hdf5_inner(hid_t mesh_group) const
363✔
1863
{
1864
  write_dataset(mesh_group, "r_grid", grid_[0]);
363✔
1865
  write_dataset(mesh_group, "phi_grid", grid_[1]);
363✔
1866
  write_dataset(mesh_group, "z_grid", grid_[2]);
363✔
1867
  write_dataset(mesh_group, "origin", origin_);
363✔
1868
}
363✔
1869

1870
double CylindricalMesh::volume(const MeshIndex& ijk) const
792✔
1871
{
1872
  double r_i = grid_[0][ijk[0] - 1];
792✔
1873
  double r_o = grid_[0][ijk[0]];
792✔
1874

1875
  double phi_i = grid_[1][ijk[1] - 1];
792✔
1876
  double phi_o = grid_[1][ijk[1]];
792✔
1877

1878
  double z_i = grid_[2][ijk[2] - 1];
792✔
1879
  double z_o = grid_[2][ijk[2]];
792✔
1880

1881
  return 0.5 * (r_o * r_o - r_i * r_i) * (phi_o - phi_i) * (z_o - z_i);
792✔
1882
}
1883

1884
//==============================================================================
1885
// SphericalMesh implementation
1886
//==============================================================================
1887

1888
SphericalMesh::SphericalMesh(pugi::xml_node node)
324✔
1889
  : PeriodicStructuredMesh {node}
324✔
1890
{
1891
  n_dimension_ = 3;
324✔
1892

1893
  grid_[0] = get_node_array<double>(node, "r_grid");
324✔
1894
  grid_[1] = get_node_array<double>(node, "theta_grid");
324✔
1895
  grid_[2] = get_node_array<double>(node, "phi_grid");
324✔
1896
  origin_ = get_node_position(node, "origin");
324✔
1897

1898
  if (int err = set_grid()) {
324!
UNCOV
1899
    fatal_error(openmc_err_msg);
×
1900
  }
1901
}
324✔
1902

NEW
UNCOV
1903
SphericalMesh::SphericalMesh(hid_t group) : PeriodicStructuredMesh {group}
×
1904
{
NEW
UNCOV
1905
  n_dimension_ = 3;
×
1906

NEW
UNCOV
1907
  read_dataset(group, "r_grid", grid_[0]);
×
NEW
UNCOV
1908
  read_dataset(group, "theta_grid", grid_[1]);
×
NEW
UNCOV
1909
  read_dataset(group, "phi_grid", grid_[2]);
×
NEW
UNCOV
1910
  read_dataset(group, "origin", origin_);
×
1911

NEW
UNCOV
1912
  if (int err = set_grid()) {
×
NEW
UNCOV
1913
    fatal_error(openmc_err_msg);
×
1914
  }
NEW
UNCOV
1915
}
×
1916

1917
const std::string SphericalMesh::mesh_type = "spherical";
1918

1919
std::string SphericalMesh::get_mesh_type() const
363✔
1920
{
1921
  return mesh_type;
363✔
1922
}
1923

1924
StructuredMesh::MeshIndex SphericalMesh::get_indices(
68,175,250✔
1925
  Position r, bool& in_mesh) const
1926
{
1927
  r = local_coords(r);
68,175,250✔
1928

1929
  Position mapped_r;
68,175,250✔
1930
  mapped_r[0] = r.norm();
68,175,250✔
1931

1932
  if (mapped_r[0] < FP_PRECISION) {
68,175,250!
UNCOV
1933
    mapped_r[1] = 0.0;
×
UNCOV
1934
    mapped_r[2] = 0.0;
×
1935
  } else {
1936
    mapped_r[1] = std::acos(r.z / mapped_r.x);
68,175,250✔
1937
    mapped_r[2] = std::atan2(r.y, r.x);
68,175,250✔
1938
    if (mapped_r[2] < 0)
68,175,250✔
1939
      mapped_r[2] += 2 * M_PI;
34,062,281✔
1940
  }
1941

1942
  MeshIndex idx = StructuredMesh::get_indices(mapped_r, in_mesh);
68,175,250✔
1943

1944
  idx[1] = sanitize_theta(idx[1]);
68,175,250✔
1945
  idx[2] = sanitize_phi(idx[2]);
68,175,250✔
1946

1947
  return idx;
68,175,250✔
1948
}
1949

1950
Position SphericalMesh::sample_element(
110✔
1951
  const MeshIndex& ijk, uint64_t* seed) const
1952
{
1953
  double r_min = this->r(ijk[0] - 1);
110✔
1954
  double r_max = this->r(ijk[0]);
110✔
1955

1956
  double theta_min = this->theta(ijk[1] - 1);
110✔
1957
  double theta_max = this->theta(ijk[1]);
110✔
1958

1959
  double phi_min = this->phi(ijk[2] - 1);
110✔
1960
  double phi_max = this->phi(ijk[2]);
110✔
1961

1962
  double cos_theta =
1963
    uniform_distribution(std::cos(theta_min), std::cos(theta_max), seed);
110✔
1964
  double sin_theta = std::sin(std::acos(cos_theta));
110✔
1965
  double phi = uniform_distribution(phi_min, phi_max, seed);
110✔
1966
  double r_min_cub = std::pow(r_min, 3);
110✔
1967
  double r_max_cub = std::pow(r_max, 3);
110✔
1968
  // might be faster to do rejection here?
1969
  double r = std::cbrt(uniform_distribution(r_min_cub, r_max_cub, seed));
110✔
1970

1971
  double x = r * std::cos(phi) * sin_theta;
110✔
1972
  double y = r * std::sin(phi) * sin_theta;
110✔
1973
  double z = r * cos_theta;
110✔
1974

1975
  return origin_ + Position(x, y, z);
110✔
1976
}
1977

1978
double SphericalMesh::find_r_crossing(
443,086,688✔
1979
  const Position& r, const Direction& u, double l, int shell) const
1980
{
1981
  if ((shell < 0) || (shell > shape_[0]))
443,086,688✔
1982
    return INFTY;
39,620,317✔
1983

1984
  // solve |r+s*u| = r0
1985
  // |r+s*u| = |r| + 2*s*r*u + s^2 (|u|==1 !)
1986
  const double r0 = grid_[0][shell];
403,466,371✔
1987
  if (r0 == 0.0)
403,466,371✔
1988
    return INFTY;
7,261,639✔
1989
  const double p = r.dot(u);
396,204,732✔
1990
  double c = r.dot(r) - r0 * r0;
396,204,732✔
1991
  double D = p * p - c;
396,204,732✔
1992

1993
  if (std::abs(c) <= RADIAL_MESH_TOL)
396,204,732✔
1994
    return INFTY;
10,598,654✔
1995

1996
  if (D >= 0.0) {
385,606,078✔
1997
    D = std::sqrt(D);
357,729,130✔
1998
    // the solution -p - D is always smaller as -p + D : Check this one first
1999
    if (-p - D > l)
357,729,130✔
2000
      return -p - D;
64,280,876✔
2001
    if (-p + D > l)
293,448,254✔
2002
      return -p + D;
176,905,300✔
2003
  }
2004

2005
  return INFTY;
144,419,902✔
2006
}
2007

2008
double SphericalMesh::find_theta_crossing(
109,327,592✔
2009
  const Position& r, const Direction& u, double l, int shell) const
2010
{
2011
  // Theta grid is [0, π], thus there is no real surface to cross
2012
  if (full_theta_ && (shape_[1] == 1))
109,327,592✔
2013
    return INFTY;
70,969,052✔
2014

2015
  shell = sanitize_theta(shell);
38,358,540✔
2016

2017
  // solving z(s) = cos/theta) * r(s) with r(s) = r+s*u
2018
  // yields
2019
  // a*s^2 + 2*b*s + c == 0 with
2020
  // a = cos(theta)^2 - u.z * u.z
2021
  // b = r*u * cos(theta)^2 - u.z * r.z
2022
  // c = r*r * cos(theta)^2 - r.z^2
2023

2024
  const double cos_t = std::cos(grid_[1][shell]);
38,358,540✔
2025
  const bool sgn = std::signbit(cos_t);
38,358,540✔
2026
  const double cos_t_2 = cos_t * cos_t;
38,358,540✔
2027

2028
  const double a = cos_t_2 - u.z * u.z;
38,358,540✔
2029
  const double b = r.dot(u) * cos_t_2 - r.z * u.z;
38,358,540✔
2030
  const double c = r.dot(r) * cos_t_2 - r.z * r.z;
38,358,540✔
2031

2032
  // if factor of s^2 is zero, direction of flight is parallel to theta
2033
  // surface
2034
  if (std::abs(a) < FP_PRECISION) {
38,358,540✔
2035
    // if b vanishes, direction of flight is within theta surface and crossing
2036
    // is not possible
2037
    if (std::abs(b) < FP_PRECISION)
482,548!
2038
      return INFTY;
482,548✔
2039

UNCOV
2040
    const double s = -0.5 * c / b;
×
2041
    // Check if solution is in positive direction of flight and has correct
2042
    // sign
UNCOV
2043
    if ((s > l) && (std::signbit(r.z + s * u.z) == sgn))
×
UNCOV
2044
      return s;
×
2045

2046
    // no crossing is possible
UNCOV
2047
    return INFTY;
×
2048
  }
2049

2050
  const double p = b / a;
37,875,992✔
2051
  double D = p * p - c / a;
37,875,992✔
2052

2053
  if (D < 0.0)
37,875,992✔
2054
    return INFTY;
10,954,988✔
2055

2056
  D = std::sqrt(D);
26,921,004✔
2057

2058
  // the solution -p-D is always smaller as -p+D : Check this one first
2059
  double s = -p - D;
26,921,004✔
2060
  // Check if solution is in positive direction of flight and has correct sign
2061
  if ((s > l) && (std::signbit(r.z + s * u.z) == sgn))
26,921,004✔
2062
    return s;
5,282,607✔
2063

2064
  s = -p + D;
21,638,397✔
2065
  // Check if solution is in positive direction of flight and has correct sign
2066
  if ((s > l) && (std::signbit(r.z + s * u.z) == sgn))
21,638,397✔
2067
    return s;
10,163,296✔
2068

2069
  return INFTY;
11,475,101✔
2070
}
2071

2072
double SphericalMesh::find_phi_crossing(
110,917,070✔
2073
  const Position& r, const Direction& u, double l, int shell) const
2074
{
2075
  // Phi grid is [0, 2Ï€], thus there is no real surface to cross
2076
  if (full_phi_ && (shape_[2] == 1))
110,917,070✔
2077
    return INFTY;
70,969,052✔
2078

2079
  shell = sanitize_phi(shell);
39,948,018✔
2080

2081
  const double p0 = grid_[2][shell];
39,948,018✔
2082

2083
  // solve y(s)/x(s) = tan(p0) = sin(p0)/cos(p0)
2084
  // => x(s) * cos(p0) = y(s) * sin(p0)
2085
  // => (y + s * v) * cos(p0) = (x + s * u) * sin(p0)
2086
  // = s * (v * cos(p0) - u * sin(p0)) = - (y * cos(p0) - x * sin(p0))
2087

2088
  const double c0 = std::cos(p0);
39,948,018✔
2089
  const double s0 = std::sin(p0);
39,948,018✔
2090

2091
  const double denominator = (u.x * s0 - u.y * c0);
39,948,018✔
2092

2093
  // Check if direction of flight is not parallel to phi surface
2094
  if (std::abs(denominator) > FP_PRECISION) {
39,948,018✔
2095
    const double s = -(r.x * s0 - r.y * c0) / denominator;
39,714,026✔
2096
    // Check if solution is in positive direction of flight and crosses the
2097
    // correct phi surface (not -phi)
2098
    if ((s > l) && ((c0 * (r.x + s * u.x) + s0 * (r.y + s * u.y)) > 0.0))
39,714,026✔
2099
      return s;
17,579,452✔
2100
  }
2101

2102
  return INFTY;
22,368,566✔
2103
}
2104

2105
StructuredMesh::MeshDistance SphericalMesh::distance_to_grid_boundary(
331,665,675✔
2106
  const MeshIndex& ijk, int i, const Position& r0, const Direction& u,
2107
  double l) const
2108
{
2109

2110
  if (i == 0) {
331,665,675✔
2111
    return std::min(
221,543,344✔
2112
      MeshDistance(ijk[i] + 1, true, find_r_crossing(r0, u, l, ijk[i])),
221,543,344✔
2113
      MeshDistance(ijk[i] - 1, false, find_r_crossing(r0, u, l, ijk[i] - 1)));
443,086,688✔
2114

2115
  } else if (i == 1) {
110,122,331✔
2116
    return std::min(MeshDistance(sanitize_theta(ijk[i] + 1), true,
54,663,796✔
2117
                      find_theta_crossing(r0, u, l, ijk[i])),
54,663,796✔
2118
      MeshDistance(sanitize_theta(ijk[i] - 1), false,
54,663,796✔
2119
        find_theta_crossing(r0, u, l, ijk[i] - 1)));
109,327,592✔
2120

2121
  } else {
2122
    return std::min(MeshDistance(sanitize_phi(ijk[i] + 1), true,
55,458,535✔
2123
                      find_phi_crossing(r0, u, l, ijk[i])),
55,458,535✔
2124
      MeshDistance(sanitize_phi(ijk[i] - 1), false,
55,458,535✔
2125
        find_phi_crossing(r0, u, l, ijk[i] - 1)));
110,917,070✔
2126
  }
2127
}
2128

2129
int SphericalMesh::set_grid()
346✔
2130
{
2131
  shape_ = {static_cast<int>(grid_[0].size()) - 1,
346✔
2132
    static_cast<int>(grid_[1].size()) - 1,
346✔
2133
    static_cast<int>(grid_[2].size()) - 1};
346✔
2134

2135
  for (const auto& g : grid_) {
1,384✔
2136
    if (g.size() < 2) {
1,038!
UNCOV
2137
      set_errmsg("x-, y-, and z- grids for spherical meshes "
×
2138
                 "must each have at least 2 points");
UNCOV
2139
      return OPENMC_E_INVALID_ARGUMENT;
×
2140
    }
2141
    if (std::adjacent_find(g.begin(), g.end(), std::greater_equal<>()) !=
1,038✔
2142
        g.end()) {
2,076!
UNCOV
2143
      set_errmsg("Values in for r-, theta-, and phi- grids for "
×
2144
                 "spherical meshes must be sorted and unique.");
UNCOV
2145
      return OPENMC_E_INVALID_ARGUMENT;
×
2146
    }
2147
    if (g.front() < 0.0) {
1,038!
UNCOV
2148
      set_errmsg("r-, theta-, and phi- grids for "
×
2149
                 "spherical meshes must start at v >= 0.");
UNCOV
2150
      return OPENMC_E_INVALID_ARGUMENT;
×
2151
    }
2152
  }
2153
  if (grid_[1].back() > PI) {
346!
UNCOV
2154
    set_errmsg("theta-grids for "
×
2155
               "spherical meshes must end with theta <= pi.");
2156

UNCOV
2157
    return OPENMC_E_INVALID_ARGUMENT;
×
2158
  }
2159
  if (grid_[2].back() > 2 * PI) {
346!
UNCOV
2160
    set_errmsg("phi-grids for "
×
2161
               "spherical meshes must end with phi <= 2*pi.");
UNCOV
2162
    return OPENMC_E_INVALID_ARGUMENT;
×
2163
  }
2164

2165
  full_theta_ = (grid_[1].front() == 0.0) && (grid_[1].back() == PI);
346!
2166
  full_phi_ = (grid_[2].front() == 0.0) && (grid_[2].back() == 2 * PI);
346✔
2167

2168
  double r = grid_[0].back();
346✔
2169
  lower_left_ = {origin_[0] - r, origin_[1] - r, origin_[2] - r};
346✔
2170
  upper_right_ = {origin_[0] + r, origin_[1] + r, origin_[2] + r};
346✔
2171

2172
  return 0;
346✔
2173
}
2174

2175
int SphericalMesh::get_index_in_direction(double r, int i) const
204,525,750✔
2176
{
2177
  return lower_bound_index(grid_[i].begin(), grid_[i].end(), r) + 1;
204,525,750✔
2178
}
2179

2180
std::pair<vector<double>, vector<double>> SphericalMesh::plot(
×
2181
  Position plot_ll, Position plot_ur) const
2182
{
UNCOV
2183
  fatal_error("Plot of spherical Mesh not implemented");
×
2184

2185
  // Figure out which axes lie in the plane of the plot.
2186
  array<vector<double>, 2> axis_lines;
2187
  return {axis_lines[0], axis_lines[1]};
2188
}
2189

2190
void SphericalMesh::to_hdf5_inner(hid_t mesh_group) const
297✔
2191
{
2192
  write_dataset(mesh_group, "r_grid", grid_[0]);
297✔
2193
  write_dataset(mesh_group, "theta_grid", grid_[1]);
297✔
2194
  write_dataset(mesh_group, "phi_grid", grid_[2]);
297✔
2195
  write_dataset(mesh_group, "origin", origin_);
297✔
2196
}
297✔
2197

2198
double SphericalMesh::volume(const MeshIndex& ijk) const
935✔
2199
{
2200
  double r_i = grid_[0][ijk[0] - 1];
935✔
2201
  double r_o = grid_[0][ijk[0]];
935✔
2202

2203
  double theta_i = grid_[1][ijk[1] - 1];
935✔
2204
  double theta_o = grid_[1][ijk[1]];
935✔
2205

2206
  double phi_i = grid_[2][ijk[2] - 1];
935✔
2207
  double phi_o = grid_[2][ijk[2]];
935✔
2208

2209
  return (1.0 / 3.0) * (r_o * r_o * r_o - r_i * r_i * r_i) *
935✔
2210
         (std::cos(theta_i) - std::cos(theta_o)) * (phi_o - phi_i);
935✔
2211
}
2212

2213
//==============================================================================
2214
// Helper functions for the C API
2215
//==============================================================================
2216

2217
int check_mesh(int32_t index)
6,259✔
2218
{
2219
  if (index < 0 || index >= model::meshes.size()) {
6,259!
UNCOV
2220
    set_errmsg("Index in meshes array is out of bounds.");
×
UNCOV
2221
    return OPENMC_E_OUT_OF_BOUNDS;
×
2222
  }
2223
  return 0;
6,259✔
2224
}
2225

2226
template<class T>
2227
int check_mesh_type(int32_t index)
1,100✔
2228
{
2229
  if (int err = check_mesh(index))
1,100!
UNCOV
2230
    return err;
×
2231

2232
  T* mesh = dynamic_cast<T*>(model::meshes[index].get());
1,100!
2233
  if (!mesh) {
1,100!
UNCOV
2234
    set_errmsg("This function is not valid for input mesh.");
×
UNCOV
2235
    return OPENMC_E_INVALID_TYPE;
×
2236
  }
2237
  return 0;
1,100✔
2238
}
2239

2240
template<class T>
2241
bool is_mesh_type(int32_t index)
2242
{
2243
  T* mesh = dynamic_cast<T*>(model::meshes[index].get());
2244
  return mesh;
2245
}
2246

2247
//==============================================================================
2248
// C API functions
2249
//==============================================================================
2250

2251
// Return the type of mesh as a C string
2252
extern "C" int openmc_mesh_get_type(int32_t index, char* type)
1,441✔
2253
{
2254
  if (int err = check_mesh(index))
1,441!
UNCOV
2255
    return err;
×
2256

2257
  std::strcpy(type, model::meshes[index].get()->get_mesh_type().c_str());
1,441✔
2258

2259
  return 0;
1,441✔
2260
}
2261

2262
//! Extend the meshes array by n elements
2263
extern "C" int openmc_extend_meshes(
253✔
2264
  int32_t n, const char* type, int32_t* index_start, int32_t* index_end)
2265
{
2266
  if (index_start)
253!
2267
    *index_start = model::meshes.size();
253✔
2268
  std::string mesh_type;
253✔
2269

2270
  for (int i = 0; i < n; ++i) {
506✔
2271
    if (RegularMesh::mesh_type == type) {
253✔
2272
      model::meshes.push_back(make_unique<RegularMesh>());
165✔
2273
    } else if (RectilinearMesh::mesh_type == type) {
88✔
2274
      model::meshes.push_back(make_unique<RectilinearMesh>());
44✔
2275
    } else if (CylindricalMesh::mesh_type == type) {
44✔
2276
      model::meshes.push_back(make_unique<CylindricalMesh>());
22✔
2277
    } else if (SphericalMesh::mesh_type == type) {
22!
2278
      model::meshes.push_back(make_unique<SphericalMesh>());
22✔
2279
    } else {
UNCOV
2280
      throw std::runtime_error {"Unknown mesh type: " + std::string(type)};
×
2281
    }
2282
  }
2283
  if (index_end)
253!
UNCOV
2284
    *index_end = model::meshes.size() - 1;
×
2285

2286
  return 0;
253✔
2287
}
253✔
2288

2289
//! Adds a new unstructured mesh to OpenMC
UNCOV
2290
extern "C" int openmc_add_unstructured_mesh(
×
2291
  const char filename[], const char library[], int* id)
2292
{
UNCOV
2293
  std::string lib_name(library);
×
UNCOV
2294
  std::string mesh_file(filename);
×
UNCOV
2295
  bool valid_lib = false;
×
2296

2297
#ifdef OPENMC_DAGMC_ENABLED
2298
  if (lib_name == MOABMesh::mesh_lib_type) {
×
2299
    model::meshes.push_back(std::move(make_unique<MOABMesh>(mesh_file)));
×
2300
    valid_lib = true;
2301
  }
2302
#endif
2303

2304
#ifdef OPENMC_LIBMESH_ENABLED
2305
  if (lib_name == LibMesh::mesh_lib_type) {
×
2306
    model::meshes.push_back(std::move(make_unique<LibMesh>(mesh_file)));
×
2307
    valid_lib = true;
2308
  }
2309
#endif
2310

UNCOV
2311
  if (!valid_lib) {
×
2312
    set_errmsg(fmt::format("Mesh library {} is not supported "
×
2313
                           "by this build of OpenMC",
2314
      lib_name));
UNCOV
2315
    return OPENMC_E_INVALID_ARGUMENT;
×
2316
  }
2317

2318
  // auto-assign new ID
UNCOV
2319
  model::meshes.back()->set_id(-1);
×
UNCOV
2320
  *id = model::meshes.back()->id_;
×
2321

2322
  return 0;
×
UNCOV
2323
}
×
2324

2325
//! Return the index in the meshes array of a mesh with a given ID
2326
extern "C" int openmc_get_mesh_index(int32_t id, int32_t* index)
429✔
2327
{
2328
  auto pair = model::mesh_map.find(id);
429✔
2329
  if (pair == model::mesh_map.end()) {
429!
UNCOV
2330
    set_errmsg("No mesh exists with ID=" + std::to_string(id) + ".");
×
UNCOV
2331
    return OPENMC_E_INVALID_ID;
×
2332
  }
2333
  *index = pair->second;
429✔
2334
  return 0;
429✔
2335
}
2336

2337
//! Return the ID of a mesh
2338
extern "C" int openmc_mesh_get_id(int32_t index, int32_t* id)
2,772✔
2339
{
2340
  if (int err = check_mesh(index))
2,772!
UNCOV
2341
    return err;
×
2342
  *id = model::meshes[index]->id_;
2,772✔
2343
  return 0;
2,772✔
2344
}
2345

2346
//! Set the ID of a mesh
2347
extern "C" int openmc_mesh_set_id(int32_t index, int32_t id)
253✔
2348
{
2349
  if (int err = check_mesh(index))
253!
UNCOV
2350
    return err;
×
2351
  model::meshes[index]->id_ = id;
253✔
2352
  model::mesh_map[id] = index;
253✔
2353
  return 0;
253✔
2354
}
2355

2356
//! Get the number of elements in a mesh
2357
extern "C" int openmc_mesh_get_n_elements(int32_t index, size_t* n)
253✔
2358
{
2359
  if (int err = check_mesh(index))
253!
UNCOV
2360
    return err;
×
2361
  *n = model::meshes[index]->n_bins();
253✔
2362
  return 0;
253✔
2363
}
2364

2365
//! Get the volume of each element in the mesh
2366
extern "C" int openmc_mesh_get_volumes(int32_t index, double* volumes)
88✔
2367
{
2368
  if (int err = check_mesh(index))
88!
UNCOV
2369
    return err;
×
2370
  for (int i = 0; i < model::meshes[index]->n_bins(); ++i) {
968✔
2371
    volumes[i] = model::meshes[index]->volume(i);
880✔
2372
  }
2373
  return 0;
88✔
2374
}
2375

2376
//! Get the bounding box of a mesh
2377
extern "C" int openmc_mesh_bounding_box(int32_t index, double* ll, double* ur)
143✔
2378
{
2379
  if (int err = check_mesh(index))
143!
UNCOV
2380
    return err;
×
2381

2382
  BoundingBox bbox = model::meshes[index]->bounding_box();
143✔
2383

2384
  // set lower left corner values
2385
  ll[0] = bbox.xmin;
143✔
2386
  ll[1] = bbox.ymin;
143✔
2387
  ll[2] = bbox.zmin;
143✔
2388

2389
  // set upper right corner values
2390
  ur[0] = bbox.xmax;
143✔
2391
  ur[1] = bbox.ymax;
143✔
2392
  ur[2] = bbox.zmax;
143✔
2393
  return 0;
143✔
2394
}
2395

2396
extern "C" int openmc_mesh_material_volumes(int32_t index, int nx, int ny,
165✔
2397
  int nz, int table_size, int32_t* materials, double* volumes)
2398
{
2399
  if (int err = check_mesh(index))
165!
UNCOV
2400
    return err;
×
2401

2402
  try {
2403
    model::meshes[index]->material_volumes(
165✔
2404
      nx, ny, nz, table_size, materials, volumes);
2405
  } catch (const std::exception& e) {
11!
2406
    set_errmsg(e.what());
11✔
2407
    if (starts_with(e.what(), "Mesh")) {
11!
2408
      return OPENMC_E_GEOMETRY;
11✔
2409
    } else {
UNCOV
2410
      return OPENMC_E_ALLOCATE;
×
2411
    }
2412
  }
11✔
2413

2414
  return 0;
154✔
2415
}
2416

2417
extern "C" int openmc_mesh_get_plot_bins(int32_t index, Position origin,
44✔
2418
  Position width, int basis, int* pixels, int32_t* data)
2419
{
2420
  if (int err = check_mesh(index))
44!
UNCOV
2421
    return err;
×
2422
  const auto& mesh = model::meshes[index].get();
44✔
2423

2424
  int pixel_width = pixels[0];
44✔
2425
  int pixel_height = pixels[1];
44✔
2426

2427
  // get pixel size
2428
  double in_pixel = (width[0]) / static_cast<double>(pixel_width);
44✔
2429
  double out_pixel = (width[1]) / static_cast<double>(pixel_height);
44✔
2430

2431
  // setup basis indices and initial position centered on pixel
2432
  int in_i, out_i;
2433
  Position xyz = origin;
44✔
2434
  enum class PlotBasis { xy = 1, xz = 2, yz = 3 };
2435
  PlotBasis basis_enum = static_cast<PlotBasis>(basis);
44✔
2436
  switch (basis_enum) {
44!
2437
  case PlotBasis::xy:
44✔
2438
    in_i = 0;
44✔
2439
    out_i = 1;
44✔
2440
    break;
44✔
UNCOV
2441
  case PlotBasis::xz:
×
2442
    in_i = 0;
×
UNCOV
2443
    out_i = 2;
×
UNCOV
2444
    break;
×
UNCOV
2445
  case PlotBasis::yz:
×
UNCOV
2446
    in_i = 1;
×
UNCOV
2447
    out_i = 2;
×
UNCOV
2448
    break;
×
UNCOV
2449
  default:
×
UNCOV
2450
    UNREACHABLE();
×
2451
  }
2452

2453
  // set initial position
2454
  xyz[in_i] = origin[in_i] - width[0] / 2. + in_pixel / 2.;
44✔
2455
  xyz[out_i] = origin[out_i] + width[1] / 2. - out_pixel / 2.;
44✔
2456

2457
#pragma omp parallel
24✔
2458
  {
2459
    Position r = xyz;
20✔
2460

2461
#pragma omp for
2462
    for (int y = 0; y < pixel_height; y++) {
420✔
2463
      r[out_i] = xyz[out_i] - out_pixel * y;
400✔
2464
      for (int x = 0; x < pixel_width; x++) {
8,400✔
2465
        r[in_i] = xyz[in_i] + in_pixel * x;
8,000✔
2466
        data[pixel_width * y + x] = mesh->get_bin(r);
8,000✔
2467
      }
2468
    }
2469
  }
2470

2471
  return 0;
44✔
2472
}
2473

2474
//! Get the dimension of a regular mesh
2475
extern "C" int openmc_regular_mesh_get_dimension(
11✔
2476
  int32_t index, int** dims, int* n)
2477
{
2478
  if (int err = check_mesh_type<RegularMesh>(index))
11!
2479
    return err;
×
2480
  RegularMesh* mesh = dynamic_cast<RegularMesh*>(model::meshes[index].get());
11!
2481
  *dims = mesh->shape_.data();
11✔
2482
  *n = mesh->n_dimension_;
11✔
2483
  return 0;
11✔
2484
}
2485

2486
//! Set the dimension of a regular mesh
2487
extern "C" int openmc_regular_mesh_set_dimension(
187✔
2488
  int32_t index, int n, const int* dims)
2489
{
2490
  if (int err = check_mesh_type<RegularMesh>(index))
187!
UNCOV
2491
    return err;
×
2492
  RegularMesh* mesh = dynamic_cast<RegularMesh*>(model::meshes[index].get());
187!
2493

2494
  // Copy dimension
2495
  mesh->n_dimension_ = n;
187✔
2496
  std::copy(dims, dims + n, mesh->shape_.begin());
187✔
2497
  return 0;
187✔
2498
}
2499

2500
//! Get the regular mesh parameters
2501
extern "C" int openmc_regular_mesh_get_params(
209✔
2502
  int32_t index, double** ll, double** ur, double** width, int* n)
2503
{
2504
  if (int err = check_mesh_type<RegularMesh>(index))
209!
UNCOV
2505
    return err;
×
2506
  RegularMesh* m = dynamic_cast<RegularMesh*>(model::meshes[index].get());
209!
2507

2508
  if (m->lower_left_.dimension() == 0) {
209!
UNCOV
2509
    set_errmsg("Mesh parameters have not been set.");
×
UNCOV
2510
    return OPENMC_E_ALLOCATE;
×
2511
  }
2512

2513
  *ll = m->lower_left_.data();
209✔
2514
  *ur = m->upper_right_.data();
209✔
2515
  *width = m->width_.data();
209✔
2516
  *n = m->n_dimension_;
209✔
2517
  return 0;
209✔
2518
}
2519

2520
//! Set the regular mesh parameters
2521
extern "C" int openmc_regular_mesh_set_params(
220✔
2522
  int32_t index, int n, const double* ll, const double* ur, const double* width)
2523
{
2524
  if (int err = check_mesh_type<RegularMesh>(index))
220!
UNCOV
2525
    return err;
×
2526
  RegularMesh* m = dynamic_cast<RegularMesh*>(model::meshes[index].get());
220!
2527

2528
  if (m->n_dimension_ == -1) {
220!
UNCOV
2529
    set_errmsg("Need to set mesh dimension before setting parameters.");
×
UNCOV
2530
    return OPENMC_E_UNASSIGNED;
×
2531
  }
2532

2533
  vector<std::size_t> shape = {static_cast<std::size_t>(n)};
220✔
2534
  if (ll && ur) {
220✔
2535
    m->lower_left_ = xt::adapt(ll, n, xt::no_ownership(), shape);
198✔
2536
    m->upper_right_ = xt::adapt(ur, n, xt::no_ownership(), shape);
198✔
2537
    m->width_ = (m->upper_right_ - m->lower_left_) / m->get_x_shape();
198✔
2538
  } else if (ll && width) {
22!
2539
    m->lower_left_ = xt::adapt(ll, n, xt::no_ownership(), shape);
11✔
2540
    m->width_ = xt::adapt(width, n, xt::no_ownership(), shape);
11✔
2541
    m->upper_right_ = m->lower_left_ + m->get_x_shape() * m->width_;
11✔
2542
  } else if (ur && width) {
11!
2543
    m->upper_right_ = xt::adapt(ur, n, xt::no_ownership(), shape);
11✔
2544
    m->width_ = xt::adapt(width, n, xt::no_ownership(), shape);
11✔
2545
    m->lower_left_ = m->upper_right_ - m->get_x_shape() * m->width_;
11✔
2546
  } else {
UNCOV
2547
    set_errmsg("At least two parameters must be specified.");
×
UNCOV
2548
    return OPENMC_E_INVALID_ARGUMENT;
×
2549
  }
2550

2551
  // Set material volumes
2552

2553
  // TODO: incorporate this into method in RegularMesh that can be called from
2554
  // here and from constructor
2555
  m->volume_frac_ = 1.0 / xt::prod(m->get_x_shape())();
220✔
2556
  m->element_volume_ = 1.0;
220✔
2557
  for (int i = 0; i < m->n_dimension_; i++) {
880✔
2558
    m->element_volume_ *= m->width_[i];
660✔
2559
  }
2560

2561
  return 0;
220✔
2562
}
220✔
2563

2564
//! Set the mesh parameters for rectilinear, cylindrical and spharical meshes
2565
template<class C>
2566
int openmc_structured_mesh_set_grid_impl(int32_t index, const double* grid_x,
88✔
2567
  const int nx, const double* grid_y, const int ny, const double* grid_z,
2568
  const int nz)
2569
{
2570
  if (int err = check_mesh_type<C>(index))
88!
UNCOV
2571
    return err;
×
2572

2573
  C* m = dynamic_cast<C*>(model::meshes[index].get());
88!
2574

2575
  m->n_dimension_ = 3;
88✔
2576

2577
  m->grid_[0].reserve(nx);
88✔
2578
  m->grid_[1].reserve(ny);
88✔
2579
  m->grid_[2].reserve(nz);
88✔
2580

2581
  for (int i = 0; i < nx; i++) {
572✔
2582
    m->grid_[0].push_back(grid_x[i]);
484✔
2583
  }
2584
  for (int i = 0; i < ny; i++) {
341✔
2585
    m->grid_[1].push_back(grid_y[i]);
253✔
2586
  }
2587
  for (int i = 0; i < nz; i++) {
319✔
2588
    m->grid_[2].push_back(grid_z[i]);
231✔
2589
  }
2590

2591
  int err = m->set_grid();
88✔
2592
  return err;
88✔
2593
}
2594

2595
//! Get the mesh parameters for rectilinear, cylindrical and spherical meshes
2596
template<class C>
2597
int openmc_structured_mesh_get_grid_impl(int32_t index, double** grid_x,
385✔
2598
  int* nx, double** grid_y, int* ny, double** grid_z, int* nz)
2599
{
2600
  if (int err = check_mesh_type<C>(index))
385!
UNCOV
2601
    return err;
×
2602
  C* m = dynamic_cast<C*>(model::meshes[index].get());
385!
2603

2604
  if (m->lower_left_.dimension() == 0) {
385!
UNCOV
2605
    set_errmsg("Mesh parameters have not been set.");
×
UNCOV
2606
    return OPENMC_E_ALLOCATE;
×
2607
  }
2608

2609
  *grid_x = m->grid_[0].data();
385✔
2610
  *nx = m->grid_[0].size();
385✔
2611
  *grid_y = m->grid_[1].data();
385✔
2612
  *ny = m->grid_[1].size();
385✔
2613
  *grid_z = m->grid_[2].data();
385✔
2614
  *nz = m->grid_[2].size();
385✔
2615

2616
  return 0;
385✔
2617
}
2618

2619
//! Get the rectilinear mesh grid
2620
extern "C" int openmc_rectilinear_mesh_get_grid(int32_t index, double** grid_x,
143✔
2621
  int* nx, double** grid_y, int* ny, double** grid_z, int* nz)
2622
{
2623
  return openmc_structured_mesh_get_grid_impl<RectilinearMesh>(
143✔
2624
    index, grid_x, nx, grid_y, ny, grid_z, nz);
143✔
2625
}
2626

2627
//! Set the rectilienar mesh parameters
2628
extern "C" int openmc_rectilinear_mesh_set_grid(int32_t index,
44✔
2629
  const double* grid_x, const int nx, const double* grid_y, const int ny,
2630
  const double* grid_z, const int nz)
2631
{
2632
  return openmc_structured_mesh_set_grid_impl<RectilinearMesh>(
44✔
2633
    index, grid_x, nx, grid_y, ny, grid_z, nz);
44✔
2634
}
2635

2636
//! Get the cylindrical mesh grid
2637
extern "C" int openmc_cylindrical_mesh_get_grid(int32_t index, double** grid_x,
121✔
2638
  int* nx, double** grid_y, int* ny, double** grid_z, int* nz)
2639
{
2640
  return openmc_structured_mesh_get_grid_impl<CylindricalMesh>(
121✔
2641
    index, grid_x, nx, grid_y, ny, grid_z, nz);
121✔
2642
}
2643

2644
//! Set the cylindrical mesh parameters
2645
extern "C" int openmc_cylindrical_mesh_set_grid(int32_t index,
22✔
2646
  const double* grid_x, const int nx, const double* grid_y, const int ny,
2647
  const double* grid_z, const int nz)
2648
{
2649
  return openmc_structured_mesh_set_grid_impl<CylindricalMesh>(
22✔
2650
    index, grid_x, nx, grid_y, ny, grid_z, nz);
22✔
2651
}
2652

2653
//! Get the spherical mesh grid
2654
extern "C" int openmc_spherical_mesh_get_grid(int32_t index, double** grid_x,
121✔
2655
  int* nx, double** grid_y, int* ny, double** grid_z, int* nz)
2656
{
2657

2658
  return openmc_structured_mesh_get_grid_impl<SphericalMesh>(
121✔
2659
    index, grid_x, nx, grid_y, ny, grid_z, nz);
121✔
2660
  ;
2661
}
2662

2663
//! Set the spherical mesh parameters
2664
extern "C" int openmc_spherical_mesh_set_grid(int32_t index,
22✔
2665
  const double* grid_x, const int nx, const double* grid_y, const int ny,
2666
  const double* grid_z, const int nz)
2667
{
2668
  return openmc_structured_mesh_set_grid_impl<SphericalMesh>(
22✔
2669
    index, grid_x, nx, grid_y, ny, grid_z, nz);
22✔
2670
}
2671

2672
#ifdef OPENMC_DAGMC_ENABLED
2673

2674
const std::string MOABMesh::mesh_lib_type = "moab";
2675

2676
MOABMesh::MOABMesh(pugi::xml_node node) : UnstructuredMesh(node)
23✔
2677
{
2678
  initialize();
23✔
2679
}
23✔
2680

2681
MOABMesh::MOABMesh(hid_t group) : UnstructuredMesh(group)
×
2682
{
2683
  initialize();
×
2684
}
2685

2686
MOABMesh::MOABMesh(const std::string& filename, double length_multiplier)
2687
  : UnstructuredMesh()
×
2688
{
2689
  n_dimension_ = 3;
2690
  filename_ = filename;
×
2691
  set_length_multiplier(length_multiplier);
×
2692
  initialize();
×
2693
}
2694

2695
MOABMesh::MOABMesh(std::shared_ptr<moab::Interface> external_mbi)
1✔
2696
{
2697
  mbi_ = external_mbi;
1✔
2698
  filename_ = "unknown (external file)";
1✔
2699
  this->initialize();
1✔
2700
}
1✔
2701

2702
void MOABMesh::initialize()
24✔
2703
{
2704

2705
  // Create the MOAB interface and load data from file
2706
  this->create_interface();
24✔
2707

2708
  // Initialise MOAB error code
2709
  moab::ErrorCode rval = moab::MB_SUCCESS;
24✔
2710

2711
  // Set the dimension
2712
  n_dimension_ = 3;
24✔
2713

2714
  // set member range of tetrahedral entities
2715
  rval = mbi_->get_entities_by_dimension(0, n_dimension_, ehs_);
24✔
2716
  if (rval != moab::MB_SUCCESS) {
24!
2717
    fatal_error("Failed to get all tetrahedral elements");
2718
  }
2719

2720
  if (!ehs_.all_of_type(moab::MBTET)) {
24!
2721
    warning("Non-tetrahedral elements found in unstructured "
×
2722
            "mesh file: " +
2723
            filename_);
2724
  }
2725

2726
  // set member range of vertices
2727
  int vertex_dim = 0;
24✔
2728
  rval = mbi_->get_entities_by_dimension(0, vertex_dim, verts_);
24✔
2729
  if (rval != moab::MB_SUCCESS) {
24!
2730
    fatal_error("Failed to get all vertex handles");
2731
  }
2732

2733
  // make an entity set for all tetrahedra
2734
  // this is used for convenience later in output
2735
  rval = mbi_->create_meshset(moab::MESHSET_SET, tetset_);
24✔
2736
  if (rval != moab::MB_SUCCESS) {
24!
2737
    fatal_error("Failed to create an entity set for the tetrahedral elements");
2738
  }
2739

2740
  rval = mbi_->add_entities(tetset_, ehs_);
24✔
2741
  if (rval != moab::MB_SUCCESS) {
24!
2742
    fatal_error("Failed to add tetrahedra to an entity set.");
2743
  }
2744

2745
  if (length_multiplier_ > 0.0) {
24!
2746
    // get the connectivity of all tets
2747
    moab::Range adj;
×
2748
    rval = mbi_->get_adjacencies(ehs_, 0, true, adj, moab::Interface::UNION);
×
2749
    if (rval != moab::MB_SUCCESS) {
×
2750
      fatal_error("Failed to get adjacent vertices of tetrahedra.");
2751
    }
2752
    // scale all vertex coords by multiplier (done individually so not all
2753
    // coordinates are in memory twice at once)
2754
    for (auto vert : adj) {
×
2755
      // retrieve coords
2756
      std::array<double, 3> coord;
2757
      rval = mbi_->get_coords(&vert, 1, coord.data());
×
2758
      if (rval != moab::MB_SUCCESS) {
×
2759
        fatal_error("Could not get coordinates of vertex.");
2760
      }
2761
      // scale coords
2762
      for (auto& c : coord) {
×
2763
        c *= length_multiplier_;
2764
      }
2765
      // set new coords
2766
      rval = mbi_->set_coords(&vert, 1, coord.data());
×
2767
      if (rval != moab::MB_SUCCESS) {
×
2768
        fatal_error("Failed to set new vertex coordinates");
2769
      }
2770
    }
2771
  }
2772

2773
  // Determine bounds of mesh
2774
  this->determine_bounds();
24✔
2775
}
24✔
2776

2777
void MOABMesh::prepare_for_point_location()
20✔
2778
{
2779
  // if the KDTree has already been constructed, do nothing
2780
  if (kdtree_)
20!
2781
    return;
2782

2783
  // build acceleration data structures
2784
  compute_barycentric_data(ehs_);
20✔
2785
  build_kdtree(ehs_);
20✔
2786
}
2787

2788
void MOABMesh::create_interface()
24✔
2789
{
2790
  // Do not create a MOAB instance if one is already in memory
2791
  if (mbi_)
24✔
2792
    return;
1✔
2793

2794
  // create MOAB instance
2795
  mbi_ = std::make_shared<moab::Core>();
23✔
2796

2797
  // load unstructured mesh file
2798
  moab::ErrorCode rval = mbi_->load_file(filename_.c_str());
23✔
2799
  if (rval != moab::MB_SUCCESS) {
23!
2800
    fatal_error("Failed to load the unstructured mesh file: " + filename_);
2801
  }
2802
}
2803

2804
void MOABMesh::build_kdtree(const moab::Range& all_tets)
20✔
2805
{
2806
  moab::Range all_tris;
20✔
2807
  int adj_dim = 2;
20✔
2808
  write_message("Getting tet adjacencies...", 7);
20✔
2809
  moab::ErrorCode rval = mbi_->get_adjacencies(
20✔
2810
    all_tets, adj_dim, true, all_tris, moab::Interface::UNION);
2811
  if (rval != moab::MB_SUCCESS) {
20!
2812
    fatal_error("Failed to get adjacent triangles for tets");
2813
  }
2814

2815
  if (!all_tris.all_of_type(moab::MBTRI)) {
20!
2816
    warning("Non-triangle elements found in tet adjacencies in "
×
2817
            "unstructured mesh file: " +
2818
            filename_);
×
2819
  }
2820

2821
  // combine into one range
2822
  moab::Range all_tets_and_tris;
20✔
2823
  all_tets_and_tris.merge(all_tets);
20✔
2824
  all_tets_and_tris.merge(all_tris);
20✔
2825

2826
  // create a kd-tree instance
2827
  write_message(
20✔
2828
    7, "Building adaptive k-d tree for tet mesh with ID {}...", id_);
20✔
2829
  kdtree_ = make_unique<moab::AdaptiveKDTree>(mbi_.get());
20✔
2830

2831
  // Determine what options to use
2832
  std::ostringstream options_stream;
20✔
2833
  if (options_.empty()) {
20✔
2834
    options_stream << "MAX_DEPTH=20;PLANE_SET=2;";
4✔
2835
  } else {
2836
    options_stream << options_;
16✔
2837
  }
2838
  moab::FileOptions file_opts(options_stream.str().c_str());
20✔
2839

2840
  // Build the k-d tree
2841
  rval = kdtree_->build_tree(all_tets_and_tris, &kdtree_root_, &file_opts);
20✔
2842
  if (rval != moab::MB_SUCCESS) {
20!
2843
    fatal_error("Failed to construct KDTree for the "
2844
                "unstructured mesh file: " +
2845
                filename_);
×
2846
  }
2847
}
20✔
2848

2849
void MOABMesh::intersect_track(const moab::CartVect& start,
1,543,564✔
2850
  const moab::CartVect& dir, double track_len, vector<double>& hits) const
2851
{
2852
  hits.clear();
1,543,564✔
2853

2854
  moab::ErrorCode rval;
2855
  vector<moab::EntityHandle> tris;
1,543,564✔
2856
  // get all intersections with triangles in the tet mesh
2857
  // (distances are relative to the start point, not the previous
2858
  // intersection)
2859
  rval = kdtree_->ray_intersect_triangles(kdtree_root_, FP_COINCIDENT,
1,543,564✔
2860
    dir.array(), start.array(), tris, hits, 0, track_len);
2861
  if (rval != moab::MB_SUCCESS) {
1,543,564!
2862
    fatal_error(
2863
      "Failed to compute intersections on unstructured mesh: " + filename_);
×
2864
  }
2865

2866
  // remove duplicate intersection distances
2867
  std::unique(hits.begin(), hits.end());
1,543,564✔
2868

2869
  // sorts by first component of std::pair by default
2870
  std::sort(hits.begin(), hits.end());
1,543,564✔
2871
}
1,543,564✔
2872

2873
void MOABMesh::bins_crossed(Position r0, Position r1, const Direction& u,
1,543,564✔
2874
  vector<int>& bins, vector<double>& lengths) const
2875
{
2876
  moab::CartVect start(r0.x, r0.y, r0.z);
1,543,564✔
2877
  moab::CartVect end(r1.x, r1.y, r1.z);
1,543,564✔
2878
  moab::CartVect dir(u.x, u.y, u.z);
1,543,564✔
2879
  dir.normalize();
1,543,564✔
2880

2881
  double track_len = (end - start).length();
1,543,564✔
2882
  if (track_len == 0.0)
1,543,564!
2883
    return;
721,692✔
2884

2885
  start -= TINY_BIT * dir;
1,543,564✔
2886
  end += TINY_BIT * dir;
1,543,564✔
2887

2888
  vector<double> hits;
1,543,564✔
2889
  intersect_track(start, dir, track_len, hits);
1,543,564✔
2890

2891
  bins.clear();
1,543,564✔
2892
  lengths.clear();
1,543,564✔
2893

2894
  // if there are no intersections the track may lie entirely
2895
  // within a single tet. If this is the case, apply entire
2896
  // score to that tet and return.
2897
  if (hits.size() == 0) {
1,543,564✔
2898
    Position midpoint = r0 + u * (track_len * 0.5);
721,692✔
2899
    int bin = this->get_bin(midpoint);
721,692✔
2900
    if (bin != -1) {
721,692✔
2901
      bins.push_back(bin);
242,866✔
2902
      lengths.push_back(1.0);
242,866✔
2903
    }
2904
    return;
721,692✔
2905
  }
2906

2907
  // for each segment in the set of tracks, try to look up a tet
2908
  // at the midpoint of the segment
2909
  Position current = r0;
821,872✔
2910
  double last_dist = 0.0;
821,872✔
2911
  for (const auto& hit : hits) {
5,516,019✔
2912
    // get the segment length
2913
    double segment_length = hit - last_dist;
4,694,147✔
2914
    last_dist = hit;
4,694,147✔
2915
    // find the midpoint of this segment
2916
    Position midpoint = current + u * (segment_length * 0.5);
4,694,147✔
2917
    // try to find a tet for this position
2918
    int bin = this->get_bin(midpoint);
4,694,147✔
2919

2920
    // determine the start point for this segment
2921
    current = r0 + u * hit;
4,694,147✔
2922

2923
    if (bin == -1) {
4,694,147✔
2924
      continue;
20,522✔
2925
    }
2926

2927
    bins.push_back(bin);
4,673,625✔
2928
    lengths.push_back(segment_length / track_len);
4,673,625✔
2929
  }
2930

2931
  // tally remaining portion of track after last hit if
2932
  // the last segment of the track is in the mesh but doesn't
2933
  // reach the other side of the tet
2934
  if (hits.back() < track_len) {
821,872!
2935
    Position segment_start = r0 + u * hits.back();
821,872✔
2936
    double segment_length = track_len - hits.back();
821,872✔
2937
    Position midpoint = segment_start + u * (segment_length * 0.5);
821,872✔
2938
    int bin = this->get_bin(midpoint);
821,872✔
2939
    if (bin != -1) {
821,872✔
2940
      bins.push_back(bin);
766,509✔
2941
      lengths.push_back(segment_length / track_len);
766,509✔
2942
    }
2943
  }
2944
};
1,543,564✔
2945

2946
moab::EntityHandle MOABMesh::get_tet(const Position& r) const
7,317,030✔
2947
{
2948
  moab::CartVect pos(r.x, r.y, r.z);
7,317,030✔
2949
  // find the leaf of the kd-tree for this position
2950
  moab::AdaptiveKDTreeIter kdtree_iter;
7,317,030✔
2951
  moab::ErrorCode rval = kdtree_->point_search(pos.array(), kdtree_iter);
7,317,030✔
2952
  if (rval != moab::MB_SUCCESS) {
7,317,030✔
2953
    return 0;
1,011,877✔
2954
  }
2955

2956
  // retrieve the tet elements of this leaf
2957
  moab::EntityHandle leaf = kdtree_iter.handle();
6,305,153✔
2958
  moab::Range tets;
6,305,153✔
2959
  rval = mbi_->get_entities_by_dimension(leaf, 3, tets, false);
6,305,153✔
2960
  if (rval != moab::MB_SUCCESS) {
6,305,153!
2961
    warning("MOAB error finding tets.");
×
2962
  }
2963

2964
  // loop over the tets in this leaf, returning the containing tet if found
2965
  for (const auto& tet : tets) {
260,209,001✔
2966
    if (point_in_tet(pos, tet)) {
260,206,154✔
2967
      return tet;
6,302,306✔
2968
    }
2969
  }
2970

2971
  // if no tet is found, return an invalid handle
2972
  return 0;
2,847✔
2973
}
7,317,030✔
2974

2975
double MOABMesh::volume(int bin) const
167,856✔
2976
{
2977
  return tet_volume(get_ent_handle_from_bin(bin));
167,856✔
2978
}
2979

2980
std::string MOABMesh::library() const
32✔
2981
{
2982
  return mesh_lib_type;
32✔
2983
}
2984

2985
// Sample position within a tet for MOAB type tets
2986
Position MOABMesh::sample_element(int32_t bin, uint64_t* seed) const
200,410✔
2987
{
2988

2989
  moab::EntityHandle tet_ent = get_ent_handle_from_bin(bin);
200,410✔
2990

2991
  // Get vertex coordinates for MOAB tet
2992
  const moab::EntityHandle* conn1;
2993
  int conn1_size;
2994
  moab::ErrorCode rval = mbi_->get_connectivity(tet_ent, conn1, conn1_size);
200,410✔
2995
  if (rval != moab::MB_SUCCESS || conn1_size != 4) {
200,410!
2996
    fatal_error(fmt::format(
×
2997
      "Failed to get tet connectivity or connectivity size ({}) is invalid.",
2998
      conn1_size));
2999
  }
3000
  moab::CartVect p[4];
1,002,050✔
3001
  rval = mbi_->get_coords(conn1, conn1_size, p[0].array());
200,410✔
3002
  if (rval != moab::MB_SUCCESS) {
200,410!
3003
    fatal_error("Failed to get tet coords");
3004
  }
3005

3006
  std::array<Position, 4> tet_verts;
200,410✔
3007
  for (int i = 0; i < 4; i++) {
1,002,050✔
3008
    tet_verts[i] = {p[i][0], p[i][1], p[i][2]};
801,640✔
3009
  }
3010
  // Samples position within tet using Barycentric stuff
3011
  return this->sample_tet(tet_verts, seed);
400,820✔
3012
}
3013

3014
double MOABMesh::tet_volume(moab::EntityHandle tet) const
167,856✔
3015
{
3016
  vector<moab::EntityHandle> conn;
167,856✔
3017
  moab::ErrorCode rval = mbi_->get_connectivity(&tet, 1, conn);
167,856✔
3018
  if (rval != moab::MB_SUCCESS) {
167,856!
3019
    fatal_error("Failed to get tet connectivity");
3020
  }
3021

3022
  moab::CartVect p[4];
839,280✔
3023
  rval = mbi_->get_coords(conn.data(), conn.size(), p[0].array());
167,856✔
3024
  if (rval != moab::MB_SUCCESS) {
167,856!
3025
    fatal_error("Failed to get tet coords");
3026
  }
3027

3028
  return 1.0 / 6.0 * (((p[1] - p[0]) * (p[2] - p[0])) % (p[3] - p[0]));
335,712✔
3029
}
167,856✔
3030

3031
int MOABMesh::get_bin(Position r) const
7,317,030✔
3032
{
3033
  moab::EntityHandle tet = get_tet(r);
7,317,030✔
3034
  if (tet == 0) {
7,317,030✔
3035
    return -1;
1,014,724✔
3036
  } else {
3037
    return get_bin_from_ent_handle(tet);
6,302,306✔
3038
  }
3039
}
3040

3041
void MOABMesh::compute_barycentric_data(const moab::Range& tets)
20✔
3042
{
3043
  moab::ErrorCode rval;
3044

3045
  baryc_data_.clear();
20✔
3046
  baryc_data_.resize(tets.size());
20✔
3047

3048
  // compute the barycentric data for each tet element
3049
  // and store it as a 3x3 matrix
3050
  for (auto& tet : tets) {
239,732✔
3051
    vector<moab::EntityHandle> verts;
239,712✔
3052
    rval = mbi_->get_connectivity(&tet, 1, verts);
239,712✔
3053
    if (rval != moab::MB_SUCCESS) {
239,712!
3054
      fatal_error("Failed to get connectivity of tet on umesh: " + filename_);
×
3055
    }
3056

3057
    moab::CartVect p[4];
1,198,560✔
3058
    rval = mbi_->get_coords(verts.data(), verts.size(), p[0].array());
239,712✔
3059
    if (rval != moab::MB_SUCCESS) {
239,712!
3060
      fatal_error("Failed to get coordinates of a tet in umesh: " + filename_);
×
3061
    }
3062

3063
    moab::Matrix3 a(p[1] - p[0], p[2] - p[0], p[3] - p[0], true);
239,712✔
3064

3065
    // invert now to avoid this cost later
3066
    a = a.transpose().inverse();
239,712✔
3067
    baryc_data_.at(get_bin_from_ent_handle(tet)) = a;
239,712✔
3068
  }
239,712✔
3069
}
20✔
3070

3071
bool MOABMesh::point_in_tet(
260,206,154✔
3072
  const moab::CartVect& r, moab::EntityHandle tet) const
3073
{
3074

3075
  moab::ErrorCode rval;
3076

3077
  // get tet vertices
3078
  vector<moab::EntityHandle> verts;
260,206,154✔
3079
  rval = mbi_->get_connectivity(&tet, 1, verts);
260,206,154✔
3080
  if (rval != moab::MB_SUCCESS) {
260,206,154!
3081
    warning("Failed to get vertices of tet in umesh: " + filename_);
×
3082
    return false;
3083
  }
3084

3085
  // first vertex is used as a reference point for the barycentric data -
3086
  // retrieve its coordinates
3087
  moab::CartVect p_zero;
260,206,154✔
3088
  rval = mbi_->get_coords(verts.data(), 1, p_zero.array());
260,206,154✔
3089
  if (rval != moab::MB_SUCCESS) {
260,206,154!
3090
    warning("Failed to get coordinates of a vertex in "
×
3091
            "unstructured mesh: " +
3092
            filename_);
×
3093
    return false;
3094
  }
3095

3096
  // look up barycentric data
3097
  int idx = get_bin_from_ent_handle(tet);
260,206,154✔
3098
  const moab::Matrix3& a_inv = baryc_data_[idx];
260,206,154✔
3099

3100
  moab::CartVect bary_coords = a_inv * (r - p_zero);
260,206,154✔
3101

3102
  return (bary_coords[0] >= 0.0 && bary_coords[1] >= 0.0 &&
421,413,584✔
3103
          bary_coords[2] >= 0.0 &&
443,101,423✔
3104
          bary_coords[0] + bary_coords[1] + bary_coords[2] <= 1.0);
281,893,993✔
3105
}
260,206,154✔
3106

3107
int MOABMesh::get_bin_from_index(int idx) const
3108
{
3109
  if (idx >= n_bins()) {
×
3110
    fatal_error(fmt::format("Invalid bin index: {}", idx));
×
3111
  }
3112
  return ehs_[idx] - ehs_[0];
3113
}
3114

3115
int MOABMesh::get_index(const Position& r, bool* in_mesh) const
3116
{
3117
  int bin = get_bin(r);
3118
  *in_mesh = bin != -1;
3119
  return bin;
3120
}
3121

3122
int MOABMesh::get_index_from_bin(int bin) const
3123
{
3124
  return bin;
3125
}
3126

3127
std::pair<vector<double>, vector<double>> MOABMesh::plot(
3128
  Position plot_ll, Position plot_ur) const
3129
{
3130
  // TODO: Implement mesh lines
3131
  return {};
3132
}
3133

3134
int MOABMesh::get_vert_idx_from_handle(moab::EntityHandle vert) const
815,424✔
3135
{
3136
  int idx = vert - verts_[0];
815,424✔
3137
  if (idx >= n_vertices()) {
815,424!
3138
    fatal_error(
3139
      fmt::format("Invalid vertex idx {} (# vertices {})", idx, n_vertices()));
×
3140
  }
3141
  return idx;
815,424✔
3142
}
3143

3144
int MOABMesh::get_bin_from_ent_handle(moab::EntityHandle eh) const
266,748,172✔
3145
{
3146
  int bin = eh - ehs_[0];
266,748,172✔
3147
  if (bin >= n_bins()) {
266,748,172!
3148
    fatal_error(fmt::format("Invalid bin: {}", bin));
×
3149
  }
3150
  return bin;
266,748,172✔
3151
}
3152

3153
moab::EntityHandle MOABMesh::get_ent_handle_from_bin(int bin) const
572,122✔
3154
{
3155
  if (bin >= n_bins()) {
572,122!
3156
    fatal_error(fmt::format("Invalid bin index: ", bin));
×
3157
  }
3158
  return ehs_[0] + bin;
572,122✔
3159
}
3160

3161
int MOABMesh::n_bins() const
267,524,219✔
3162
{
3163
  return ehs_.size();
267,524,219✔
3164
}
3165

3166
int MOABMesh::n_surface_bins() const
3167
{
3168
  // collect all triangles in the set of tets for this mesh
3169
  moab::Range tris;
×
3170
  moab::ErrorCode rval;
3171
  rval = mbi_->get_entities_by_type(0, moab::MBTRI, tris);
×
3172
  if (rval != moab::MB_SUCCESS) {
×
3173
    warning("Failed to get all triangles in the mesh instance");
×
3174
    return -1;
3175
  }
3176
  return 2 * tris.size();
×
3177
}
3178

3179
Position MOABMesh::centroid(int bin) const
3180
{
3181
  moab::ErrorCode rval;
3182

3183
  auto tet = this->get_ent_handle_from_bin(bin);
×
3184

3185
  // look up the tet connectivity
3186
  vector<moab::EntityHandle> conn;
3187
  rval = mbi_->get_connectivity(&tet, 1, conn);
×
3188
  if (rval != moab::MB_SUCCESS) {
×
3189
    warning("Failed to get connectivity of a mesh element.");
×
3190
    return {};
3191
  }
3192

3193
  // get the coordinates
3194
  vector<moab::CartVect> coords(conn.size());
×
3195
  rval = mbi_->get_coords(conn.data(), conn.size(), coords[0].array());
×
3196
  if (rval != moab::MB_SUCCESS) {
×
3197
    warning("Failed to get the coordinates of a mesh element.");
×
3198
    return {};
3199
  }
3200

3201
  // compute the centroid of the element vertices
3202
  moab::CartVect centroid(0.0, 0.0, 0.0);
3203
  for (const auto& coord : coords) {
×
3204
    centroid += coord;
3205
  }
3206
  centroid /= double(coords.size());
3207

3208
  return {centroid[0], centroid[1], centroid[2]};
3209
}
3210

3211
int MOABMesh::n_vertices() const
845,761✔
3212
{
3213
  return verts_.size();
845,761✔
3214
}
3215

3216
Position MOABMesh::vertex(int id) const
86,199✔
3217
{
3218

3219
  moab::ErrorCode rval;
3220

3221
  moab::EntityHandle vert = verts_[id];
86,199✔
3222

3223
  moab::CartVect coords;
86,199✔
3224
  rval = mbi_->get_coords(&vert, 1, coords.array());
86,199✔
3225
  if (rval != moab::MB_SUCCESS) {
86,199!
3226
    fatal_error("Failed to get the coordinates of a vertex.");
3227
  }
3228

3229
  return {coords[0], coords[1], coords[2]};
172,398✔
3230
}
3231

3232
std::vector<int> MOABMesh::connectivity(int bin) const
203,856✔
3233
{
3234
  moab::ErrorCode rval;
3235

3236
  auto tet = get_ent_handle_from_bin(bin);
203,856✔
3237

3238
  // look up the tet connectivity
3239
  vector<moab::EntityHandle> conn;
203,856✔
3240
  rval = mbi_->get_connectivity(&tet, 1, conn);
203,856✔
3241
  if (rval != moab::MB_SUCCESS) {
203,856!
3242
    fatal_error("Failed to get connectivity of a mesh element.");
3243
    return {};
3244
  }
3245

3246
  std::vector<int> verts(4);
203,856✔
3247
  for (int i = 0; i < verts.size(); i++) {
1,019,280✔
3248
    verts[i] = get_vert_idx_from_handle(conn[i]);
815,424✔
3249
  }
3250

3251
  return verts;
203,856✔
3252
}
203,856✔
3253

3254
std::pair<moab::Tag, moab::Tag> MOABMesh::get_score_tags(
3255
  std::string score) const
3256
{
3257
  moab::ErrorCode rval;
3258
  // add a tag to the mesh
3259
  // all scores are treated as a single value
3260
  // with an uncertainty
3261
  moab::Tag value_tag;
3262

3263
  // create the value tag if not present and get handle
3264
  double default_val = 0.0;
3265
  auto val_string = score + "_mean";
×
3266
  rval = mbi_->tag_get_handle(val_string.c_str(), 1, moab::MB_TYPE_DOUBLE,
×
3267
    value_tag, moab::MB_TAG_DENSE | moab::MB_TAG_CREAT, &default_val);
3268
  if (rval != moab::MB_SUCCESS) {
×
3269
    auto msg =
3270
      fmt::format("Could not create or retrieve the value tag for the score {}"
3271
                  " on unstructured mesh {}",
3272
        score, id_);
×
3273
    fatal_error(msg);
3274
  }
3275

3276
  // create the std dev tag if not present and get handle
3277
  moab::Tag error_tag;
3278
  std::string err_string = score + "_std_dev";
×
3279
  rval = mbi_->tag_get_handle(err_string.c_str(), 1, moab::MB_TYPE_DOUBLE,
×
3280
    error_tag, moab::MB_TAG_DENSE | moab::MB_TAG_CREAT, &default_val);
3281
  if (rval != moab::MB_SUCCESS) {
×
3282
    auto msg =
3283
      fmt::format("Could not create or retrieve the error tag for the score {}"
3284
                  " on unstructured mesh {}",
3285
        score, id_);
×
3286
    fatal_error(msg);
3287
  }
3288

3289
  // return the populated tag handles
3290
  return {value_tag, error_tag};
3291
}
3292

3293
void MOABMesh::add_score(const std::string& score)
3294
{
3295
  auto score_tags = get_score_tags(score);
×
3296
  tag_names_.push_back(score);
×
3297
}
3298

3299
void MOABMesh::remove_scores()
3300
{
3301
  for (const auto& name : tag_names_) {
×
3302
    auto value_name = name + "_mean";
×
3303
    moab::Tag tag;
3304
    moab::ErrorCode rval = mbi_->tag_get_handle(value_name.c_str(), tag);
×
3305
    if (rval != moab::MB_SUCCESS)
×
3306
      return;
3307

3308
    rval = mbi_->tag_delete(tag);
×
3309
    if (rval != moab::MB_SUCCESS) {
×
3310
      auto msg = fmt::format("Failed to delete mesh tag for the score {}"
3311
                             " on unstructured mesh {}",
3312
        name, id_);
×
3313
      fatal_error(msg);
3314
    }
3315

3316
    auto std_dev_name = name + "_std_dev";
×
3317
    rval = mbi_->tag_get_handle(std_dev_name.c_str(), tag);
×
3318
    if (rval != moab::MB_SUCCESS) {
×
3319
      auto msg =
3320
        fmt::format("Std. Dev. mesh tag does not exist for the score {}"
3321
                    " on unstructured mesh {}",
3322
          name, id_);
×
3323
    }
3324

3325
    rval = mbi_->tag_delete(tag);
×
3326
    if (rval != moab::MB_SUCCESS) {
×
3327
      auto msg = fmt::format("Failed to delete mesh tag for the score {}"
3328
                             " on unstructured mesh {}",
3329
        name, id_);
×
3330
      fatal_error(msg);
3331
    }
3332
  }
×
3333
  tag_names_.clear();
3334
}
3335

3336
void MOABMesh::set_score_data(const std::string& score,
3337
  const vector<double>& values, const vector<double>& std_dev)
3338
{
3339
  auto score_tags = this->get_score_tags(score);
×
3340

3341
  moab::ErrorCode rval;
3342
  // set the score value
3343
  rval = mbi_->tag_set_data(score_tags.first, ehs_, values.data());
×
3344
  if (rval != moab::MB_SUCCESS) {
×
3345
    auto msg = fmt::format("Failed to set the tally value for score '{}' "
3346
                           "on unstructured mesh {}",
3347
      score, id_);
×
3348
    warning(msg);
×
3349
  }
3350

3351
  // set the error value
3352
  rval = mbi_->tag_set_data(score_tags.second, ehs_, std_dev.data());
×
3353
  if (rval != moab::MB_SUCCESS) {
×
3354
    auto msg = fmt::format("Failed to set the tally error for score '{}' "
3355
                           "on unstructured mesh {}",
3356
      score, id_);
×
3357
    warning(msg);
×
3358
  }
3359
}
3360

3361
void MOABMesh::write(const std::string& base_filename) const
3362
{
3363
  // add extension to the base name
3364
  auto filename = base_filename + ".vtk";
×
3365
  write_message(5, "Writing unstructured mesh {}...", filename);
×
3366
  filename = settings::path_output + filename;
×
3367

3368
  // write the tetrahedral elements of the mesh only
3369
  // to avoid clutter from zero-value data on other
3370
  // elements during visualization
3371
  moab::ErrorCode rval;
3372
  rval = mbi_->write_mesh(filename.c_str(), &tetset_, 1);
×
3373
  if (rval != moab::MB_SUCCESS) {
×
3374
    auto msg = fmt::format("Failed to write unstructured mesh {}", id_);
×
3375
    warning(msg);
×
3376
  }
3377
}
3378

3379
#endif
3380

3381
#ifdef OPENMC_LIBMESH_ENABLED
3382

3383
const std::string LibMesh::mesh_lib_type = "libmesh";
3384

3385
LibMesh::LibMesh(pugi::xml_node node) : UnstructuredMesh(node)
23✔
3386
{
3387
  // filename_ and length_multiplier_ will already be set by the
3388
  // UnstructuredMesh constructor
3389
  set_mesh_pointer_from_filename(filename_);
23✔
3390
  set_length_multiplier(length_multiplier_);
23✔
3391
  initialize();
23✔
3392
}
23✔
3393

3394
LibMesh::LibMesh(hid_t group) : UnstructuredMesh(group)
×
3395
{
3396
  // filename_ and length_multiplier_ will already be set by the
3397
  // UnstructuredMesh constructor
3398
  set_mesh_pointer_from_filename(filename_);
×
3399
  set_length_multiplier(length_multiplier_);
×
3400
  initialize();
×
3401
}
3402

3403
// create the mesh from a pointer to a libMesh Mesh
3404
LibMesh::LibMesh(libMesh::MeshBase& input_mesh, double length_multiplier)
×
3405
{
3406
  if (!dynamic_cast<libMesh::ReplicatedMesh*>(&input_mesh)) {
×
3407
    fatal_error("At present LibMesh tallies require a replicated mesh. Please "
3408
                "ensure 'input_mesh' is a libMesh::ReplicatedMesh.");
3409
  }
3410

3411
  m_ = &input_mesh;
3412
  set_length_multiplier(length_multiplier);
×
3413
  initialize();
×
3414
}
3415

3416
// create the mesh from an input file
3417
LibMesh::LibMesh(const std::string& filename, double length_multiplier)
×
3418
{
3419
  n_dimension_ = 3;
3420
  set_mesh_pointer_from_filename(filename);
×
3421
  set_length_multiplier(length_multiplier);
×
3422
  initialize();
×
3423
}
3424

3425
void LibMesh::set_mesh_pointer_from_filename(const std::string& filename)
23✔
3426
{
3427
  filename_ = filename;
23✔
3428
  unique_m_ =
3429
    make_unique<libMesh::ReplicatedMesh>(*settings::libmesh_comm, n_dimension_);
23✔
3430
  m_ = unique_m_.get();
23✔
3431
  m_->read(filename_);
23✔
3432
}
23✔
3433

3434
// build a libMesh equation system for storing values
3435
void LibMesh::build_eqn_sys()
15✔
3436
{
3437
  eq_system_name_ = fmt::format("mesh_{}_system", id_);
30✔
3438
  equation_systems_ = make_unique<libMesh::EquationSystems>(*m_);
15✔
3439
  libMesh::ExplicitSystem& eq_sys =
3440
    equation_systems_->add_system<libMesh::ExplicitSystem>(eq_system_name_);
15✔
3441
}
15✔
3442

3443
// intialize from mesh file
3444
void LibMesh::initialize()
23✔
3445
{
3446
  if (!settings::libmesh_comm) {
23!
3447
    fatal_error("Attempting to use an unstructured mesh without a libMesh "
3448
                "communicator.");
3449
  }
3450

3451
  // assuming that unstructured meshes used in OpenMC are 3D
3452
  n_dimension_ = 3;
23✔
3453

3454
  if (length_multiplier_ > 0.0) {
23!
3455
    libMesh::MeshTools::Modification::scale(*m_, length_multiplier_);
×
3456
  }
3457
  // if OpenMC is managing the libMesh::MeshBase instance, prepare the mesh.
3458
  // Otherwise assume that it is prepared by its owning application
3459
  if (unique_m_) {
23!
3460
    m_->prepare_for_use();
23✔
3461
  }
3462

3463
  // ensure that the loaded mesh is 3 dimensional
3464
  if (m_->mesh_dimension() != n_dimension_) {
23!
3465
    fatal_error(fmt::format("Mesh file {} specified for use in an unstructured "
3466
                            "mesh is not a 3D mesh.",
3467
      filename_));
3468
  }
3469

3470
  for (int i = 0; i < num_threads(); i++) {
69✔
3471
    pl_.emplace_back(m_->sub_point_locator());
46✔
3472
    pl_.back()->set_contains_point_tol(FP_COINCIDENT);
46✔
3473
    pl_.back()->enable_out_of_mesh_mode();
46✔
3474
  }
3475

3476
  // store first element in the mesh to use as an offset for bin indices
3477
  auto first_elem = *m_->elements_begin();
23✔
3478
  first_element_id_ = first_elem->id();
23✔
3479

3480
  // bounding box for the mesh for quick rejection checks
3481
  bbox_ = libMesh::MeshTools::create_bounding_box(*m_);
23✔
3482
  libMesh::Point ll = bbox_.min();
23✔
3483
  libMesh::Point ur = bbox_.max();
23✔
3484
  lower_left_ = {ll(0), ll(1), ll(2)};
23✔
3485
  upper_right_ = {ur(0), ur(1), ur(2)};
23✔
3486
}
23✔
3487

3488
// Sample position within a tet for LibMesh type tets
3489
Position LibMesh::sample_element(int32_t bin, uint64_t* seed) const
400,820✔
3490
{
3491
  const auto& elem = get_element_from_bin(bin);
400,820✔
3492
  // Get tet vertex coordinates from LibMesh
3493
  std::array<Position, 4> tet_verts;
400,820✔
3494
  for (int i = 0; i < elem.n_nodes(); i++) {
2,004,100✔
3495
    auto node_ref = elem.node_ref(i);
1,603,280✔
3496
    tet_verts[i] = {node_ref(0), node_ref(1), node_ref(2)};
1,603,280✔
3497
  }
1,603,280✔
3498
  // Samples position within tet using Barycentric coordinates
3499
  return this->sample_tet(tet_verts, seed);
801,640✔
3500
}
3501

3502
Position LibMesh::centroid(int bin) const
3503
{
3504
  const auto& elem = this->get_element_from_bin(bin);
×
3505
  auto centroid = elem.vertex_average();
×
3506
  return {centroid(0), centroid(1), centroid(2)};
3507
}
3508

3509
int LibMesh::n_vertices() const
39,978✔
3510
{
3511
  return m_->n_nodes();
39,978✔
3512
}
3513

3514
Position LibMesh::vertex(int vertex_id) const
39,942✔
3515
{
3516
  const auto node_ref = m_->node_ref(vertex_id);
39,942✔
3517
  return {node_ref(0), node_ref(1), node_ref(2)};
79,884✔
3518
}
39,942✔
3519

3520
std::vector<int> LibMesh::connectivity(int elem_id) const
265,856✔
3521
{
3522
  std::vector<int> conn;
265,856✔
3523
  const auto* elem_ptr = m_->elem_ptr(elem_id);
265,856✔
3524
  for (int i = 0; i < elem_ptr->n_nodes(); i++) {
1,337,280✔
3525
    conn.push_back(elem_ptr->node_id(i));
1,071,424✔
3526
  }
3527
  return conn;
265,856✔
3528
}
3529

3530
std::string LibMesh::library() const
33✔
3531
{
3532
  return mesh_lib_type;
33✔
3533
}
3534

3535
int LibMesh::n_bins() const
1,784,287✔
3536
{
3537
  return m_->n_elem();
1,784,287✔
3538
}
3539

3540
int LibMesh::n_surface_bins() const
3541
{
3542
  int n_bins = 0;
3543
  for (int i = 0; i < this->n_bins(); i++) {
×
3544
    const libMesh::Elem& e = get_element_from_bin(i);
3545
    n_bins += e.n_faces();
3546
    // if this is a boundary element, it will only be visited once,
3547
    // the number of surface bins is incremented to
3548
    for (auto neighbor_ptr : e.neighbor_ptr_range()) {
×
3549
      // null neighbor pointer indicates a boundary face
3550
      if (!neighbor_ptr) {
×
3551
        n_bins++;
3552
      }
3553
    }
3554
  }
3555
  return n_bins;
3556
}
3557

3558
void LibMesh::add_score(const std::string& var_name)
15✔
3559
{
3560
  if (!equation_systems_) {
15!
3561
    build_eqn_sys();
15✔
3562
  }
3563

3564
  // check if this is a new variable
3565
  std::string value_name = var_name + "_mean";
15✔
3566
  if (!variable_map_.count(value_name)) {
15!
3567
    auto& eqn_sys = equation_systems_->get_system(eq_system_name_);
15✔
3568
    auto var_num =
3569
      eqn_sys.add_variable(value_name, libMesh::CONSTANT, libMesh::MONOMIAL);
15✔
3570
    variable_map_[value_name] = var_num;
15✔
3571
  }
3572

3573
  std::string std_dev_name = var_name + "_std_dev";
15✔
3574
  // check if this is a new variable
3575
  if (!variable_map_.count(std_dev_name)) {
15!
3576
    auto& eqn_sys = equation_systems_->get_system(eq_system_name_);
15✔
3577
    auto var_num =
3578
      eqn_sys.add_variable(std_dev_name, libMesh::CONSTANT, libMesh::MONOMIAL);
15✔
3579
    variable_map_[std_dev_name] = var_num;
15✔
3580
  }
3581
}
15✔
3582

3583
void LibMesh::remove_scores()
15✔
3584
{
3585
  if (equation_systems_) {
15!
3586
    auto& eqn_sys = equation_systems_->get_system(eq_system_name_);
15✔
3587
    eqn_sys.clear();
15✔
3588
    variable_map_.clear();
15✔
3589
  }
3590
}
15✔
3591

3592
void LibMesh::set_score_data(const std::string& var_name,
15✔
3593
  const vector<double>& values, const vector<double>& std_dev)
3594
{
3595
  if (!equation_systems_) {
15!
3596
    build_eqn_sys();
×
3597
  }
3598

3599
  auto& eqn_sys = equation_systems_->get_system(eq_system_name_);
15✔
3600

3601
  if (!eqn_sys.is_initialized()) {
15!
3602
    equation_systems_->init();
15✔
3603
  }
3604

3605
  const libMesh::DofMap& dof_map = eqn_sys.get_dof_map();
15✔
3606

3607
  // look up the value variable
3608
  std::string value_name = var_name + "_mean";
15✔
3609
  unsigned int value_num = variable_map_.at(value_name);
15✔
3610
  // look up the std dev variable
3611
  std::string std_dev_name = var_name + "_std_dev";
15✔
3612
  unsigned int std_dev_num = variable_map_.at(std_dev_name);
15✔
3613

3614
  for (auto it = m_->local_elements_begin(); it != m_->local_elements_end();
97,871✔
3615
       it++) {
3616
    if (!(*it)->active()) {
97,856!
3617
      continue;
3618
    }
3619

3620
    auto bin = get_bin_from_element(*it);
97,856✔
3621

3622
    // set value
3623
    vector<libMesh::dof_id_type> value_dof_indices;
97,856✔
3624
    dof_map.dof_indices(*it, value_dof_indices, value_num);
97,856✔
3625
    assert(value_dof_indices.size() == 1);
3626
    eqn_sys.solution->set(value_dof_indices[0], values.at(bin));
97,856✔
3627

3628
    // set std dev
3629
    vector<libMesh::dof_id_type> std_dev_dof_indices;
97,856✔
3630
    dof_map.dof_indices(*it, std_dev_dof_indices, std_dev_num);
97,856✔
3631
    assert(std_dev_dof_indices.size() == 1);
3632
    eqn_sys.solution->set(std_dev_dof_indices[0], std_dev.at(bin));
97,856✔
3633
  }
97,871✔
3634
}
15✔
3635

3636
void LibMesh::write(const std::string& filename) const
15✔
3637
{
3638
  write_message(fmt::format(
15✔
3639
    "Writing file: {}.e for unstructured mesh {}", filename, this->id_));
15✔
3640
  libMesh::ExodusII_IO exo(*m_);
15✔
3641
  std::set<std::string> systems_out = {eq_system_name_};
45✔
3642
  exo.write_discontinuous_exodusII(
15✔
3643
    filename + ".e", *equation_systems_, &systems_out);
30✔
3644
}
15✔
3645

3646
void LibMesh::bins_crossed(Position r0, Position r1, const Direction& u,
3647
  vector<int>& bins, vector<double>& lengths) const
3648
{
3649
  // TODO: Implement triangle crossings here
3650
  fatal_error("Tracklength tallies on libMesh instances are not implemented.");
3651
}
3652

3653
int LibMesh::get_bin(Position r) const
2,340,484✔
3654
{
3655
  // look-up a tet using the point locator
3656
  libMesh::Point p(r.x, r.y, r.z);
2,340,484✔
3657

3658
  // quick rejection check
3659
  if (!bbox_.contains_point(p)) {
2,340,484✔
3660
    return -1;
918,796✔
3661
  }
3662

3663
  const auto& point_locator = pl_.at(thread_num());
1,421,688✔
3664

3665
  const auto elem_ptr = (*point_locator)(p);
1,421,688✔
3666
  return elem_ptr ? get_bin_from_element(elem_ptr) : -1;
1,421,688✔
3667
}
2,340,484✔
3668

3669
int LibMesh::get_bin_from_element(const libMesh::Elem* elem) const
1,518,314✔
3670
{
3671
  int bin = elem->id() - first_element_id_;
1,518,314✔
3672
  if (bin >= n_bins() || bin < 0) {
1,518,314!
3673
    fatal_error(fmt::format("Invalid bin: {}", bin));
3674
  }
3675
  return bin;
1,518,314✔
3676
}
3677

3678
std::pair<vector<double>, vector<double>> LibMesh::plot(
3679
  Position plot_ll, Position plot_ur) const
3680
{
3681
  return {};
3682
}
3683

3684
const libMesh::Elem& LibMesh::get_element_from_bin(int bin) const
765,460✔
3685
{
3686
  return m_->elem_ref(bin);
765,460✔
3687
}
3688

3689
double LibMesh::volume(int bin) const
364,640✔
3690
{
3691
  return this->get_element_from_bin(bin).volume();
364,640✔
3692
}
3693

3694
AdaptiveLibMesh::AdaptiveLibMesh(
3695
  libMesh::MeshBase& input_mesh, double length_multiplier)
3696
  : LibMesh(input_mesh, length_multiplier), num_active_(m_->n_active_elem())
×
3697
{
3698
  // if the mesh is adaptive elements aren't guaranteed by libMesh to be
3699
  // contiguous in ID space, so we need to map from bin indices (defined over
3700
  // active elements) to global dof ids
3701
  bin_to_elem_map_.reserve(num_active_);
×
3702
  elem_to_bin_map_.resize(m_->n_elem(), -1);
×
3703
  for (auto it = m_->active_elements_begin(); it != m_->active_elements_end();
×
3704
       it++) {
3705
    auto elem = *it;
×
3706

3707
    bin_to_elem_map_.push_back(elem->id());
×
3708
    elem_to_bin_map_[elem->id()] = bin_to_elem_map_.size() - 1;
×
3709
  }
3710
}
3711

3712
int AdaptiveLibMesh::n_bins() const
3713
{
3714
  return num_active_;
3715
}
3716

3717
void AdaptiveLibMesh::add_score(const std::string& var_name)
3718
{
3719
  warning(fmt::format(
×
3720
    "Exodus output cannot be provided as unstructured mesh {} is adaptive.",
3721
    this->id_));
3722
}
3723

3724
void AdaptiveLibMesh::set_score_data(const std::string& var_name,
3725
  const vector<double>& values, const vector<double>& std_dev)
3726
{
3727
  warning(fmt::format(
×
3728
    "Exodus output cannot be provided as unstructured mesh {} is adaptive.",
3729
    this->id_));
3730
}
3731

3732
void AdaptiveLibMesh::write(const std::string& filename) const
3733
{
3734
  warning(fmt::format(
×
3735
    "Exodus output cannot be provided as unstructured mesh {} is adaptive.",
3736
    this->id_));
3737
}
3738

3739
int AdaptiveLibMesh::get_bin_from_element(const libMesh::Elem* elem) const
3740
{
3741
  int bin = elem_to_bin_map_[elem->id()];
×
3742
  if (bin >= n_bins() || bin < 0) {
×
3743
    fatal_error(fmt::format("Invalid bin: {}", bin));
3744
  }
3745
  return bin;
3746
}
3747

3748
const libMesh::Elem& AdaptiveLibMesh::get_element_from_bin(int bin) const
3749
{
3750
  return m_->elem_ref(bin_to_elem_map_.at(bin));
3751
}
3752

3753
#endif // OPENMC_LIBMESH_ENABLED
3754

3755
//==============================================================================
3756
// Non-member functions
3757
//==============================================================================
3758

3759
void read_meshes(pugi::xml_node root)
11,334✔
3760
{
3761
  std::unordered_set<int> mesh_ids;
11,334✔
3762

3763
  for (auto node : root.children("mesh")) {
14,045✔
3764
    // Check to make sure multiple meshes in the same file don't share IDs
3765
    int id = std::stoi(get_node_value(node, "id"));
2,711✔
3766
    if (contains(mesh_ids, id)) {
2,711!
UNCOV
3767
      fatal_error(fmt::format("Two or more meshes use the same unique ID "
×
3768
                              "'{}' in the same input file",
3769
        id));
3770
    }
3771
    mesh_ids.insert(id);
2,711✔
3772

3773
    // If we've already read a mesh with the same ID in a *different* file,
3774
    // assume it is the same here
3775
    if (model::mesh_map.find(id) != model::mesh_map.end()) {
2,711!
UNCOV
3776
      warning(fmt::format("Mesh with ID={} appears in multiple files.", id));
×
UNCOV
3777
      continue;
×
3778
    }
3779

3780
    std::string mesh_type;
2,711✔
3781
    if (check_for_node(node, "type")) {
2,711✔
3782
      mesh_type = get_node_value(node, "type", true, true);
944✔
3783
    } else {
3784
      mesh_type = "regular";
1,767✔
3785
    }
3786

3787
    // determine the mesh library to use
3788
    std::string mesh_lib;
2,711✔
3789
    if (check_for_node(node, "library")) {
2,711✔
3790
      mesh_lib = get_node_value(node, "library", true, true);
46!
3791
    }
3792

3793
    // Read mesh and add to vector
3794
    if (mesh_type == RegularMesh::mesh_type) {
2,711✔
3795
      model::meshes.push_back(make_unique<RegularMesh>(node));
1,837✔
3796
    } else if (mesh_type == RectilinearMesh::mesh_type) {
874✔
3797
      model::meshes.push_back(make_unique<RectilinearMesh>(node));
114✔
3798
    } else if (mesh_type == CylindricalMesh::mesh_type) {
760✔
3799
      model::meshes.push_back(make_unique<CylindricalMesh>(node));
390✔
3800
    } else if (mesh_type == SphericalMesh::mesh_type) {
370✔
3801
      model::meshes.push_back(make_unique<SphericalMesh>(node));
324✔
3802
#ifdef OPENMC_DAGMC_ENABLED
3803
    } else if (mesh_type == UnstructuredMesh::mesh_type &&
46!
3804
               mesh_lib == MOABMesh::mesh_lib_type) {
23✔
3805
      model::meshes.push_back(make_unique<MOABMesh>(node));
23✔
3806
#endif
3807
#ifdef OPENMC_LIBMESH_ENABLED
3808
    } else if (mesh_type == UnstructuredMesh::mesh_type &&
46!
3809
               mesh_lib == LibMesh::mesh_lib_type) {
23✔
3810
      model::meshes.push_back(make_unique<LibMesh>(node));
23✔
3811
#endif
UNCOV
3812
    } else if (mesh_type == UnstructuredMesh::mesh_type) {
×
UNCOV
3813
      fatal_error("Unstructured mesh support is not enabled or the mesh "
×
3814
                  "library is invalid.");
3815
    } else {
UNCOV
3816
      fatal_error("Invalid mesh type: " + mesh_type);
×
3817
    }
3818

3819
    // Map ID to position in vector
3820
    model::mesh_map[model::meshes.back()->id_] = model::meshes.size() - 1;
2,711✔
3821
  }
2,711✔
3822
}
11,334✔
3823

3824
void read_meshes(hid_t group)
11✔
3825
{
3826
  std::unordered_set<int> mesh_ids;
11✔
3827

3828
  std::array<int, 1> ids;
3829
  read_attribute(group, "ids", ids);
11✔
3830

3831
  for (auto id : ids) {
22✔
3832

3833
    // Check to make sure multiple meshes in the same file don't share IDs
3834
    if (contains(mesh_ids, id)) {
11!
NEW
3835
      fatal_error(fmt::format("Two or more meshes use the same unique ID "
×
3836
                              "'{}' in the same HDF5 input file",
3837
        id));
3838
    }
3839
    mesh_ids.insert(id);
11✔
3840

3841
    // If we've already read a mesh with the same ID in a *different* file,
3842
    // assume it is the same here
3843
    if (model::mesh_map.find(id) != model::mesh_map.end()) {
11!
3844
      warning(fmt::format("Mesh with ID={} appears in multiple files.", id));
11✔
3845
      continue;
11✔
3846
    }
3847

NEW
3848
    std::string name = fmt::format("mesh {}", id);
×
NEW
3849
    hid_t mesh_group = open_group(group, name.c_str());
×
3850

NEW
3851
    std::string mesh_type;
×
NEW
3852
    if (object_exists(mesh_group, "type")) {
×
NEW
3853
      read_dataset(mesh_group, "type", mesh_type);
×
3854
    } else {
NEW
3855
      mesh_type = "regular";
×
3856
    }
3857

3858
    // determine the mesh library to use
NEW
3859
    std::string mesh_lib;
×
NEW
3860
    if (object_exists(mesh_group, "library")) {
×
NEW
3861
      read_dataset(mesh_group, "library", mesh_lib);
×
3862
    }
3863

3864
    // Read mesh and add to vector
NEW
3865
    if (mesh_type == RegularMesh::mesh_type) {
×
NEW
3866
      model::meshes.push_back(make_unique<RegularMesh>(mesh_group));
×
NEW
UNCOV
3867
    } else if (mesh_type == RectilinearMesh::mesh_type) {
×
NEW
UNCOV
3868
      model::meshes.push_back(make_unique<RectilinearMesh>(mesh_group));
×
NEW
3869
    } else if (mesh_type == CylindricalMesh::mesh_type) {
×
NEW
UNCOV
3870
      model::meshes.push_back(make_unique<CylindricalMesh>(mesh_group));
×
NEW
UNCOV
3871
    } else if (mesh_type == SphericalMesh::mesh_type) {
×
NEW
UNCOV
3872
      model::meshes.push_back(make_unique<SphericalMesh>(mesh_group));
×
3873
#ifdef OPENMC_DAGMC_ENABLED
3874
    } else if (mesh_type == UnstructuredMesh::mesh_type &&
×
3875
               mesh_lib == MOABMesh::mesh_lib_type) {
3876
      model::meshes.push_back(make_unique<MOABMesh>(mesh_group));
×
3877
#endif
3878
#ifdef OPENMC_LIBMESH_ENABLED
3879
    } else if (mesh_type == UnstructuredMesh::mesh_type &&
×
3880
               mesh_lib == LibMesh::mesh_lib_type) {
3881
      model::meshes.push_back(make_unique<LibMesh>(mesh_group));
×
3882
#endif
NEW
UNCOV
3883
    } else if (mesh_type == UnstructuredMesh::mesh_type) {
×
NEW
UNCOV
3884
      fatal_error("Unstructured mesh support is not enabled or the mesh "
×
3885
                  "library is invalid.");
3886
    } else {
NEW
UNCOV
3887
      fatal_error("Invalid mesh type: " + mesh_type);
×
3888
    }
3889

3890
    // Map ID to position in vector
NEW
UNCOV
3891
    model::mesh_map[model::meshes.back()->id_] = model::meshes.size() - 1;
×
NEW
UNCOV
3892
  }
×
3893
}
11✔
3894

3895
void meshes_to_hdf5(hid_t group)
6,383✔
3896
{
3897
  // Write number of meshes
3898
  hid_t meshes_group = create_group(group, "meshes");
6,383✔
3899
  int32_t n_meshes = model::meshes.size();
6,383✔
3900
  write_attribute(meshes_group, "n_meshes", n_meshes);
6,383✔
3901

3902
  if (n_meshes > 0) {
6,383✔
3903
    // Write IDs of meshes
3904
    vector<int> ids;
1,959✔
3905
    for (const auto& m : model::meshes) {
4,475✔
3906
      m->to_hdf5(meshes_group);
2,516✔
3907
      ids.push_back(m->id_);
2,516✔
3908
    }
3909
    write_attribute(meshes_group, "ids", ids);
1,959✔
3910
  }
1,959✔
3911

3912
  close_group(meshes_group);
6,383✔
3913
}
6,383✔
3914

3915
void free_memory_mesh()
7,469✔
3916
{
3917
  model::meshes.clear();
7,469✔
3918
  model::mesh_map.clear();
7,469✔
3919
}
7,469✔
3920

3921
extern "C" int n_meshes()
308✔
3922
{
3923
  return model::meshes.size();
308✔
3924
}
3925

3926
} // 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

© 2025 Coveralls, Inc