20469813151

Committed 23 Dec 2025 07:31PM UTC coverage: 82.153% (+0.01%) from 82.139%

Build # 20469813151

Build Type

Pull #3692

github

Committed by

web-flow

Commit Message

Merge 9404f8888 into a2fd6cc57

Pull Request Pull Request #3692: Fix a bug in rotational periodic boundary conditions

Run Details

17037 of 23611 branches covered (72.16%)

Branch coverage included in aggregate %.

29 of 33 new or added lines in 3 files covered. (87.88%)

157 existing lines in 4 files now uncovered.

55201 of 64320 relevant lines covered (85.82%)

43432878.37 hits per line

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

75.87

/src/tallies/filter_mesh.cpp

#include "openmc/tallies/filter_mesh.h"

#include <fmt/core.h>

#include "openmc/capi.h"
#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/mesh.h"
#include "openmc/position.h"
#include "openmc/xml_interface.h"

namespace openmc {

void MeshFilter::from_xml(pugi::xml_node node)
{
  auto bins_ = get_node_array<int32_t>(node, "bins");
  if (bins_.size() != 1) {
    fatal_error(
      "Only one mesh can be specified per " + type_str() + " mesh filter.");
  }

  auto id = bins_[0];
  auto search = model::mesh_map.find(id);
  if (search != model::mesh_map.end()) {
    set_mesh(search->second);
  } else {
    fatal_error(
      fmt::format("Could not find mesh {} specified on tally filter.", id));
  }

  if (check_for_node(node, "translation")) {
    set_translation(get_node_array<double>(node, "translation"));
  }
  // Read the rotation transform.
  if (check_for_node(node, "rotation")) {
    set_rotation(get_node_array<double>(node, "rotation"));
  }
}

void MeshFilter::get_all_bins(
  const Particle& p, TallyEstimator estimator, FilterMatch& match) const
{

  Position last_r = p.r_last();
  Position r = p.r();
  Position u = p.u();

  // apply translation if present
  if (translated_) {
    last_r -= translation();
    r -= translation();
  }
  // apply rotation if present
  if (!rotation_.empty()) {
    last_r = last_r.rotate(rotation_);
    r = r.rotate(rotation_);
    u = u.rotate(rotation_);
  }

  if (estimator != TallyEstimator::TRACKLENGTH) {
    auto bin = model::meshes[mesh_]->get_bin(r);
    if (bin >= 0) {
      match.bins_.push_back(bin);
      match.weights_.push_back(1.0);
    }
  } else {
    model::meshes[mesh_]->bins_crossed(
      last_r, r, u, match.bins_, match.weights_);
  }
}

void MeshFilter::to_statepoint(hid_t filter_group) const
{
  Filter::to_statepoint(filter_group);
  write_dataset(filter_group, "bins", model::meshes[mesh_]->id_);
  if (translated_) {
    write_dataset(filter_group, "translation", translation_);
  }
  if (rotated_) {
    write_dataset(filter_group, "rotation", rotation_);
  }
}

std::string MeshFilter::text_label(int bin) const
{
  auto& mesh = *model::meshes.at(mesh_);
  std::string label = mesh.bin_label(bin);
  return label;
}

void MeshFilter::set_mesh(int32_t mesh)
{
  // perform any additional perparation for mesh tallies here
  mesh_ = mesh;
  n_bins_ = model::meshes[mesh_]->n_bins();
  model::meshes[mesh_]->prepare_for_point_location();
}

void MeshFilter::set_translation(const Position& translation)
{
  translated_ = true;
  translation_ = translation;
}

void MeshFilter::set_translation(const double translation[3])
{
  this->set_translation({translation[0], translation[1], translation[2]});
}

void MeshFilter::set_rotation(const vector<double>& rot)
{
  rotated_ = true;

  // Compute and store the inverse rotation matrix for the angles given.
  rotation_.clear();
  rotation_.reserve(rot.size() == 9 ? 9 : 12);
  if (rot.size() == 3) {
    double phi = -rot[0] * PI / 180.0;
    double theta = -rot[1] * PI / 180.0;
    double psi = -rot[2] * PI / 180.0;
    rotation_.push_back(std::cos(theta) * std::cos(psi));
    rotation_.push_back(-std::cos(phi) * std::sin(psi) +
                        std::sin(phi) * std::sin(theta) * std::cos(psi));
    rotation_.push_back(std::sin(phi) * std::sin(psi) +
                        std::cos(phi) * std::sin(theta) * std::cos(psi));
    rotation_.push_back(std::cos(theta) * std::sin(psi));
    rotation_.push_back(std::cos(phi) * std::cos(psi) +
                        std::sin(phi) * std::sin(theta) * std::sin(psi));
    rotation_.push_back(-std::sin(phi) * std::cos(psi) +
                        std::cos(phi) * std::sin(theta) * std::sin(psi));
    rotation_.push_back(-std::sin(theta));
    rotation_.push_back(std::sin(phi) * std::cos(theta));
    rotation_.push_back(std::cos(phi) * std::cos(theta));

    // When user specifies angles, write them at end of vector
    rotation_.push_back(rot[0]);
    rotation_.push_back(rot[1]);
    rotation_.push_back(rot[2]);
  } else {
    std::copy(rot.begin(), rot.end(), std::back_inserter(rotation_));
  }
}

//==============================================================================
// C-API functions
//==============================================================================

extern "C" int openmc_mesh_filter_get_mesh(int32_t index, int32_t* index_mesh)
{
  if (!index_mesh) {
    set_errmsg("Mesh index argument is a null pointer.");
    return OPENMC_E_INVALID_ARGUMENT;
  }

  // Make sure this is a valid index to an allocated filter.
  if (int err = verify_filter(index))
    return err;

  // Get a pointer to the filter and downcast.
  const auto& filt_base = model::tally_filters[index].get();
  auto* filt = dynamic_cast<MeshFilter*>(filt_base);

  // Check the filter type.
  if (!filt) {
    set_errmsg("Tried to get mesh on a non-mesh filter.");
    return OPENMC_E_INVALID_TYPE;
  }

  // Output the mesh.
  *index_mesh = filt->mesh();
  return 0;
}

extern "C" int openmc_mesh_filter_set_mesh(int32_t index, int32_t index_mesh)
{
  // Make sure this is a valid index to an allocated filter.
  if (int err = verify_filter(index))
    return err;

  // Get a pointer to the filter and downcast.
  const auto& filt_base = model::tally_filters[index].get();
  auto* filt = dynamic_cast<MeshFilter*>(filt_base);

  // Check the filter type.
  if (!filt) {
    set_errmsg("Tried to set mesh on a non-mesh filter.");
    return OPENMC_E_INVALID_TYPE;
  }

  // Check the mesh index.
  if (index_mesh < 0 || index_mesh >= model::meshes.size()) {
    set_errmsg("Index in 'meshes' array is out of bounds.");
    return OPENMC_E_OUT_OF_BOUNDS;
  }

  // Update the filter.
  filt->set_mesh(index_mesh);
  return 0;
}

extern "C" int openmc_mesh_filter_get_translation(
  int32_t index, double translation[3])
{
  // Make sure this is a valid index to an allocated filter
  if (int err = verify_filter(index))
    return err;

  // Check the filter type
  const auto& filter = model::tally_filters[index];
  if (filter->type() != FilterType::MESH &&
      filter->type() != FilterType::MESHBORN &&
      filter->type() != FilterType::MESH_SURFACE) {
    set_errmsg("Tried to get a translation from a non-mesh-based filter.");
    return OPENMC_E_INVALID_TYPE;
  }

  // Get translation from the mesh filter and set value
  auto mesh_filter = dynamic_cast<MeshFilter*>(filter.get());
  const auto& t = mesh_filter->translation();
  for (int i = 0; i < 3; i++) {
    translation[i] = t[i];
  }

  return 0;
}

extern "C" int openmc_mesh_filter_set_translation(
  int32_t index, double translation[3])
{
  // Make sure this is a valid index to an allocated filter
  if (int err = verify_filter(index))
    return err;

  const auto& filter = model::tally_filters[index];
  // Check the filter type
  if (filter->type() != FilterType::MESH &&
      filter->type() != FilterType::MESHBORN &&
      filter->type() != FilterType::MESH_SURFACE) {
    set_errmsg("Tried to set mesh on a non-mesh-based filter.");
    return OPENMC_E_INVALID_TYPE;
  }

  // Get a pointer to the filter and downcast
  auto mesh_filter = dynamic_cast<MeshFilter*>(filter.get());

  // Set the translation
  mesh_filter->set_translation(translation);

  return 0;
}

//! Return the rotation matrix of a mesh filter
extern "C" int openmc_mesh_filter_get_rotation(
  int32_t index, double rot[], size_t* n)
{
  // Make sure this is a valid index to an allocated filter
  if (int err = verify_filter(index))
    return err;

  // Check the filter type
  const auto& filter = model::tally_filters[index];
  if (filter->type() != FilterType::MESH) {
    set_errmsg("Tried to get a rotation from a non-mesh filter.");
    return OPENMC_E_INVALID_TYPE;
  }
  // Get rotation from the mesh filter and set value
  auto mesh_filter = dynamic_cast<MeshFilter*>(filter.get());
  *n = mesh_filter->rotation().size();
  std::memcpy(rot, mesh_filter->rotation().data(),
    *n * sizeof(mesh_filter->rotation()[0]));
  return 0;
}

//! Set the flattened rotation matrix of a mesh filter
extern "C" int openmc_mesh_filter_set_rotation(
  int32_t index, const double rot[], size_t rot_len)
{
  // Make sure this is a valid index to an allocated filter
  if (int err = verify_filter(index))
    return err;

  const auto& filter = model::tally_filters[index];
  // Check the filter type
  if (filter->type() != FilterType::MESH) {
    set_errmsg("Tried to set a rotation from a non-mesh filter.");
    return OPENMC_E_INVALID_TYPE;
  }

  // Get a pointer to the filter and downcast
  auto mesh_filter = dynamic_cast<MeshFilter*>(filter.get());
  std::vector<double> vec_rot(rot, rot + rot_len);
  mesh_filter->set_rotation(vec_rot);
  return 0;
}

} // namespace openmc

