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

openmc-dev / openmc / 5979019501

25 Aug 2023 06:04PM UTC coverage: 84.434% (-0.003%) from 84.437%
5979019501

push

github

web-flow
Always set html_theme in readthedocs configuration (#2667)

46790 of 55416 relevant lines covered (84.43%)

78088697.44 hits per line

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

92.23
/src/universe.cpp
1
#include "openmc/universe.h"
2

3
#include <set>
4

5
#include "openmc/hdf5_interface.h"
6

7
namespace openmc {
8

9
namespace model {
10

11
std::unordered_map<int32_t, int32_t> universe_map;
12
vector<unique_ptr<Universe>> universes;
13

14
} // namespace model
15

16
//==============================================================================
17
// Universe implementation
18
//==============================================================================
19

20
void Universe::to_hdf5(hid_t universes_group) const
13,775✔
21
{
22
  // Create a group for this universe.
23
  auto group = create_group(universes_group, fmt::format("universe {}", id_));
27,550✔
24

25
  // Write the geometry representation type.
26
  write_string(group, "geom_type", "csg", false);
13,775✔
27

28
  // Write the contained cells.
29
  if (cells_.size() > 0) {
13,775✔
30
    vector<int32_t> cell_ids;
13,775✔
31
    for (auto i_cell : cells_)
37,675✔
32
      cell_ids.push_back(model::cells[i_cell]->id_);
23,900✔
33
    write_dataset(group, "cells", cell_ids);
13,775✔
34
  }
13,775✔
35

36
  close_group(group);
13,775✔
37
}
13,775✔
38

39
bool Universe::find_cell(Particle& p) const
867,590,181✔
40
{
41
  const auto& cells {
42
    !partitioner_ ? cells_ : partitioner_->get_cells(p.r_local(), p.u_local())};
867,590,181✔
43

44
  Position r {p.r_local()};
867,590,181✔
45
  Position u {p.u_local()};
867,590,181✔
46
  auto surf = p.surface();
867,590,181✔
47
  int32_t i_univ = p.lowest_coord().universe;
867,590,181✔
48

49
  for (auto i_cell : cells) {
1,541,554,142✔
50
    if (model::cells[i_cell]->universe_ != i_univ)
1,332,425,152✔
51
      continue;
×
52
    // Check if this cell contains the particle
53
    if (model::cells[i_cell]->contains(r, u, surf)) {
1,332,425,152✔
54
      p.lowest_coord().cell = i_cell;
658,461,191✔
55
      return true;
658,461,191✔
56
    }
57
  }
58
  return false;
209,128,990✔
59
}
60

61
BoundingBox Universe::bounding_box() const
14✔
62
{
63
  BoundingBox bbox = {INFTY, -INFTY, INFTY, -INFTY, INFTY, -INFTY};
14✔
64
  if (cells_.size() == 0) {
14✔
65
    return {};
×
66
  } else {
67
    for (const auto& cell : cells_) {
56✔
68
      auto& c = model::cells[cell];
42✔
69
      bbox |= c->bounding_box();
42✔
70
    }
71
  }
72
  return bbox;
14✔
73
}
74

75
//==============================================================================
76
// UniversePartitioner implementation
77
//==============================================================================
78

79
UniversePartitioner::UniversePartitioner(const Universe& univ)
114✔
80
{
81
  // Define an ordered set of surface indices that point to z-planes.  Use a
82
  // functor to to order the set by the z0_ values of the corresponding planes.
83
  struct compare_surfs {
84
    bool operator()(const int32_t& i_surf, const int32_t& j_surf) const
11,723✔
85
    {
86
      const auto* surf = model::surfaces[i_surf].get();
11,723✔
87
      const auto* zplane = dynamic_cast<const SurfaceZPlane*>(surf);
11,723✔
88
      double zi = zplane->z0_;
11,723✔
89
      surf = model::surfaces[j_surf].get();
11,723✔
90
      zplane = dynamic_cast<const SurfaceZPlane*>(surf);
11,723✔
91
      double zj = zplane->z0_;
11,723✔
92
      return zi < zj;
11,723✔
93
    }
94
  };
95
  std::set<int32_t, compare_surfs> surf_set;
114✔
96

97
  // Find all of the z-planes in this universe.  A set is used here for the
98
  // O(log(n)) insertions that will ensure entries are not repeated.
99
  for (auto i_cell : univ.cells_) {
1,596✔
100
    for (auto token : model::cells[i_cell]->surfaces()) {
6,555✔
101
      auto i_surf = std::abs(token) - 1;
5,073✔
102
      const auto* surf = model::surfaces[i_surf].get();
5,073✔
103
      if (const auto* zplane = dynamic_cast<const SurfaceZPlane*>(surf))
5,073✔
104
        surf_set.insert(i_surf);
2,888✔
105
    }
1,482✔
106
  }
107

108
  // Populate the surfs_ vector from the ordered set.
109
  surfs_.insert(surfs_.begin(), surf_set.begin(), surf_set.end());
114✔
110

111
  // Populate the partition lists.
112
  partitions_.resize(surfs_.size() + 1);
114✔
113
  for (auto i_cell : univ.cells_) {
1,596✔
114
    // It is difficult to determine the bounds of a complex cell, so add complex
115
    // cells to all partitions.
116
    if (!model::cells[i_cell]->is_simple()) {
1,482✔
117
      for (auto& p : partitions_)
×
118
        p.push_back(i_cell);
×
119
      continue;
×
120
    }
121

122
    // Find the tokens for bounding z-planes.
123
    int32_t lower_token = 0, upper_token = 0;
1,482✔
124
    double min_z, max_z;
125
    for (auto token : model::cells[i_cell]->surfaces()) {
6,555✔
126
      const auto* surf = model::surfaces[std::abs(token) - 1].get();
5,073✔
127
      if (const auto* zplane = dynamic_cast<const SurfaceZPlane*>(surf)) {
5,073✔
128
        if (lower_token == 0 || zplane->z0_ < min_z) {
2,888✔
129
          lower_token = token;
1,482✔
130
          min_z = zplane->z0_;
1,482✔
131
        }
132
        if (upper_token == 0 || zplane->z0_ > max_z) {
2,888✔
133
          upper_token = token;
2,888✔
134
          max_z = zplane->z0_;
2,888✔
135
        }
136
      }
137
    }
1,482✔
138

139
    // If there are no bounding z-planes, add this cell to all partitions.
140
    if (lower_token == 0) {
1,482✔
141
      for (auto& p : partitions_)
×
142
        p.push_back(i_cell);
×
143
      continue;
×
144
    }
145

146
    // Find the first partition this cell lies in.  If the lower_token indicates
147
    // a negative halfspace, then the cell is unbounded in the lower direction
148
    // and it lies in the first partition onward.  Otherwise, it is bounded by
149
    // the positive halfspace given by the lower_token.
150
    int first_partition = 0;
1,482✔
151
    if (lower_token > 0) {
1,482✔
152
      for (int i = 0; i < surfs_.size(); ++i) {
5,491✔
153
        if (lower_token == surfs_[i] + 1) {
5,491✔
154
          first_partition = i + 1;
1,444✔
155
          break;
1,444✔
156
        }
157
      }
158
    }
159

160
    // Find the last partition this cell lies in.  The logic is analogous to the
161
    // logic for first_partition.
162
    int last_partition = surfs_.size();
1,482✔
163
    if (upper_token < 0) {
1,482✔
164
      for (int i = first_partition; i < surfs_.size(); ++i) {
3,116✔
165
        if (upper_token == -(surfs_[i] + 1)) {
3,116✔
166
          last_partition = i;
1,444✔
167
          break;
1,444✔
168
        }
169
      }
170
    }
171

172
    // Add the cell to all relevant partitions.
173
    for (int i = first_partition; i <= last_partition; ++i) {
4,655✔
174
      partitions_[i].push_back(i_cell);
3,173✔
175
    }
176
  }
177
}
114✔
178

179
const vector<int32_t>& UniversePartitioner::get_cells(
316,363✔
180
  Position r, Direction u) const
181
{
182
  // Perform a binary search for the partition containing the given coordinates.
183
  int left = 0;
316,363✔
184
  int middle = (surfs_.size() - 1) / 2;
316,363✔
185
  int right = surfs_.size() - 1;
316,363✔
186
  while (true) {
187
    // Check the sense of the coordinates for the current surface.
188
    const auto& surf = *model::surfaces[surfs_[middle]];
949,122✔
189
    if (surf.sense(r, u)) {
949,122✔
190
      // The coordinates lie in the positive halfspace.  Recurse if there are
191
      // more surfaces to check.  Otherwise, return the cells on the positive
192
      // side of this surface.
193
      int right_leaf = right - (right - middle) / 2;
367,942✔
194
      if (right_leaf != middle) {
367,942✔
195
        left = middle + 1;
317,636✔
196
        middle = right_leaf;
317,636✔
197
      } else {
198
        return partitions_[middle + 1];
50,306✔
199
      }
200

201
    } else {
202
      // The coordinates lie in the negative halfspace.  Recurse if there are
203
      // more surfaces to check.  Otherwise, return the cells on the negative
204
      // side of this surface.
205
      int left_leaf = left + (middle - left) / 2;
581,180✔
206
      if (left_leaf != middle) {
581,180✔
207
        right = middle - 1;
315,123✔
208
        middle = left_leaf;
315,123✔
209
      } else {
210
        return partitions_[middle];
266,057✔
211
      }
212
    }
213
  }
632,759✔
214
}
215

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

© 2026 Coveralls, Inc