1	#include "openmc/tallies/filter_mesh.h"
2
3	#include <fmt/core.h>
4
5	#include "openmc/capi.h"
6	#include "openmc/constants.h"
7	#include "openmc/error.h"
8	#include "openmc/mesh.h"
9	#include "openmc/position.h"
10	#include "openmc/xml_interface.h"
11
12	namespace openmc {
13
14	void MeshFilter::from_xml(pugi::xml_node node)	2,172✔
15	{
16	auto bins_ = get_node_array<int32_t>(node, "bins");	2,172✔
17	if (bins_.size() != 1) {	2,172!
18	fatal_error(	×
UNCOV 19	"Only one mesh can be specified per " + type_str() + " mesh filter.");	×
20	}
21
22	auto id = bins_[0];	2,172✔
23	auto search = model::mesh_map.find(id);	2,172✔
24	if (search != model::mesh_map.end()) {	2,172!
25	set_mesh(search->second);	2,172✔
26	} else {
27	fatal_error(	×
UNCOV 28	fmt::format("Could not find mesh {} specified on tally filter.", id));	×
29	}
30
31	if (check_for_node(node, "translation")) {	2,172✔
32	set_translation(get_node_array<double>(node, "translation"));	32✔
33	}
34	// Read the rotation transform.
35	if (check_for_node(node, "rotation")) {	2,172✔
36	set_rotation(get_node_array<double>(node, "rotation"));	16✔
37	}
38	}	2,172✔
39
40	void MeshFilter::get_all_bins(	443,075,473✔
41	const Particle& p, TallyEstimator estimator, FilterMatch& match) const
42	{
43
44	Position last_r = p.r_last();	443,075,473✔
45	Position r = p.r();	443,075,473✔
46	Position u = p.u();	443,075,473✔
47
48	// apply translation if present
49	if (translated_) {	443,075,473✔
50	last_r -= translation();	757,966✔
51	r -= translation();	757,966✔
52	}
53	// apply rotation if present
54	if (!rotation_.empty()) {	443,075,473✔
55	last_r = last_r.rotate(rotation_);	378,983✔
56	r = r.rotate(rotation_);	378,983✔
57	u = u.rotate(rotation_);	378,983✔
58	}
59
60	if (estimator != TallyEstimator::TRACKLENGTH) {	443,075,473✔
61	auto bin = model::meshes[mesh_]->get_bin(r);	110,088,560✔
62	if (bin >= 0) {	110,088,560✔
63	match.bins_.push_back(bin);	87,722,629✔
64	match.weights_.push_back(1.0);	87,722,629✔
65	}
66	} else {
67	model::meshes[mesh_]->bins_crossed(	332,986,913✔
68	last_r, r, u, match.bins_, match.weights_);	332,986,913✔
69	}
70	}	443,075,473✔
71
72	void MeshFilter::to_statepoint(hid_t filter_group) const	2,383✔
73	{
74	Filter::to_statepoint(filter_group);	2,383✔
75	write_dataset(filter_group, "bins", model::meshes[mesh_]->id_);	2,383✔
76	if (translated_) {	2,383✔
77	write_dataset(filter_group, "translation", translation_);	22✔
78	}
79	if (rotated_) {	2,383✔
80	write_dataset(filter_group, "rotation", rotation_);	11✔
81	}
82	}	2,383✔
83
84	std::string MeshFilter::text_label(int bin) const	5,334,479✔
85	{
86	auto& mesh = *model::meshes.at(mesh_);	5,334,479✔
87	std::string label = mesh.bin_label(bin);	5,334,479✔
88	return label;	5,334,479✔
89	}
90
91	void MeshFilter::set_mesh(int32_t mesh)	2,293✔
92	{
93	// perform any additional perparation for mesh tallies here
94	mesh_ = mesh;	2,293✔
95	n_bins_ = model::meshes[mesh_]->n_bins();	2,293✔
96	model::meshes[mesh_]->prepare_for_point_location();	2,293✔
97	}	2,293✔
98
99	void MeshFilter::set_translation(const Position& translation)	54✔
100	{
101	translated_ = true;	54✔
102	translation_ = translation;	54✔
103	}	54✔
104
105	void MeshFilter::set_translation(const double translation[3])	22✔
106	{
107	this->set_translation({translation[0], translation[1], translation[2]});	22✔
108	}	22✔
109
110	void MeshFilter::set_rotation(const vector<double>& rot)	27✔
111	{
112	rotated_ = true;	27✔
113
114	// Compute and store the inverse rotation matrix for the angles given.
115	rotation_.clear();	27✔
116	rotation_.reserve(rot.size() == 9 ? 9 : 12);	27!
117	if (rot.size() == 3) {	27!
118	double phi = -rot[0] * PI / 180.0;	27✔
119	double theta = -rot[1] * PI / 180.0;	27✔
120	double psi = -rot[2] * PI / 180.0;	27✔
121	rotation_.push_back(std::cos(theta) * std::cos(psi));	27✔
UNCOV 122	rotation_.push_back(-std::cos(phi) * std::sin(psi) +	×
123	std::sin(phi) * std::sin(theta) * std::cos(psi));	27✔
UNCOV 124	rotation_.push_back(std::sin(phi) * std::sin(psi) +	×
125	std::cos(phi) * std::sin(theta) * std::cos(psi));	27✔
126	rotation_.push_back(std::cos(theta) * std::sin(psi));	27✔
UNCOV 127	rotation_.push_back(std::cos(phi) * std::cos(psi) +	×
128	std::sin(phi) * std::sin(theta) * std::sin(psi));	27✔
UNCOV 129	rotation_.push_back(-std::sin(phi) * std::cos(psi) +	×
130	std::cos(phi) * std::sin(theta) * std::sin(psi));	27✔
131	rotation_.push_back(-std::sin(theta));	27✔
132	rotation_.push_back(std::sin(phi) * std::cos(theta));	27✔
133	rotation_.push_back(std::cos(phi) * std::cos(theta));	27✔
134
135	// When user specifies angles, write them at end of vector
136	rotation_.push_back(rot[0]);	27✔
137	rotation_.push_back(rot[1]);	27✔
138	rotation_.push_back(rot[2]);	27✔
139	} else {
UNCOV 140	std::copy(rot.begin(), rot.end(), std::back_inserter(rotation_));	×
141	}
142	}	27✔
143
144	//==============================================================================
145	// C-API functions
146	//==============================================================================
147
148	extern "C" int openmc_mesh_filter_get_mesh(int32_t index, int32_t* index_mesh)	550✔
149	{
150	if (!index_mesh) {	550!
UNCOV 151	set_errmsg("Mesh index argument is a null pointer.");	×
UNCOV 152	return OPENMC_E_INVALID_ARGUMENT;	×
153	}
154
155	// Make sure this is a valid index to an allocated filter.
156	if (int err = verify_filter(index))	550!
UNCOV 157	return err;	×
158
159	// Get a pointer to the filter and downcast.
160	const auto& filt_base = model::tally_filters[index].get();	550✔
161	auto* filt = dynamic_cast<MeshFilter*>(filt_base);	550!
162
163	// Check the filter type.
164	if (!filt) {	550!
165	set_errmsg("Tried to get mesh on a non-mesh filter.");	×
166	return OPENMC_E_INVALID_TYPE;	×
167	}
168
169	// Output the mesh.
170	*index_mesh = filt->mesh();	550✔
171	return 0;	550✔
172	}
173
174	extern "C" int openmc_mesh_filter_set_mesh(int32_t index, int32_t index_mesh)	500✔
175	{
176	// Make sure this is a valid index to an allocated filter.
177	if (int err = verify_filter(index))	500!
UNCOV 178	return err;	×
179
180	// Get a pointer to the filter and downcast.
181	const auto& filt_base = model::tally_filters[index].get();	500✔
182	auto* filt = dynamic_cast<MeshFilter*>(filt_base);	500!
183
184	// Check the filter type.
185	if (!filt) {	500!
UNCOV 186	set_errmsg("Tried to set mesh on a non-mesh filter.");	×
UNCOV 187	return OPENMC_E_INVALID_TYPE;	×
188	}
189
190	// Check the mesh index.
191	if (index_mesh < 0 \|\| index_mesh >= model::meshes.size()) {	500!
192	set_errmsg("Index in 'meshes' array is out of bounds.");	×
UNCOV 193	return OPENMC_E_OUT_OF_BOUNDS;	×
194	}
195
196	// Update the filter.
197	filt->set_mesh(index_mesh);	500✔
198	return 0;	500✔
199	}
200
201	extern "C" int openmc_mesh_filter_get_translation(	22✔
202	int32_t index, double translation[3])
203	{
204	// Make sure this is a valid index to an allocated filter
205	if (int err = verify_filter(index))	22!
UNCOV 206	return err;	×
207
208	// Check the filter type
209	const auto& filter = model::tally_filters[index];	22✔
210	if (filter->type() != FilterType::MESH &&	22✔
211	filter->type() != FilterType::MESHBORN &&	33!
212	filter->type() != FilterType::MESH_SURFACE) {	11!
UNCOV 213	set_errmsg("Tried to get a translation from a non-mesh-based filter.");	×
UNCOV 214	return OPENMC_E_INVALID_TYPE;	×
215	}
216
217	// Get translation from the mesh filter and set value
218	auto mesh_filter = dynamic_cast<MeshFilter*>(filter.get());	22!
219	const auto& t = mesh_filter->translation();	22✔
220	for (int i = 0; i < 3; i++) {	88✔
221	translation[i] = t[i];	66✔
222	}
223
224	return 0;	22✔
225	}
226
227	extern "C" int openmc_mesh_filter_set_translation(	22✔
228	int32_t index, double translation[3])
229	{
230	// Make sure this is a valid index to an allocated filter
231	if (int err = verify_filter(index))	22!
UNCOV 232	return err;	×
233
234	const auto& filter = model::tally_filters[index];	22✔
235	// Check the filter type
236	if (filter->type() != FilterType::MESH &&	22✔
237	filter->type() != FilterType::MESHBORN &&	33!
238	filter->type() != FilterType::MESH_SURFACE) {	11!
UNCOV 239	set_errmsg("Tried to set mesh on a non-mesh-based filter.");	×
UNCOV 240	return OPENMC_E_INVALID_TYPE;	×
241	}
242
243	// Get a pointer to the filter and downcast
244	auto mesh_filter = dynamic_cast<MeshFilter*>(filter.get());	22!
245
246	// Set the translation
247	mesh_filter->set_translation(translation);	22✔
248
249	return 0;	22✔
250	}
251
252	//! Return the rotation matrix of a mesh filter
253	extern "C" int openmc_mesh_filter_get_rotation(	11✔
254	int32_t index, double rot[], size_t* n)
255	{
256	// Make sure this is a valid index to an allocated filter
257	if (int err = verify_filter(index))	11!
UNCOV 258	return err;	×
259
260	// Check the filter type
261	const auto& filter = model::tally_filters[index];	11✔
262	if (filter->type() != FilterType::MESH) {	11!
UNCOV 263	set_errmsg("Tried to get a rotation from a non-mesh filter.");	×
UNCOV 264	return OPENMC_E_INVALID_TYPE;	×
265	}
266	// Get rotation from the mesh filter and set value
267	auto mesh_filter = dynamic_cast<MeshFilter*>(filter.get());	11!
268	*n = mesh_filter->rotation().size();	11✔
269	std::memcpy(rot, mesh_filter->rotation().data(),	11✔
270	n sizeof(mesh_filter->rotation()[0]));	11✔
271	return 0;	11✔
272	}
273
274	//! Set the flattened rotation matrix of a mesh filter
275	extern "C" int openmc_mesh_filter_set_rotation(	11✔
276	int32_t index, const double rot[], size_t rot_len)
277	{
278	// Make sure this is a valid index to an allocated filter
279	if (int err = verify_filter(index))	11!
UNCOV 280	return err;	×
281
282	const auto& filter = model::tally_filters[index];	11✔
283	// Check the filter type
284	if (filter->type() != FilterType::MESH) {	11!
UNCOV 285	set_errmsg("Tried to set a rotation from a non-mesh filter.");	×
UNCOV 286	return OPENMC_E_INVALID_TYPE;	×
287	}
288
289	// Get a pointer to the filter and downcast
290	auto mesh_filter = dynamic_cast<MeshFilter*>(filter.get());	11!
291	std::vector<double> vec_rot(rot, rot + rot_len);	11✔
292	mesh_filter->set_rotation(vec_rot);	11✔
293	return 0;	11✔
294	}	11✔
295
296	} // namespace openmc

openmc-dev / openmc / 20469813151

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

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